Add quadratic and cubic shape functions

.github/workflows/python-app.yml

@@ -19,10 +19,10 @@ jobs:
 
     steps:
     - uses: actions/checkout@v4
-    - name: Set up Python 3.10
+    - name: Set up Python 3.12.4
       uses: actions/setup-python@v3
       with:
-        python-version: "3.10"
+        python-version: "3.12.4"
     - name: Install dependencies
       run: |
         sudo apt-get update

plutho/mesh/custom_mesh_parser.py

@@ -10,6 +10,10 @@
 import numpy as np
 import numpy.typing as npt
 
+# Local libraries
+from .gmsh_parser import element_to_nodes_map, \
+    element_to_boundary_nodes_map
+
 
 class TokenType(Enum):
     SectionStart = "SectionStart",
@@ -140,10 +144,14 @@ class CustomParser:
     nodes: npt.NDArray
     elements: List[Element]
 
-    def __init__(self, file_path: str):
+    element_order: int
+
+    def __init__(self, file_path: str, element_order: int):
         if not os.path.isfile(file_path):
             raise IOError(f"Mesh file {file_path} not found.")
 
+        self.element_order = element_order
+
         # Load text
         mesh_text = ""
         with open(file_path, "r", encoding="UTF-8") as fd:
@@ -167,9 +175,10 @@ def get_mesh_nodes_and_elements(self) -> Tuple[npt.NDArray, npt.NDArray]:
             if element.type == 2:
                 triangle_elements.append(element)
 
+        nodes_per_element = element_to_nodes_map[2, self.element_order]
         number_of_triangle_elements = len(triangle_elements)
         elements = np.zeros(
-            shape=(number_of_triangle_elements, 3),
+            shape=(number_of_triangle_elements, nodes_per_element),
             dtype=np.int64
         )
 
@@ -235,11 +244,13 @@ def get_triangle_elements_by_physical_group(
         elements_pg = self.get_elements_by_physical_group(physical_group_name)
         _, elements = self.get_mesh_nodes_and_elements()
 
+        nodes_per_element = element_to_nodes_map[2, self.element_order]
+
         # Find triangle elements, which share the same nodes
         triangle_elements = []
         for element_pg in elements_pg:
             for element in elements:
-                if len(element) != 3:
+                if len(element) != nodes_per_element:
                     continue
 
                 found = True
@@ -391,10 +402,12 @@ def _parse_elements(self):
                 tags.append(int(token.value))
 
             # Parse nodes
-            # element_type = 1 -> Line -> 2 Points
-            # element_type = 2 -> Triangle -> 3 Points
+            number_of_nodes = element_to_nodes_map[
+                element_type,
+                self.element_order
+            ]
             nodes = []
-            for _ in range(element_type+1):
+            for _ in range(number_of_nodes):
                 token = self._next_token()
                 expect_token(token, TokenType.TokenInt)
                 # Subtract 1 because in gmsh the indices start at 1

plutho/mesh/gmsh_parser.py

@@ -9,11 +9,28 @@
 import numpy as np
 import numpy.typing as npt
 
+
+# First index represents the element type:
+# 1: Line
+# 2: Triangle
+# Second index represents the element order
+element_to_nodes_map = np.array([
+    [-1, -1, -1, -1],
+    [-1, 2, 3, 4],
+    [-1, 3, 6, 10]
+])
+
+# Since its the same for lines and triangles, this map only has one row
+# and the index represents the element order
+element_to_boundary_nodes_map = np.array([-1, 2, 3, 4])
+
+
 class GmshParser:
     """Class to use the gmsh python interface to read the mesh files.
     """
+    element_order: int
 
-    def __init__(self, file_path):
+    def __init__(self, file_path, element_order):
         if not gmsh.is_initialized():
             gmsh.initialize()
 
@@ -22,6 +39,8 @@ def __init__(self, file_path):
         else:
             raise IOError(f"Mesh file {file_path} not found.")
 
+        self.element_order = element_order
+
     @staticmethod
     def _get_nodes(node_tags, node_coords, _=None) -> npt.NDArray:
         """Internal function to return a list of the node coordinates based
@@ -43,44 +62,106 @@ def _get_nodes(node_tags, node_coords, _=None) -> npt.NDArray:
 
         return nodes
 
-    @staticmethod
     def _get_elements(
+        self,
         element_types,
         element_tags,
         element_node_tags
-    ) -> npt.NDArray:
-        """Returns a list of elements used in the simulation from the gmsh
-        getElements() api function.
+    ) -> Dict[int, npt.NDArray]:
+        """
+        For each given element type a list of elements is returned.
 
         Parameters:
-            element_types: Element types from gmsh api. Similar to the element
-                dimension so 2 equals triangles.
-            element_tags: Tags of the element per element type.
-            element_node_tags: Tags of the nodes per element tag.
+            element_types: List of element types for which the elements lists
+                are extracted.
+            element_tags: Tags of the elements per element type.
+            element_node_tags: Tags of the nodes of the elements per element
+                type.
 
         Returns:
-            A list of lists where each inner list contains the node indices for
-            the triangle at the outer list index.
+            A list of elements for each given element type:
+            List of elements [e1, e2, e3, ...] where each element consists of
+            its node tags: e1: [n1, n2, ...].
         """
