[WIP] Implemented get_tet_mesh and fixed dual volumes bug. Fixing tes…

…t_simplex.py

src/dctkit/mesh/simplex.py

@@ -107,7 +107,9 @@ def get_primal_volumes(self):
         """Compute all the primal volumes."""
         # loop over all p-simplices (1..dim + 1)
         # (volume of 0-simplices is 1, we do not store it)
-        self.primal_volumes = sl.ShiftedList([None] * (self.dim), -1)
+        # self.primal_volumes = sl.ShiftedList([None] * (self.dim), -1)
+        self.primal_volumes = [None]*(self.dim + 1)
+        self.primal_volumes[0] = np.ones(self.num_nodes, dtype=self.float_dtype)
         for p in range(1, self.dim + 1):
             S = self.S[p]
             num_p_simplices, _ = S.shape
@@ -120,7 +122,8 @@ def get_primal_volumes(self):
 
     def get_dual_volumes(self):
         """Compute all the dual volumes."""
-        self.dual_volumes = sl.ShiftedList([None] * (self.dim), -1)
+        # self.dual_volumes = sl.ShiftedList([None] * (self.dim), -1)
+        self.dual_volumes = [None] * (self.dim+1)
         self.dual_volumes[self.dim] = np.ones(self.S[self.dim].
                                               shape[0], dtype=self.float_dtype)
         # loop over simplices at all dimensions
@@ -169,18 +172,18 @@ def get_hodge_star(self):
         if self.is_well_centered:
             self.hodge_star_inverse = [None]*(n + 1)
         for p in range(n + 1):
-            if p == 0:
-                # volumes of vertices are 1 by definition
-                pv = 1
-                dv = self.dual_volumes[n - p]
-            elif p == n:
-                pv = self.primal_volumes[p]
-                # volumes of vertices are 1 by definition
-                dv = 1
-            else:
-                pv = self.primal_volumes[p]
-                dv = self.dual_volumes[p]
-            self.hodge_star[p] = dv/pv
+            # if p == 0:
+            #    # volumes of vertices are 1 by definition
+            #    pv = 1
+            #    dv = self.dual_volumes[p]
+            # elif p == n:
+            #    pv = self.primal_volumes[p]
+            #    # volumes of vertices are 1 by definition
+            #    dv = 1
+            # else:
+            #    pv = self.primal_volumes[p]
+            #    dv = self.dual_volumes[p]
+            self.hodge_star[p] = self.dual_volumes[p]/self.primal_volumes[p]
             if self.is_well_centered:
                 self.hodge_star_inverse[p] = 1.0/self.hodge_star[p]
                 # adjust the sign in order to have star_inv*star = (-1)^(p*(n-p))

src/dctkit/mesh/util.py

@@ -2,11 +2,12 @@
 import gmsh  # type: ignore
 import numpy as np
 import numpy.typing as npt
-from typing import Tuple
+from typing import Tuple, Dict
 
 
-def read_mesh(filename: str = None, format: str = "gmsh") -> Tuple[
-        int, int, npt.NDArray, npt.NDArray, npt.NDArray]:
+def read_mesh(tet_type: Dict, filename: str = None,
+              format: str = "gmsh") -> Tuple[int, int, npt.NDArray,
+                                             npt.NDArray, npt.NDArray]:
     """Reads a mesh from file.
 
     Args:
@@ -27,46 +28,53 @@ def read_mesh(filename: str = None, format: str = "gmsh") -> Tuple[
     if filename is not None:
         gmsh.open(filename)
 
+    # FIXME: fix docs.
+    familyName = tet_type["familyName"]
+    numNodesPerTet = tet_type["numNodesPerTet"]
+    elementType = gmsh.model.mesh.getElementType(familyName, 1)
+
     # Get nodes and corresponding coordinates
     nodeTags, coords, _ = gmsh.model.mesh.getNodes()
     numNodes = len(nodeTags)
 
     coords = np.array(coords, dtype=float_dtype)
     # Get 2D elements and associated node tags
-    elemTags, nodeTagsPerElem = gmsh.model.mesh.getElementsByType(2)
+    # print(gmsh.model.mesh.getElementsByType(4))
+    elemTags, nodeTagsPerElem = gmsh.model.mesh.getElementsByType(elementType)
 
     # Decrease element IDs by 1 to have node indices starting from 0
     nodeTagsPerElem = np.array(nodeTagsPerElem, dtype=int_dtype) - 1
-    nodeTagsPerElem = nodeTagsPerElem.reshape(len(nodeTagsPerElem) // 3, 3)
+    nodeTagsPerElem = nodeTagsPerElem.reshape(
+        len(nodeTagsPerElem) // numNodesPerTet, numNodesPerTet)
     # Get number of TRIANGLES
     numElements = len(elemTags)
 
     # Position vectors of mesh points
     node_coords = coords.reshape(len(coords)//3, 3)
 
-    # get node tags per boundary elements
-    _, nodeTagsPerBElem = gmsh.model.mesh.getElementsByType(1)
+    # get node tags per boundary element
+    if elementType == 4:
+        bndElementType = 2
+        numNodesPerBndElem = 3
+    else:
+        bndElementType = elementType - 1
+        numNodesPerBndElem = elementType
+    _, nodeTagsPerBElem = gmsh.model.mesh.getElementsByType(bndElementType)
     nodeTagsPerBElem = np.array(nodeTagsPerBElem, dtype=int_dtype) - 1
-    nodeTagsPerBElem = nodeTagsPerBElem.reshape(len(nodeTagsPerBElem) // 2, 2)
+    nodeTagsPerBElem = nodeTagsPerBElem.reshape(
+        len(nodeTagsPerBElem) // numNodesPerBndElem, numNodesPerBndElem)
     # we sort every row to get the orientation used for our simulations
     nodeTagsPerBElem = np.sort(nodeTagsPerBElem)
 
     return numNodes, numElements, nodeTagsPerElem, node_coords, nodeTagsPerBElem
 
 
-def generate_square_mesh(lc: float) -> Tuple[int, int, npt.NDArray,
-                                             npt.NDArray, npt.NDArray]:
+def generate_square_mesh(tet_type: Dict, lc: float) -> None:
     """ Generate a simple square mesh.
 
     Args:
         lc: target mesh size (lc) close to a given point.
 
