Simplicial complex

src/dctkit/mesh/simplex.py

@@ -6,7 +6,7 @@
 import numpy.typing as npt
 from jax import Array
 import jax.numpy as jnp
-from typing import Tuple, Any
+from typing import Tuple, Any, List
 
 
 class SimplicialComplex:
@@ -39,12 +39,12 @@ class SimplicialComplex:
     """
 
     def __init__(self, tet_node_tags: npt.NDArray, node_coords: npt.NDArray,
-                 is_well_centered: bool = False):
+                 space_dim: int = 3, is_well_centered: bool = False):
 
-        self.node_coords = node_coords.astype(dctkit.float_dtype)
+        self.node_coords = node_coords.astype(dctkit.float_dtype)[:, :space_dim]
         tet_node_tags = tet_node_tags.astype(dctkit.int_dtype)
         self.num_nodes = node_coords.shape[0]
-        self.space_dim = node_coords.shape[1]
+        self.space_dim = space_dim
         self.float_dtype = dctkit.float_dtype
         self.int_dtype = dctkit.int_dtype
         self.is_well_centered = is_well_centered
@@ -92,7 +92,8 @@ def get_tets_containing_a_boundary_face(self):
         if not hasattr(self, "bnd_faces_indices"):
             self.get_complex_boundary_faces_indices()
         dim = self.dim - 1
-        self.tets_cont_bnd_face = get_cofaces(self.bnd_faces_indices, dim, self)
+        self.tets_cont_bnd_face = get_cofaces(
+            self.bnd_faces_indices, dim, self)
 
     def get_circumcenters(self):
         """Compute all the circumcenters."""
@@ -214,7 +215,11 @@ def get_dual_edge_vectors(self):
         # boundary faces arranged by rows, padded with zeros for the non-boundary edges
         if not hasattr(self, "bnd_faces_indices"):
             self.get_complex_boundary_faces_indices()
-        circ_faces = self.circ[dim-1]
+        if dim == 1:
+            # in this case faces = nodes
+            circ_faces = self.node_coords
+        else:
+            circ_faces = self.circ[dim-1]
         circ_bnd_faces = np.zeros(circ_faces.shape, dtype=dctkit.float_dtype)
         circ_bnd_faces[self.bnd_faces_indices] = circ_faces[self.bnd_faces_indices]
 
@@ -424,7 +429,7 @@ def get_deformation_gradient(self, node_coords: npt.NDArray) -> Array:
 
 
 def get_cofaces(faces_ids: list[int] | npt.NDArray, faces_dim: int,
-                S: SimplicialComplex) -> list[npt.NDArray]:
+                S: SimplicialComplex) -> List[npt.NDArray]:
     """Get the IDs of the cofaces of a simplex, i.e. the neighour simplices of one
     higher dimension.
 

src/dctkit/mesh/util.py

@@ -5,7 +5,8 @@
 from meshio import Mesh
 
 
-def build_complex_from_mesh(mesh: Mesh, is_well_centered=True) -> SimplicialComplex:
+def build_complex_from_mesh(mesh: Mesh, space_dim=3,
+                            is_well_centered=True) -> SimplicialComplex:
     """Build a SimplicialComplex object from a meshio.Mesh object.
 
     Args:
@@ -26,7 +27,8 @@ def build_complex_from_mesh(mesh: Mesh, is_well_centered=True) -> SimplicialComp
     elif "line" in cell_types:
         tet_node_tags = mesh.cells_dict["line"]
 
-    S = SimplicialComplex(tet_node_tags, node_coords, is_well_centered=is_well_centered)
+    S = SimplicialComplex(tet_node_tags, node_coords, space_dim,
+                          is_well_centered=is_well_centered)
     return S