Reading and writing meshtags on 2rd order mesh with adios4dolfinx

dolfinx
,
1

I changed the mesh used here https://github.com/jorgensd/adios4dolfinx/blob/main/docs/meshtags.py
from dolfinx buidin mesh to Gmsh mesh.

Here is the code:

import logging
from pathlib import Path
from mpi4py import MPI
import dolfinx
import numpy as np
import adios4dolfinx

assert MPI.COMM_WORLD.size == 1, "This example should only be run with 1 MPI process"

# mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, nx=3, ny=4, nz=5)

import gmsh
gdim = 3
marker = 3
mesh_size = 0.25
mesh_order = 2
gmsh.initialize()
factory = gmsh.model.occ
# add box
tag1 = factory.add_box(0, 0, 0, 1, 1, 1)
factory.synchronize()
ov = gmsh.model.get_entities(gdim)
gmsh.model.addPhysicalGroup(gdim, [ov[0][1]], marker, 'fluid')
gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size)
gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size)

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(mesh_order)

mesh, cell_tags, facet_tags = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim)


entity_midpoints = {}
meshtags = {}
for i in range(mesh.topology.dim + 1):
    mesh.topology.create_entities(i)
    e_map = mesh.topology.index_map(i)

    # Compute midpoints of entities
    entities = np.arange(e_map.size_local, dtype=np.int32)
    mesh.topology.create_connectivity(i, mesh.topology.dim)
    entity_midpoints[i] = dolfinx.mesh.compute_midpoints(mesh, i, entities)
    # Associate each local index with its global index
    values = np.arange(e_map.size_local, dtype=np.int32) + e_map.local_range[0]
    meshtags[i] = dolfinx.mesh.meshtags(mesh, i, entities, values)

filename = Path("mesh_with_meshtags.bp")
adios4dolfinx.write_mesh(filename, mesh)
for i, tag in meshtags.items():
    adios4dolfinx.write_meshtags(filename, mesh, tag, meshtag_name=f"meshtags_{i}")

def verify_meshtags(filename: Path):
    # We assume that entity_midpoints have been sent to the engine
    from mpi4py import MPI

    import dolfinx
    import numpy as np

    import adios4dolfinx

    read_mesh = adios4dolfinx.read_mesh(filename, MPI.COMM_WORLD)
    prefix = f"{read_mesh.comm.rank + 1}/{read_mesh.comm.size}: "
    for i in range(read_mesh.topology.dim + 1):
        # Read mesh from file
        meshtags = adios4dolfinx.read_meshtags(filename, read_mesh, meshtag_name=f"meshtags_{i}")

        # Compute midpoints for all local entities on process
        read_mesh.topology.create_connectivity(i, read_mesh.topology.dim)
        midpoints = dolfinx.mesh.compute_midpoints(read_mesh, i, meshtags.indices)
        # Compare locally computed midpoint with reference data
        for global_pos, midpoint in zip(meshtags.values, midpoints):
            np.testing.assert_allclose(
                entity_midpoints[i][global_pos],
                midpoint,
                err_msg=f"{prefix}: Midpoint ({i , global_pos}) do not match",
            )
        print(f"{prefix} Matching of all entities of dimension {i} successful")

verify_meshtags(filename)

If mesh_order=1 (line 16), the code runs fine. However, if mesh_order=2, it shows some error

image
image1102×201 23.5 KB

Is this an issue with dolfinx v0.8.0 and adios4dolifnx v0.8.1?

2

Please upgrade to v0.9.0, as the issue seems to be resolved there.

3

Thanks you. I will try it in 0.9.0 later.