Vertex submesh support (#3455)

* Add tests for submesh of vertices

* Redirect basix branch

* Ruff formatting

* Add fixes and simplifications for submesh of vertices

* Tidy up test

* Ignore when edim becomes negative by setting to zero in parametrized test

* Revert branches

* Remove formatting error

cpp/dolfinx/mesh/utils.h

@@ -1250,7 +1250,12 @@ create_subgeometry(const Mesh<T>& mesh, int dim,
   // Create sub-geometry coordinate element
   CellType sub_coord_cell
       = cell_entity_type(geometry.cmap().cell_shape(), dim, 0);
-  fem::CoordinateElement<T> sub_cmap(sub_coord_cell, geometry.cmap().degree(),
+  // Special handling if point meshes, as they only support constant basis
+  // functions
+  int degree = geometry.cmap().degree();
+  if (sub_coord_cell == CellType::point)
+    degree = 0;
+  fem::CoordinateElement<T> sub_cmap(sub_coord_cell, degree,
                                      geometry.cmap().variant());
 
   // Sub-geometry input_global_indices

python/dolfinx/mesh.py

@@ -675,11 +675,10 @@ def create_submesh(
     submsh, entity_map, vertex_map, geom_map = _cpp.mesh.create_submesh(
         msh._cpp_object, dim, entities
     )
-    submsh_ufl_cell = ufl.Cell(to_string(submsh.topology.cell_type))
     submsh_domain = ufl.Mesh(
         basix.ufl.element(
             "Lagrange",
-            submsh_ufl_cell.cellname(),
+            to_string(submsh.topology.cell_type),
             submsh.geometry.cmap.degree,
             basix.LagrangeVariant(submsh.geometry.cmap.variant),
             shape=(submsh.geometry.dim,),

python/test/unit/mesh/test_mesh.py

@@ -503,35 +503,37 @@ def boundary_2(x):
 
 
 # TODO Test that submesh of full mesh is a copy of the mesh
-@pytest.mark.parametrize("d", [2, 3])
+@pytest.mark.parametrize("d", [1, 2, 3])
 @pytest.mark.parametrize("n", [3, 6])
 @pytest.mark.parametrize("codim", [0, 1, 2])
 @pytest.mark.parametrize("marker", [lambda x: x[0] >= 0.5, lambda x: x[0] >= -1])
 @pytest.mark.parametrize("ghost_mode", [GhostMode.none, GhostMode.shared_facet])
 @pytest.mark.parametrize("simplex", [True, False])
 def test_submesh_full(d, n, codim, marker, ghost_mode, simplex):
-    if d == codim:
-        pytest.xfail("Cannot create vertex submesh")
-    if d == 2:
+    if d == 1:
+        mesh = create_unit_interval(MPI.COMM_WORLD, n, ghost_mode=ghost_mode)
+    elif d == 2:
         ct = CellType.triangle if simplex else CellType.quadrilateral
         mesh = create_unit_square(MPI.COMM_WORLD, n, n, ghost_mode=ghost_mode, cell_type=ct)
     else:
         ct = CellType.tetrahedron if simplex else CellType.hexahedron
         mesh = create_unit_cube(MPI.COMM_WORLD, n, n, n, ghost_mode=ghost_mode, cell_type=ct)
 
-    edim = mesh.topology.dim - codim
+    edim = max(mesh.topology.dim - codim, 0)
     entities = locate_entities(mesh, edim, marker)
     submesh, entity_map, vertex_map, geom_map = create_submesh(mesh, edim, entities)
     submesh_topology_test(mesh, submesh, entity_map, vertex_map, edim)
     submesh_geometry_test(mesh, submesh, entity_map, geom_map, edim)
 
 
-@pytest.mark.parametrize("d", [2, 3])
+@pytest.mark.parametrize("d", [1, 2, 3])
 @pytest.mark.parametrize("n", [3, 6])
 @pytest.mark.parametrize("boundary", [boundary_0, boundary_1, boundary_2])
 @pytest.mark.parametrize("ghost_mode", [GhostMode.none, GhostMode.shared_facet])
 def test_submesh_boundary(d, n, boundary, ghost_mode):
-    if d == 2:
+    if d == 1:
+        mesh = create_unit_interval(MPI.COMM_WORLD, n, ghost_mode=ghost_mode)
+    elif d == 2:
         mesh = create_unit_square(MPI.COMM_WORLD, n, n, ghost_mode=ghost_mode)
     else:
         mesh = create_unit_cube(MPI.COMM_WORLD, n, n, n, ghost_mode=ghost_mode)