Enable plots of ElementalNodal data

src/ansys/dpf/core/elements.py

@@ -677,7 +677,7 @@ def map_scoping(self, external_scope):
 
         Examples
         --------
-        Return the indices that map a field to an elements collection.
+        Return the indices that map a field to an Elements collection.
 
         >>> import ansys.dpf.core as dpf
         >>> from ansys.dpf.core import examples
@@ -689,7 +689,7 @@ def map_scoping(self, external_scope):
 
         """
         if external_scope.location in ["Nodal", "NodalElemental"]:
-            raise ValueError('Input scope location must be "Nodal"')
+            raise ValueError('Input scope location must be "Elemental"')
         arr = np.array(list(map(self.mapping_id_to_index.get, external_scope.ids)))
         mask = arr != None
         ind = arr[mask].astype(np.int32)

src/ansys/dpf/core/meshed_region.py

@@ -32,6 +32,7 @@
 
 import traceback
 import warnings
+import numpy as np
 
 from ansys.dpf.core import field, property_field, scoping, server as server_module
 from ansys.dpf.core.cache import class_handling_cache
@@ -713,3 +714,35 @@ def is_empty(self) -> bool:
         no_elements = self.elements.n_elements == 0
         no_nodes = self.nodes.n_nodes == 0
         return no_nodes and no_faces and no_elements
+
+    def location_data_len(self, location):
+        """
+
+        Parameters
+        ----------
+        location
+
+        Returns
+        -------
+        data_size : int
+        """
+        if location == locations.nodal:
+            return len(self.nodes)
+        elif location == locations.elemental:
+            return len(self.elements)
+        elif location == locations.elemental_nodal:
+            return np.sum(self.get_elemental_nodal_size_list())
+        else:
+            raise TypeError(f"Location {location} is not recognized.")
+
+    def get_elemental_nodal_size_list(self):
+        # Get the field of element types
+        element_types_field = self.elements.element_types_field
+        # get the number of nodes for each possible element type
+        size_map = dict([(e_type.value, element_types.descriptor(e_type).n_nodes)
+                         for e_type in element_types])
+        keys = list(size_map.keys())
+        sort_idx = np.argsort(keys)
+        idx = np.searchsorted(keys, element_types_field.data, sorter=sort_idx)
+        return np.asarray(list(size_map.values()))[sort_idx][idx]
+

src/ansys/dpf/core/plotter.py

@@ -317,8 +317,16 @@ def add_field(
                 show_min = False
         elif location == locations.overall:
             mesh_location = meshed_region.elements
+        elif location == locations.elemental_nodal:
+            mesh_location = meshed_region.elements
+            # If ElementalNodal, first extend results to mid-nodes
+            field = dpf.core.operators.averaging.extend_to_mid_nodes(
+                field=field
+            ).eval()
         else:
-            raise ValueError("Only elemental, nodal or faces location are supported for plotting.")
+            raise ValueError(
+                "Only elemental, elemental nodal, nodal or faces location are supported for plotting."
+            )
 
         # Treat multilayered shells
         if not isinstance(shell_layer, eshell_layers):
@@ -333,13 +341,25 @@ def add_field(
             )
             field = change_shell_layer_op.get_output(0, core.types.field)
 
+        location_data_len = meshed_region.location_data_len(location)
         component_count = field.component_count
         if component_count > 1:
-            overall_data = np.full((len(mesh_location), component_count), np.nan)
+            overall_data = np.full((location_data_len, component_count), np.nan)
         else:
-            overall_data = np.full(len(mesh_location), np.nan)
+            overall_data = np.full(location_data_len, np.nan)
         if location != locations.overall:
             ind, mask = mesh_location.map_scoping(field.scoping)
+
+            # Rework ind and mask to take into account n_nodes per element if ElementalNodal
+            if location == locations.elemental_nodal:
+                n_nodes_list = meshed_region.get_elemental_nodal_size_list().astype(np.int32)
+                first_index = np.insert(np.cumsum(n_nodes_list)[:-1], 0, 0).astype(np.int32)
+                mask_2 = np.asarray([mask_i for i, mask_i in enumerate(mask)
+                                     for _ in range(n_nodes_list[ind[i]])])
+                ind_2 = np.asarray([first_index[ind_i]+j for ind_i in ind
+                                    for j in range(n_nodes_list[ind_i])])  # OK
+                mask = mask_2
+                ind = ind_2
             overall_data[ind] = field.data[mask]
         else:
             overall_data[:] = field.data[0]
@@ -354,6 +374,8 @@ def add_field(
             grid = meshed_region._as_vtk(
                 meshed_region.deform_by(deform_by, scale_factor), as_linear
             )
+        if location == locations.elemental_nodal:
+            grid = grid.shrink(1.0)
         grid.set_active_scalars(None)
         self._plotter.add_mesh(grid, scalars=overall_data, **kwargs_in)
 
@@ -1022,14 +1044,23 @@ def plot_contour(
                 unit = field.unit
                 break
 
+        # If ElementalNodal, first extend results to mid-nodes
+        if location == locations.elemental_nodal:
+            fields_container = dpf.core.operators.averaging.extend_to_mid_nodes_fc(
+                fields_container=fields_container
+            ).eval()
+
+        location_data_len = mesh.location_data_len(location)  # 3751 nodes  # 3000 elements  # 24000 elemental nodal  # field.data 72000
         if location == locations.nodal:
             mesh_location = mesh.nodes
         elif location == locations.elemental:
             mesh_location = mesh.elements
         elif location == locations.faces:
             mesh_location = mesh.faces
+        elif location == locations.elemental_nodal:
+            mesh_location = mesh.elements
         else:
-            raise ValueError("Only elemental, nodal or faces location are supported for plotting.")
+            raise ValueError("Only elemental, elemental nodal, nodal or faces location are supported for plotting.")
 
         # pre-loop: check if shell layers for each field, if yes, set the shell layers
         changeOp = core.Operator("change_shellLayers")
@@ -1039,6 +1070,8 @@ def plot_contour(
                 eshell_layers.topbottom,
                 eshell_layers.topbottommid,
             ]:
+                if location == locations.elemental_nodal:
+                    raise TypeError("Trying to plot ElementalNodal values for shells.")
                 changeOp.inputs.fields_container.connect(fields_container)
                 sl = eshell_layers.top
                 if shell_layers is not None:
@@ -1053,13 +1086,27 @@ def plot_contour(
 
         # Merge field data into a single array
         if component_count > 1:
-            overall_data = np.full((len(mesh_location), component_count), np.nan)
+            overall_data = np.full((location_data_len, component_count), np.nan)
         else:
-            overall_data = np.full(len(mesh_location), np.nan)
+            overall_data = np.full(location_data_len, np.nan)
+
+        # field._data_pointer gives the first index of each entity data
+        # (should be of size nb_elements)
 
         for field in fields_container:
             ind, mask = mesh_location.map_scoping(field.scoping)
-            overall_data[ind] = field.data[mask]
+            if location == locations.elemental_nodal:
+                # Rework ind and mask to take into account n_nodes per element
+                # entity_index_map = field._data_pointer
+                n_nodes_list = mesh.get_elemental_nodal_size_list().astype(np.int32)  # OK
+                first_index = np.insert(np.cumsum(n_nodes_list)[:-1], 0, 0).astype(np.int32)
+                mask_2 = np.asarray([mask_i for i, mask_i in enumerate(mask)
+                                     for _ in range(n_nodes_list[ind[i]])])
+                ind_2 = np.asarray([first_index[ind_i]+j for ind_i in ind
+                                    for j in range(n_nodes_list[ind_i])])  # OK
+                mask = mask_2
+                ind = ind_2
+            overall_data[ind] = field.data[mask]  # (24000,3)[:3000] = (24000,3)[:3000]
 
         # create the plotter and add the meshes
 
@@ -1096,6 +1143,8 @@ def plot_contour(
                 mesh.as_linear = as_linear
             else:
                 grid = mesh.grid
+        if location == locations.elemental_nodal:
+            grid = grid.shrink(1.0)
         grid.clear_data()
         self._internal_plotter._plotter.add_mesh(grid, scalars=overall_data, **kwargs_in)
 

tests/test_plotter.py

@@ -129,7 +129,7 @@ def test_plotter_on_fields_container_elemental(allkindofcomplexity):
     avg_op.inputs.fields_container.connect(stress.outputs.fields_container)
     fc = avg_op.outputs.fields_container()
     pl = Plotter(model.metadata.meshed_region)
-    cpos = pl.plot_contour(fc)
+    _ = pl.plot_contour(fc)
 
 
 @pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista")
@@ -141,7 +141,7 @@ def test_plotter_on_fields_container_nodal(allkindofcomplexity):
     avg_op.inputs.fields_container.connect(stress.outputs.fields_container)
     fc = avg_op.outputs.fields_container()
     pl = Plotter(model.metadata.meshed_region)
-    cpos = pl.plot_contour(fc)
+    _ = pl.plot_contour(fc)
 
 
 @pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista")
@@ -199,6 +199,53 @@ def test_field_nodal_plot(allkindofcomplexity):
     remove_picture(picture)
 
 
+@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista")
+def test_field_elemental_nodal_plot_simple(simple_bar):
+    model = Model(simple_bar)
+    stress = model.results.element_nodal_forces()
+    fc = stress.outputs.fields_container()
+    f = fc[0]
+    print(f.data.shape)
+    f.plot()
+    picture = 'test_plotter1.png'
+    remove_picture(picture)
+    f.plot(off_screen=True, screenshot=picture)
+    assert os.path.exists(os.path.join(os.getcwd(), picture))
+    remove_picture(picture)
+
+
+@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista")
+def test_field_elemental_nodal_plot_scoped(simple_bar):
+    model = Model(simple_bar)
+    mesh_scoping = dpf.core.mesh_scoping_factory.elemental_scoping(element_ids=list(range(1501,
+                                                                                          3001)))
+    stress = model.results.element_nodal_forces.on_mesh_scoping(mesh_scoping)
+    fc = stress.eval()
+    f = fc[0]
+    print(f.data.shape)
+    f.plot()
+    picture = 'test_plotter2.png'
+    remove_picture(picture)
+    f.plot(off_screen=True, screenshot=picture)
+    assert os.path.exists(os.path.join(os.getcwd(), picture))
+    remove_picture(picture)
+
+
+@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista")
+def test_field_elemental_nodal_plot_multiple_solid_types():
+    from ansys.dpf.core import examples
+    model = dpf.core.Model(examples.download_hemisphere())
+    stress = model.results.stress()
+    fc = stress.outputs.fields_container()
+    f = fc[0]
+    f.plot()
+    picture = 'test_plotter3.png'
+    remove_picture(picture)
+    f.plot(off_screen=True, screenshot=picture)
+    assert os.path.exists(os.path.join(os.getcwd(), picture))
+    remove_picture(picture)
+
+
 @pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista")
 def test_field_solid_plot(allkindofcomplexity):
     model = Model(allkindofcomplexity)