-        elements = np.zeros(shape=(len(element_tags[1]), 3), dtype=int)
+        elements_per_type = {}
+
+        # Types of the given elements: 2D -> Triangle or 1D -> Line
+        # Iterate over all given types
         for i, element_type in enumerate(element_types):
-            if element_type == 2:
-                # Only looking for 3-node triangle elements
-                current_element_tags = element_tags[i]
-                current_node_tags = element_node_tags[i]
-                for j, _ in enumerate(current_element_tags):
-                    # 1 is subtracted because the indices in gmsh start with 1.
-                    elements[j] = current_node_tags[3*j:3*(j+1)] - np.ones(3)
+            elements = []
+            element_indices = []
+
+            # For each type get the order to calculate the number of nodes
+            # it consists of
+            _, dim, order, _, _, _ = \
+                gmsh.model.mesh.getElementProperties(
+                    element_type
+                )
 
-        return elements
+            match dim:
+                case 0:
+                    nodes_per_element = 1
+                case 1:
+                    # Line element
+                    nodes_per_element = order+1
+                case 2:
+                    # Triangle element
+                    nodes_per_element = int(1/2*(order+1)*(order+2))
+                case _:
+                    raise ValueError("Only supporting 2D meshes")
+
+            # The element tags and element_node_tags lists are nested lists
+            # which contain the element tags and node tags of the elements of
+            # the current type at the respective index
+            current_element_tags = element_tags[i]
+            current_node_tags = element_node_tags[i]
+
+            # Now iterate over every element of the current type and extract
+            # the node indices
+            for j, _ in enumerate(current_element_tags):
+                # 1 is subtracted because the indices in gmsh start with 1.
+                elements.append(current_node_tags[
+                        nodes_per_element*j:nodes_per_element*(j+1)
+                    ] - np.ones(nodes_per_element)
+                )
+                element_indices.append(j)
+
+            if len(elements) == 0:
+                raise ValueError(
+                    "Couldn't return elements because the list is empty."
+                )
+
+            elements_np = np.zeros(
+                shape=(len(current_element_tags), nodes_per_element),
+                dtype=int
+            )
+            for index, element in zip(element_indices, elements):
+                elements_np[index] = element
+
+            elements_per_type[element_type] = elements_np
+
+        return elements_per_type
 
     def get_mesh_nodes_and_elements(self) -> Tuple[npt.NDArray, npt.NDArray]:
         """Creates the nodes and elements lists as used in the simulation.
 
         Returns:
             List of nodes and elements"""
         nodes = self._get_nodes(*gmsh.model.mesh.getNodes())
-        elements = self._get_elements(*gmsh.model.mesh.getElements())
+
+        el_types, el_tags, el_node_tags = gmsh.model.mesh.getElements(dim=2)
+        elements_per_type = self._get_elements(
+            el_types, el_tags, el_node_tags
+        )
+
+        if len(elements_per_type) != 1:
+            raise ValueError(
+                "The given mesh as more (or less) than one element type for "
+                "dim=2"
+            )
+
+        # There should be only one element type for dim=2
+        elements = elements_per_type[el_types[0]]
 
         return nodes, elements
 