-    Returns:
-        a tuple containing the number of mesh nodes; the number of faces; the matrix
-        containing the IDs of the nodes (cols) belonging to each face (rows); the node
-        coordinates; the matrix containing the IDs of the nodes (cols) belonging to
-        each boundary element (rows).
-
     """
     if not gmsh.is_initialized():
         gmsh.initialize()
@@ -89,13 +97,56 @@ def generate_square_mesh(lc: float) -> Tuple[int, int, npt.NDArray,
     gmsh.model.addPhysicalGroup(1, [4], 3, name="right")
     gmsh.model.mesh.generate(2)
 
-    numNodes, numElements, nodeTagsPerElem, node_coords, nodeTagsPerBElem = read_mesh()
+    numNodes, numElements, nodeTagsPerElem, node_coords, nodeTagsPerBElem = read_mesh(tet_type)
 
     return numNodes, numElements, nodeTagsPerElem, node_coords, nodeTagsPerBElem
 
 
-def generate_hexagon_mesh(a: float, lc: float) -> Tuple[int, int, npt.NDArray,
-                                                        npt.NDArray, npt.NDArray]:
+def generate_tet_mesh(lc: float) -> None:
+    """ Generate a simple square mesh.
+
+    Args:
+        lc: target mesh size (lc) close to a given point.
+
+    Returns:
+        a tuple containing the number of mesh nodes; the number of faces; the matrix
+        containing the IDs of the nodes (cols) belonging to each face (rows); the node
+        coordinates; the matrix containing the IDs of the nodes (cols) belonging to
+        each boundary element (rows).
+
+    """
+    # FIXME: FIX DOCS.
+    if not gmsh.is_initialized():
+        gmsh.initialize()
+
+    gmsh.model.add("tet")
+    gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
+    gmsh.model.geo.addPoint(1, 0, 0, lc, 2)
+    gmsh.model.geo.addPoint(1/2, 1, 0, lc, 3)
+    gmsh.model.geo.addPoint(0, 0, 1, lc, 4)
+
+    gmsh.model.geo.addLine(1, 2, 1)
+    gmsh.model.geo.addLine(2, 3, 2)
+    gmsh.model.geo.addLine(3, 1, 3)
+    gmsh.model.geo.addLine(1, 4, 4)
+    gmsh.model.geo.addLine(2, 4, 5)
+    gmsh.model.geo.addLine(3, 4, 6)
+    gmsh.model.geo.addCurveLoop([1, 2, 3], 1)
+    gmsh.model.geo.addCurveLoop([1, 5, -4], 2)
+    gmsh.model.geo.addCurveLoop([-3, 6, -4], 3)
+    gmsh.model.geo.addCurveLoop([2, 6, -5], 4)
+    gmsh.model.geo.addPlaneSurface([1], 1)
+    gmsh.model.geo.addPlaneSurface([2], 2)
+    gmsh.model.geo.addPlaneSurface([3], 3)
+    gmsh.model.geo.addPlaneSurface([4], 4)
+    gmsh.model.geo.addSurfaceLoop([1, 2, 3, 4], 1)
+    gmsh.model.geo.addVolume([1], 1)
+    gmsh.model.geo.synchronize()
+    gmsh.model.mesh.generate(3)
+
+
+def generate_hexagon_mesh(tet_type: Dict, a: float, lc: float) -> Tuple[int, int, npt.NDArray,
+                                                                        npt.NDArray, npt.NDArray]:
     """Generate a regular hexagonal mesh
 
     Args:
