Fixed the labelling of cells in the FaceCellularDomain class.

src/bvpy/templates/domains/cellular.py

@@ -81,21 +81,18 @@ class FaceCellularDomain(AbstractDomain, OccModel):
                                     for ind_p, p in enumerate(c)]), ind_c)
          for ind_c, c in enumerate(cells)]
 
-        # for ind_c, c in enumerate(cells):
-        #     size = len(c)
-        #     lines = [self.factory.addLine(self.points[c[ind_p]], self.points[c[(ind_p+1)%size]]) for ind_p, p in enumerate(c)]
-        #     s = self.factory.addCurveLoop(lines)
-        #     self.factory.addSurfaceFilling([s], ind_c)
-
         self.factory.removeAllDuplicates()
 
         self.synchronize()
 
         if markers is None:
-            [self.model.addPhysicalGroup(2, [ind], ind) for ind in range(len(cells))]
+            [self.model.addPhysicalGroup(2, [ind], ind)
+             for ind in range(len(cells))]
         else:
-            [self.model.addPhysicalGroup(2, [ind], marker) for ind, marker in enumerate(markers)]
-
+            [self.model.addPhysicalGroup(2, [ind for ind, marker
+                                             in enumerate(markers)
+                                             if marker == val], val)
+             for val in list(set(markers))]
 
 
 if __name__ == "__main__":

src/bvpy/templates/vforms/linear_elasticity.py

@@ -36,7 +36,6 @@ class ElasticForm(AbstractVform):
         self.traction_force = None  # useful if Neumann condition are inforced
 
         if isinstance(source, (int, float)):
-            # self._source_term = fe.as_vector([0, 0, source])
             self._source_term = fe.Constant([0, 0, source])
         else:
             self._source_term = source

src/bvpy/utils/visu.py

@@ -78,25 +78,24 @@ def _surface_plot_function(function, colorscale):
 
     return surface
 
+
 def _surface_plot_mesh(mesh, color, opacity=0.9):
     coord = mesh.coordinates()
     triangle = mesh.cells().T
 
-    if len(coord[0,:])==2:
-        coord = np.c_[coord, np.zeros(len(coord[:,0]))]
+    if len(coord[0, :]) == 2:
+        coord = np.c_[coord, np.zeros(len(coord[:, 0]))]
 
-    surface = go.Mesh3d(
-        x=coord[:,0],
-        y=coord[:,1],
-        z=coord[:,2],
-        i=triangle[0, :],
-        j=triangle[1, :],
-        k=triangle[2, :],
-        flatshading=True,
-        color=color,
-        opacity=opacity,
-        lighting=dict(ambient=1)
-        )
+    surface = go.Mesh3d(x=coord[:, 0],
+                        y=coord[:, 1],
+                        z=coord[:, 2],
+                        i=triangle[0, :],
+                        j=triangle[1, :],
+                        k=triangle[2, :],
+                        flatshading=True,
+                        color=color,
+                        opacity=opacity,
+                        lighting=dict(ambient=1))
 
     return surface
 
@@ -179,16 +178,18 @@ def plot(obj, colorscale="inferno", wireframe=True, size=10, norm=False,
         if obj.ufl_shape == ():
             surface = _surface_plot_function(obj, colorscale=colorscale)
             data.append(surface)
-        elif obj.ufl_shape == (3,):
+        elif obj.ufl_shape == (3, ):
             if norm:
                 V = obj.function_space().split()[0].collapse()
-                norm = fe.project(fe.sqrt(fe.inner(obj,obj)), V)
+                norm = fe.project(fe.sqrt(fe.inner(obj, obj)), V)
                 surface = _surface_plot_function(norm, colorscale=colorscale)
                 data.append(surface)
 
             func = FenicsFunctionConverter(obj)
-            point_val = np.array([[i['position'], i['value']] for i in func.to_vertex_values().values()])
-            cones = _cone_plot(point_val[:,0], point_val[:,1], size=size)
+            point_val = np.array([[i['position'],
+                                   i['value']]
+                                  for i in func.to_vertex_values().values()])
+            cones = _cone_plot(point_val[:, 0], point_val[:, 1], size=size)
             data.append(cones)
 
         if wireframe: