Enable plots of ElementalNodal data (#511)

Author: PProfizi

examples/06-plotting/00-basic_plotting.py

@@ -52,9 +52,8 @@
 # - The "off_screen" keyword only works when "notebook=False" to prevent the GUI from appearing.
 
 
-# Plot a field on its supporting mesh (field location must be Elemental or Nodal)
+# Plot a field on its supporting mesh
 stress = model.results.stress()
-stress.inputs.requested_location.connect(dpf.locations.nodal)
 fc = stress.outputs.fields_container()
 field = fc[0]
 field.plot(notebook=False, shell_layers=None, show_axes=True, title="Field", text="Field plot")

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

@@ -33,6 +33,8 @@
 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
 from ansys.dpf.core.check_version import meets_version, server_meet_version, version_requires
@@ -713,3 +715,48 @@ 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: locations) -> int:
+        """Return the data length for a given mesh location.
+
+        Accepted mesh locations are nodal, elemental, faces, and elemental_nodal.
+
+        Parameters
+        ----------
+        location:
+            The mesh location to compute data length for.
+            Can be nodal, elemental, faces, or elemental_nodal.
+
+        Returns
+        -------
+        data_size:
+            If location is nodal, return the number of nodes.
+            If location is elemental, return the number of elements.
+            If location is faces, return the number of faces.
+            If location is elemental nodal, return the sum of the number of nodes per element.
+        """
+        if location == locations.nodal:
+            return len(self.nodes)
+        elif location == locations.elemental:
+            return len(self.elements)
+        elif location == locations.faces:
+            return len(self.faces)
+        elif location == locations.elemental_nodal:
+            return np.sum(self.get_elemental_nodal_size_list())
+        elif location == locations.overall:
+            return len(self.elements)
+        else:
+            raise TypeError(f"Location {location} is not recognized.")
+
+    def get_elemental_nodal_size_list(self) -> np.ndarray:
+        """Return the array of number of nodes per element in the mesh."""
+        # 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

@@ -147,6 +147,7 @@ def add_mesh(self, meshed_region, deform_by=None, scale_factor=1.0, as_linear=Tr
                 grid = meshed_region._as_vtk(
                     meshed_region.nodes.coordinates_field, as_linear=as_linear
                 )
+                meshed_region._full_grid = grid
                 meshed_region.as_linear = as_linear
             else:
                 grid = meshed_region.grid
@@ -317,8 +318,14 @@ 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, faces, or overall location are supported for plotting."
+            )
 
         # Treat multilayered shells
         if not isinstance(shell_layer, eshell_layers):
@@ -333,13 +340,27 @@ 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])]
+                )
+                mask = mask_2
+                ind = ind_2
             overall_data[ind] = field.data[mask]
         else:
             overall_data[:] = field.data[0]
@@ -348,12 +369,22 @@ def add_field(
         # Have to remove any active scalar field from the pre-existing grid object,
         # otherwise we get two scalar bars when calling several plot_contour on the same mesh
         # but not for the same field. The PyVista UnstructuredGrid keeps memory of it.
-        if not deform_by:
-            grid = meshed_region.grid
-        else:
+        if location == locations.elemental_nodal:
+            as_linear = False
+        if deform_by:
             grid = meshed_region._as_vtk(
-                meshed_region.deform_by(deform_by, scale_factor), as_linear
+                meshed_region.deform_by(deform_by, scale_factor), as_linear=as_linear
             )
+        else:
+            if as_linear != meshed_region.as_linear:
+                grid = meshed_region._as_vtk(
+                    meshed_region.nodes.coordinates_field, as_linear=as_linear
+                )
+                meshed_region.as_linear = as_linear
+            else:
+                grid = meshed_region.grid
+        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)
 
