Debug most of adding input face

src/zmsh/delaunay.py

@@ -121,7 +121,7 @@ def __init__(self, geometry: Geometry):
 
         # Store which numeric IDs can still be assigned to cells of each
         # dimension
-        self._free_cell_ids = [set() for k in range(dimension + 1)]
+        self._free_cell_ids = [set() for k in range(geometry.dimension + 1)]
 
     @property
     def face_queue(self):
@@ -152,6 +152,9 @@ def _get_new_cell_ids(self, dimension: int, num_new_cells: int):
             self.geometry.topology.cells(dimension).resize(2**exp * num_cells)
             free_cell_ids.update(set(range(num_cells, 2**exp * num_cells)))
 
+        if num_new_cells == 1:
+            return free_cell_ids.pop()
+
         return [free_cell_ids.pop() for i in range(num_new_cells)]
 
     def is_done(self):
@@ -198,29 +201,51 @@ def find_crossing_faces(self, input_face_id):
         return result
 
     def add_input_face(self, input_face_id):
-        # Find all the faces that cross the input face
+        # Get the incidence matrix for the input face
+        # NOTE: This only works in 2D, we need to be more careful about the
+        # correspondence between cell IDs in the input and result geometries
+        # in higher dimensions
         dimension = self.geometry.dimension
+        faces = self._input_geometry.topology.cells(dimension - 1)
+        vertex_ids, signs = faces[input_face_id]
+
+        # Find all the faces that cross the input face
         cofaces = self.geometry.topology.cocells(dimension - 1)
         crossing_face_ids = self.find_crossing_faces(input_face_id)
+
+        # Find all the cells containing the crossing faces
         cell_ids = cofaces[crossing_face_ids][0]
         cells = self.geometry.topology.cells(dimension)
-        face_ids = cells[cell_id][0]
+        face_ids = cells[cell_ids][0]
 
         # Merge these cells
-        raise NotImplementedError("Woops haven't done this yet")
-        self.free_cell_ids(k).update(crossing_face_ids)
+        face_ids, D = cells[cell_ids]
+        cells[face_ids, cell_ids[0]] = np.sum(D, axis=1)
+        cells[face_ids, cell_ids[1:]] = np.zeros_like(D[:, 1:], dtype=np.int8)
 
-        # Add the new face
+        # Removing the crossing faces
         faces = self.geometry.topology.cells(dimension - 1)
-        face_id = self._get_new_cell_ids(dimension - 1, 1)
-        faces[subface_ids, face_id] = signs
-        return face_id
+        subface_ids, E = faces[crossing_face_ids]
+        faces[subface_ids, crossing_face_ids] = np.zeros_like(E, dtype=np.int8)
+        self.free_cell_ids(dimension - 1).update(crossing_face_ids)
 
-    def retriangulate_polygon(self, cell_id, face_ids):
+        # Add the new face
+        new_face_id = self._get_new_cell_ids(dimension - 1, 1)
+        faces[vertex_ids, new_face_id] = signs
+        boundary_face_ids = np.setdiff1d(face_ids, crossing_face_ids)
+        return new_face_id, boundary_face_ids
+
+    def retriangulate_polygon(self, base_face_id, face_ids):
         pass
 
     def step(self):
-        raise NotImplementedError("Haven't got here yet!")
+        if self.retriangulation_queue:
+            face_id, boundary_face_ids = self.retriangulation_queue.pop()
+            return self.retriangulate_polygon(self, face_id, boundary_face_ids)
+        elif self.face_queue:
+            input_face_id = self.face_queue.pop()
+            new_face_id, boundary_face_ids = self.add_input_face(input_face_id)
+            self.retriangulation_queue.push((new_face_id, boundary_face_ids))
 
     def finalize(self):
         raise NotImplementedError("Haven't got here yet!")