@@ -122,41 +203,42 @@ def get_elements_by_physical_groups(
         self,
         needed_pg_names: List[str]
     ) -> Dict[str, npt.NDArray]:
-        """Returns the elements inside the given physical groups.
+        """Get all triangle elements for the given physical group.
 
         Parameters:
             needed_pg_names: List of names of physical groups for which the
-                triangles are returned.
+                triangle elements shall be returned.
 
         Returns:
-            Dictionary where keys are the pg names and the values are a list
-            of triangles of this physical group."""
-        pg_tags = self.get_nodes_by_physical_groups(needed_pg_names)
-        elements = self._get_elements(*gmsh.model.mesh.getElements())
-        triangle_elements = {}
-        for pg_name, nodes in pg_tags.items():
-            current_triangle_elements = []
-            for check_element in elements:
-                # If at least 2 nodes of the check_element are inside of the
-                # nodes list of the current physical group, then the element
-                # is also part of the physical group.
-                found_count = 0
-
-                if check_element[0] in nodes:
-                    found_count += 1
-                if check_element[1] in nodes:
-                    found_count += 1
-                if check_element[2] in nodes:
-                    found_count += 1
-
-                if found_count > 1:
-                    current_triangle_elements.append(check_element)
-
-            triangle_elements[pg_name] = np.array(
-                current_triangle_elements,
-                dtype=int)
-
-        return triangle_elements
+            dictionary where the key is the physical group name and the value
+            is a list of elements for this pg.
+        """
+        _, elements = self.get_mesh_nodes_and_elements()
+        pg_nodes = self.get_nodes_by_physical_groups(needed_pg_names)
+        element_order = self.element_order
+        nodes_per_line = element_order+1
+        pg_elements = {}
+
+        for pg_name in needed_pg_names:
+            nodes = pg_nodes[pg_name]
+
+            # For every possible element, check if it contains 'enought' nodes
+            # from the physical group -> Then add it to the elements list
+            # TODO Rework this: Check if its also possible to return
+            # element lists containing line elements.
+            current_elements = []
+            for element in elements:
+                count = 0
+                for node_index in element:
+                    if node_index in nodes:
+                        count += 1
+
+                if count >= nodes_per_line:
+                    current_elements.append(element)
+
+            pg_elements[pg_name] = np.array(current_elements)
+
+        return pg_elements
 
     def create_element_post_processing_view(
         self,

plutho/mesh/mesh.py

@@ -6,6 +6,7 @@
 
 # Third party libraries
 import gmsh
+import numpy as np
 import numpy.typing as npt
 
 # Local libraries
@@ -22,15 +23,18 @@ class Mesh:
     mesh_file_path: str
     file_version: str
     parser: Union[GmshParser, CustomParser]
+    element_order: int
 
     def __init__(
         self,
         file_path: str,
+        element_order: int
     ):
         if not os.path.isfile(file_path):
             raise IOError(f"Mesh file {file_path} not found.")
 
         self.mesh_file_path = file_path
+        self.element_order = element_order
 
         # Check gmsh version
         # If version2  -> custom gmsh parser
@@ -46,9 +50,9 @@ def __init__(
             self.file_version = version
 
         if version.startswith("2"):
-            self.parser = CustomParser(file_path)
+            self.parser = CustomParser(file_path, element_order)
         else:
-            self.parser = GmshParser(file_path)
+            self.parser = GmshParser(file_path, element_order)
 
     def get_mesh_nodes_and_elements(self) -> Tuple[npt.NDArray, npt.NDArray]:
         return self.parser.get_mesh_nodes_and_elements()
@@ -102,7 +106,8 @@ def generate_rectangular_mesh(
         width: float = 0.005,
         height: float = 0.001,
         mesh_size: float = 0.00015,
-        x_offset: float = 0
+        x_offset: float = 0,
+        element_order: int = 1
     ):
         """Creates a gmsh rectangular mesh given the width, height, the mesh
         size and the x_offset.
@@ -121,6 +126,7 @@ def generate_rectangular_mesh(
             gmsh.initialize()
 
         gmsh.clear()
+        gmsh.option.setNumber("Mesh.ElementOrder", element_order)
 
         corner_points = [
             [x_offset, 0],