@@ -976,15 +1007,14 @@ def plot_contour(
             import warnings
 
             warnings.simplefilter("ignore")
-
         if isinstance(field_or_fields_container, (dpf.core.Field, dpf.core.FieldsContainer)):
             fields_container = None
             if isinstance(field_or_fields_container, dpf.core.Field):
                 fields_container = dpf.core.FieldsContainer(
                     server=field_or_fields_container._server
                 )
-                fields_container.add_label(DefinitionLabels.time)
-                fields_container.add_field({DefinitionLabels.time: 1}, field_or_fields_container)
+                fields_container.add_label("id")
+                fields_container.add_field({"id": 1}, field_or_fields_container)
             elif isinstance(field_or_fields_container, dpf.core.FieldsContainer):
                 fields_container = field_or_fields_container
         else:
@@ -1022,43 +1052,96 @@ 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)
         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")
-        for field in fields_container:
-            shell_layer_check = field.shell_layers
-            if shell_layer_check in [
-                eshell_layers.topbottom,
-                eshell_layers.topbottommid,
-            ]:
-                changeOp.inputs.fields_container.connect(fields_container)
-                sl = eshell_layers.top
-                if shell_layers is not None:
-                    if not isinstance(shell_layers, eshell_layers):
-                        raise TypeError(
-                            "shell_layer attribute must be a core.shell_layers instance."
-                        )
-                    sl = shell_layers
-                changeOp.inputs.e_shell_layer.connect(sl.value)  # top layers taken
-                fields_container = changeOp.get_output(0, core.types.fields_container)
-                break
+        changeOp = core.operators.utility.change_shell_layers()
+        if location == locations.elemental_nodal:
+            # change_shell_layers does not support elemental_nodal when given a fields_container
+            new_fields_container = dpf.core.FieldsContainer()
+            for l in fields_container.labels:
+                new_fields_container.add_label(l)
+            for i, field in enumerate(fields_container):
+                label_space_i = fields_container.get_label_space(i)
+                shell_layer_check = field.shell_layers
+                if shell_layer_check in [
+                    eshell_layers.topbottom,
+                    eshell_layers.topbottommid,
+                ]:
+                    changeOp.inputs.fields_container.connect(field)
+                    changeOp.inputs.merge.connect(True)
+                    sl = eshell_layers.top
+                    if shell_layers is not None:
+                        if not isinstance(shell_layers, eshell_layers):
+                            raise TypeError(
+                                "shell_layer attribute must be a core.shell_layers instance."
+                            )
+                        sl = shell_layers
+                    changeOp.inputs.e_shell_layer.connect(sl.value)  # top layers taken
+                    field = changeOp.get_output(0, core.types.field)
+                new_fields_container.add_field(label_space=label_space_i, field=field)
+            fields_container = new_fields_container
+        else:
+            for field in fields_container:
+                shell_layer_check = field.shell_layers
+                if shell_layer_check in [
+                    eshell_layers.topbottom,
+                    eshell_layers.topbottommid,
+                ]:
+                    changeOp.inputs.fields_container.connect(fields_container)
+                    sl = eshell_layers.top
+                    if shell_layers is not None:
+                        if not isinstance(shell_layers, eshell_layers):
+                            raise TypeError(
+                                "shell_layer attribute must be a core.shell_layers instance."
+                            )
+                        sl = shell_layers
+                    changeOp.inputs.e_shell_layer.connect(sl.value)  # top layers taken
+                    fields_container = changeOp.get_output(0, core.types.fields_container)
+                    break
 
         # 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)
+            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)
+                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])]
+                )
+                mask = mask_2
+                ind = ind_2
             overall_data[ind] = field.data[mask]
 
         # create the plotter and add the meshes
@@ -1087,15 +1170,20 @@ def plot_contour(
             bound_method=self._internal_plotter._plotter.add_mesh, **kwargs
         )
         as_linear = True
+        if location == locations.elemental_nodal:
+            as_linear = False
         if deform_by:
             grid = mesh._as_vtk(mesh.deform_by(deform_by, scale_factor), as_linear=as_linear)
             self._internal_plotter.add_scale_factor_legend(scale_factor, **kwargs)
         else:
             if as_linear != mesh.as_linear:
                 grid = mesh._as_vtk(mesh.nodes.coordinates_field, as_linear=as_linear)
+                mesh._full_grid = grid
                 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)