It is tempting to store submeshes as independent meshes.
However, this will cause quite a few issues for codim-1 problems down the line, as:
The way we should go is by store the submesh entities in a way that it relates to a sub entity of the parent mesh.
The idea would be to permute the function data so that it relates to a restriction of the dofs of a reference element of the parent.
I.e. if we have a second order Lagrange function on a triangular submesh, we create a second order Lagrange on the parent cell and use the subentity closure permutation to transfer the data to the reference of the parent mesh.
comm = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_unit_cube(
comm, 4,2,2, ghost_mode=dolfinx.mesh.GhostMode.shared_facet
)
def locate_cells(x):
return np.isclose(x[0], 0.0)#(x[0] - 0.2) ** 2 + (x[1] - 0.2) ** 2 < 0.2
dim = mesh.topology.dim-1
cells = dolfinx.mesh.locate_entities(mesh, dim, locate_cells)
submesh, cell_map, _, node_map = dolfinx.mesh.create_submesh(mesh, dim, cells)
el = ("Lagrange", 2)
V_sub = dolfinx.fem.functionspace(submesh, el)
import basix.ufl
sub_e = V_sub.element.basix_element
p_e = basix.ufl.element(sub_e.family, mesh.ufl_domain().ufl_coordinate_element().basix_element.cell_type, sub_e.degree, lagrange_variant=sub_e.lagrange_variant, discontinuous=False)
num_submesh_cells = submesh.topology.index_map(submesh.topology.dim).size_local
mesh.topology.create_entity_permutations()
cell_info = mesh.topology.get_cell_permutation_info()
parent_facets = cell_map.sub_topology_to_topology(np.arange(num_submesh_cells, dtype=np.int64), False)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)
itg_data = np.full((num_submesh_cells, 2),-1, dtype=np.int32)
for i, facet in enumerate(parent_facets):
cell = f_to_c.links(facet)[0]
facets = c_to_f.links(cell)
pos = np.where(facets == facet)[0][0]
itg_data[i, :] = (cell, pos)
for cell, pos in itg_data:
permuting = np.arange(V_sub.element.basix_element.dim, dtype=np.int32)
p_e.basix_element.permute_subentity_closure(permuting, cell_info[cell], V_sub.element.basix_element.cell_type, int(pos))
print(cell, pos, permuting)
It is tempting to store submeshes as independent meshes.
However, this will cause quite a few issues for codim-1 problems down the line, as:
The way we should go is by store the submesh entities in a way that it relates to a sub entity of the parent mesh.
The idea would be to permute the function data so that it relates to a restriction of the dofs of a reference element of the parent.
I.e. if we have a second order Lagrange function on a triangular submesh, we create a second order Lagrange on the parent cell and use the subentity closure permutation to transfer the data to the reference of the parent mesh.