Strange behavior after using create_mesh()

dolfinx
,
1

Hello,

I am using the Docker container of DOLFINx-v0.8.0.
I have a global mesh that is divided into two regions by cell markers. For visualization, I have to store the two sub-meshes. I wrote a function that loops over the connectivity and creates two new meshes, which works fine.
The problem: When creating the new meshes with create_mesh(), I get strange errors in the code later on. I created a minimal example from the Poisson demo. This is just copy&paste, and I added the one line random_mesh = create_mesh(...) after the generation of the computational mesh.

from mpi4py import MPI
from petsc4py.PETSc import ScalarType  # type: ignore
import numpy as np
import ufl
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import ds, dx, grad, inner

msh = mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=((0.0, 0.0), (2.0, 1.0)),
    n=(32, 16),
    cell_type=mesh.CellType.triangle,
)

# Comment this in or out to see the error
random_mesh = mesh.create_mesh(msh.comm, [[0, 1, 2]], msh.geometry.x[:, 0:2], msh.ufl_domain())

V = fem.functionspace(msh, ("Lagrange", 1))
facets = mesh.locate_entities_boundary(
    msh,
    dim=(msh.topology.dim - 1),
    marker=lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], 2.0),
)
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx + inner(g, v) * ds
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
with io.XDMFFile(msh.comm, "out_poisson/poisson.xdmf", "w") as file:
    file.write_mesh(msh)
    file.write_function(uh)

This new mesh is just a single triangle using the first 3 points of the computational mesh. I am not doing anything with this newly created mesh, but the LinearProblem class throws RuntimeError: Incompatible mesh. entity_maps must be provided. Comment it out and errors are gone. Is there something I am overlooking?

Thank you!

2

I woudl guess an issue here is that the msh.geometry.x has way more nodes than [0,1,2], making it inconsistent.

Note that you can create submesh with dolfinx.mesh.create_submesh dolfinx/python/dolfinx/mesh.py at bd5ba2ce6d3ac9ecfe3324d83dacd2c7be4de5ed · FEniCS/dolfinx · GitHub that does not have additional nodes in them.

1 Like
3

Thank you for the hint with create_submesh(), I will have a look.

Regarding the create_mesh(): Having a larger points array does not influence the mesh (it only picks according to the connectivity), the problem is the usage of msh.ufl_domain(). If I change the code to

random_mesh = mesh.create_mesh(msh.comm, [[0, 1, 2]], msh.geometry.x[:, 0:2], 
   ufl.Mesh(basix.ufl.element("Lagrange", 'triangle', 1, shape=(2,)))
)

it works. It seems as the create_mesh() routine changes something to the underlying ufl.Mesh element from the parent mesh. I do not know if this is intended behavior though.

4

The problem is that ufl_domain() has ufl_cargo() attached (circular reference):

print(id(msh.ufl_domain().ufl_cargo()))
random_mesh = mesh.create_mesh(msh.comm, [[0, 1, 2]], msh.geometry.x[:, 0:2], msh.ufl_domain())
print(id(random_mesh.ufl_domain().ufl_cargo()))
print(id(msh.ufl_domain().ufl_cargo()))

which is set at