Mesh and MeshTag connection

mesh
1

Hello Fenics community,

I am trying to do some changes in my mesh. Therefore i have to change some connections of the cells. First i wanted to create a mesh and meshtag out of the existing mesh without changes.

The mesh was initially created with gmsh and is loaded with following command:

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "r") as xdmf:
    self.mesh = xdmf.read_mesh(name="Grid")
    self.subdomains = xdmf.read_meshtags(msh, name='Grid')

From the connections and meshtag information i wanted to create the same mesh again:

cells = self.subdomains.topology.connectivity(2, 0)

        con_vec = []
        entity_vec = []
        value_vec = []
        # Get dof layout for each cell
        num_entities = self.mesh.topology.index_map(2).size_local + self.mesh.topology.index_map(2).num_ghosts
        for entity in range(num_entities):
            entities = cells.links(entity)
            con_vec.append(entities)
            ind = self.subdomains.indices[entity]
            entity_vec.append(ind)
            value_vec.append(self.subdomains.values[ind])

        entity_vec = np.array(entity_vec, dtype=np.int32)
        con_vec = np.array(con_vec, dtype=np.int64)

        xpos = []
        for i, xi in enumerate(self.mesh.geometry.x[:]):
            xpos.append([xi[0], xi[1]])

        self.mesh = create_mesh(MPI.COMM_WORLD, con_vec, xpos, self.mesh._ufl_domain)
        self.subdomains = meshtags(self.mesh, 2, entity_vec, value_vec)

The mesh itself looks good and identical to before. But the meshtag are not connected to the right cell positions.

I am sure i am doing something wrong, but can’t understand what and why it is not working. Sry in the case i am doing something complete wrong, Still a beginner in this topics :wink:
I am using DOLFINx version: 0.7.2

Thanks for your help.

2

As your code is not a reproducible example, it is hard to give guidance. However, I think you should consider Ordering of cells via mesh.topology.original_cell_index - #2 by dokken using the original cell index map to transfer the entity tags from its input index.

3

Thanks for the quick respond. Unfortunately that wasn’t effective.

I have uploaded an example of my code here:

Code_Example

Normally I would expect to get the same result by activating Line 58.

4

Consider the following code:

import dolfinx
import ufl
from mpi4py import MPI
from dolfinx.fem import Function, VectorFunctionSpace
from petsc4py import PETSc
from dolfinx.mesh import create_mesh, meshtags
import numpy as np
import dolfinx.io

print(f"DOLFINx version: {dolfinx.__version__} is installed")
print(f"ufl version: {ufl.__version__} is installed")

filename = '00_ModelMesh_MAG.xdmf'

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    subdomains = xdmf.read_meshtags(mesh, name='Grid')

CG_B = VectorFunctionSpace(mesh, ("DG", 0))
M = Function(CG_B)

cells_magnet = subdomains.indices[subdomains.values == 3]
M.x.array[cells_magnet *
          2] = np.full_like(cells_magnet, 1000000, dtype=PETSc.ScalarType)
M.x.array[cells_magnet * 2 +
          1] = np.full_like(cells_magnet, 0, dtype=PETSc.ScalarType)

with dolfinx.io.XDMFFile(mesh.comm, "m.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(M)

con_vec = []
value_vec = []
# Get dof layout for each cell
num_entities = mesh.topology.index_map(
    2).size_local + mesh.topology.index_map(2).num_ghosts

xpos = mesh.geometry.x
geometry_dofs = dolfinx.cpp.mesh.entities_to_geometry(
    mesh._cpp_object, 2, np.arange(num_entities, dtype=np.int32), False)
for cell in range(num_entities):
    # Map vertices to geometry nodes
    con_vec.append(geometry_dofs[cell])
    value_vec.append(subdomains.values[cell])

con_vec = np.array(con_vec, dtype=np.int64)


mesh2 = create_mesh(MPI.COMM_WORLD, con_vec,
                    xpos[:, :mesh.geometry.dim], mesh._ufl_domain)
o_index = mesh2.topology.original_cell_index

cells = np.arange(mesh2.topology.index_map(2).size_local, dtype=np.int32)
subdomains = meshtags(mesh2, 2, cells, np.asarray(value_vec)[o_index])

CG_B = VectorFunctionSpace(mesh2, ("DG", 0))
M2 = Function(CG_B)

cells_magnet = subdomains.indices[subdomains.values == 3]
M2.x.array[cells_magnet *
           2] = np.full_like(cells_magnet, 1000000, dtype=PETSc.ScalarType)
M2.x.array[cells_magnet * 2 +
           1] = np.full_like(cells_magnet, 0, dtype=PETSc.ScalarType)

with dolfinx.io.XDMFFile(mesh2.comm, "M2.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh2)
    xdmf.write_function(M2)

You were adding underlying assumptions between the mesh geometry and mesh topology.

5

Hi Dokken,

thank you very much, with this my problem is solved.