From 99e29fbb1b97481bd23bb0e647883327adeb29ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Mon, 7 Oct 2024 13:23:00 +0200 Subject: [PATCH] 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 | 7 ++++++- python/dolfinx/mesh.py | 3 +-- python/test/unit/mesh/test_mesh.py | 16 +++++++++------- 3 files changed, 16 insertions(+), 10 deletions(-) diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 7479ef67bcd..1fa5724e39d 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1250,7 +1250,12 @@ create_subgeometry(const Mesh& mesh, int dim, // Create sub-geometry coordinate element CellType sub_coord_cell = cell_entity_type(geometry.cmap().cell_shape(), dim, 0); - fem::CoordinateElement 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 sub_cmap(sub_coord_cell, degree, geometry.cmap().variant()); // Sub-geometry input_global_indices diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 2203ca8d961..e58563a2bf1 100644 --- a/python/dolfinx/mesh.py +++ b/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,), diff --git a/python/test/unit/mesh/test_mesh.py b/python/test/unit/mesh/test_mesh.py index bc9bfeb68bb..e620e1c0569 100644 --- a/python/test/unit/mesh/test_mesh.py +++ b/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)