@@ -132,7 +183,7 @@ def generate_hexagon_mesh(a: float, lc: float) -> Tuple[int, int, npt.NDArray,
     gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4, 5, 6], 1)
     gmsh.model.mesh.generate(2)
 
-    numNodes, numElements, nodeTagsPerElem, node_coords, nodeTagsPerBElem = read_mesh()
+    numNodes, numElements, nodeTagsPerElem, node_coords, nodeTagsPerBElem = read_mesh(tet_type)
 
     return numNodes, numElements, nodeTagsPerElem, node_coords, nodeTagsPerBElem
 
@@ -148,6 +199,7 @@ def generate_1_D_mesh(num_nodes: int, L: float) -> Tuple[npt.NDArray, npt.NDArra
             x,y coords) and a matrix containing the IDs of the nodes belonging to
             each 1-simplex.
     """
+    # FIXME: FIX THIS.
     node_coords = np.linspace(0, L, num=num_nodes)
     x = np.zeros((num_nodes, 2))
     x[:, 0] = node_coords

tests/test_circumcenter.py

@@ -11,7 +11,9 @@ def test_circumcenter(setup_test):
 
     filename = "data/test1.msh"
     full_path = os.path.join(cwd, filename)
-    numNodes, numElements, S_2, node_coord, _ = util.read_mesh(full_path)
+    tet_type = {"familyName": "Triangle",
+                "numNodesPerTet": 3}
+    numNodes, numElements, S_2, node_coord, _ = util.read_mesh(tet_type, full_path)
 
     print(f"The number of nodes in the mesh is {numNodes}")
     print(f"The number of faces in the mesh is {numElements}")

tests/test_cochain.py

@@ -14,7 +14,9 @@
 def test_cochain(setup_test):
     filename = "data/test1.msh"
     full_path = os.path.join(cwd, filename)
-    numNodes, numElements, S_2, node_coord, _ = util.read_mesh(full_path)
+    tet_type = {"familyName": "Triangle",
+                "numNodesPerTet": 3}
+    numNodes, numElements, S_2, node_coord, _ = util.read_mesh(tet_type, full_path)
 
     print(f"The number of nodes in the mesh is {numNodes}")
     print(f"The number of faces in the mesh is {numElements}")
@@ -59,7 +61,7 @@ def test_cochain(setup_test):
 
     # primal codifferential test (we need a well-centered mesh)
 
-    _, _, S_2_new, node_coords_new, _ = util.generate_square_mesh(0.4)
+    _, _, S_2_new, node_coords_new, _ = util.generate_square_mesh(tet_type, 0.4)
     triang = tri.Triangulation(node_coords_new[:, 0], node_coords_new[:, 1])
 
     plt.triplot(triang, 'ko-')

tests/test_linear_elasticity.py

@@ -7,8 +7,10 @@
 
 
 def test_linear_elasticity(setup_test):
+    tet_type = {"familyName": "Triangle",
+                "numNodesPerTet": 3}
     lc = 0.5
-    _, _, S_2, node_coords, bnd_faces_tags = util.generate_square_mesh(lc)
+    _, _, S_2, node_coords, bnd_faces_tags = util.generate_square_mesh(tet_type, lc)
     bnodes, _ = gmsh.model.mesh.getNodesForPhysicalGroup(1, 1)
     bnodes -= 1
     bnodes = bnodes.astype(dt.int_dtype)

tests/test_optctrl.py

@@ -87,8 +87,9 @@ def test_optimal_control_poisson(setup_test):
     np.random.seed(42)
 
     lc = 0.5
-
-    _, _, S_2, node_coord, _ = util.generate_square_mesh(lc)
+    tet_type = {"familyName": "Triangle",
+                "numNodesPerTet": 3}
+    _, _, S_2, node_coord, _ = util.generate_square_mesh(tet_type, lc)
 
     S, bnodes, _ = get_complex(S_2, node_coord)
 

tests/test_poisson.py

@@ -40,10 +40,11 @@ def test_poisson(setup_test, optimizer, energy_formulation):
         assert dt.float_dtype == "float64"
 
     np.random.seed(42)
-
+    tet_type = {"familyName": "Triangle",
+                "numNodesPerTet": 3}
     lc = 0.5
 
-    _, _, S_2, node_coord, _ = util.generate_square_mesh(lc)
+    _, _, S_2, node_coord, _ = util.generate_square_mesh(tet_type, lc)
 
     S, bnodes, _ = get_complex(S_2, node_coord)