Mesh refinement using dolfinx.mesh.refine_plaza()

mesh
1

Dear all, while using the above mentioned function in python interface of latest dolfinx v0.8.0 (conda), it seems to me that the returned parent_cells aren’t what they should be. Here’s a MWE to reproduce the results:

import dolfinx.mesh
import pyvista as pv
from dolfinx.cpp.refinement import RefinementOption
from mpi4py import MPI


def main():
    msh_c = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
    msh_c.topology.create_entities(1)
    msh_f, parent_cells, _ = dolfinx.mesh.refine_plaza(msh_c, option=RefinementOption.parent_cell)

    for child, parent in enumerate(parent_cells):
        print(f"Child {child} has parent {parent}")
        child_vertices = dolfinx.mesh.compute_incident_entities(msh_f.topology, child, 2, 0)
        child_vertex_coords = msh_f.geometry.x[child_vertices]
        print(f"Child {child} has vertices {child_vertices} at {child_vertex_coords}")
        parent_vertices = dolfinx.mesh.compute_incident_entities(msh_c.topology, parent, 2, 0)
        parent_vertex_coords = msh_c.geometry.x[parent_vertices]
        print(f"Parent {parent} has vertices {parent_vertices} at {parent_vertex_coords}")

    # Plot fine mesh
    grid = pv.UnstructuredGrid(*dolfinx.plot.vtk_mesh(msh_f, msh_f.topology.dim))
    plotter = pv.Plotter(window_size=[1000, 1000])
    plotter.show_axes()
    plotter.show_grid()
    plotter.add_mesh(grid, show_edges=True)
    plotter.show()

if __name__ == "__main__":
    main()

If reproducing the result, be sure to first manually locally apply the fix to refine_plaza(): Input for mesh refinement with refine_plaza - #4 by francesco-ballarin.
Now, as you can see after running the example, child 6 for example does not have the right parent, as the output contains:

…
Child 6 has parent 1
Child 6 has vertices [4 6 8] at [[0.5 0. 0. ]
[0.5 0.25 0. ]
[0.25 0.25 0. ]]
Parent 1 has vertices [0 2 3] at [[0.5 0. 0. ]
[1. 0.5 0. ]
[0.5 0.5 0. ]]
…

From the vertices we can see that parent 1 does not contain child 6.
Is this a bug or am I interpreting the results wrong? If it’s a bug, is there any other way I can get the mesh cells’ parent-child relationships, or can it maybe be fixed easily? And if it’s the latter, how should they be interpreted?

2

I think what you want is:

from mpi4py import MPI
import numpy as np
import dolfinx

msh_c = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)
msh_c.topology.create_entities(1)
msh_f, parent_cells, _ = dolfinx.mesh.refine_plaza(msh_c, option=dolfinx.mesh.RefinementOption.parent_cell)

child_vertices = dolfinx.mesh.entities_to_geometry(msh_f, 2, np.arange(len(parent_cells), dtype=np.int32))
parent_vertices = dolfinx.mesh.entities_to_geometry(msh_c, 2, parent_cells)


for i, (child_vertex, parent_vertex) in enumerate(zip(child_vertices, parent_vertices)):
    print(f"Child {i} has parent {parent_cells[i]}")
    child_vertex_coords = msh_f.geometry.x[child_vertex]
    print(f"Child {i} has vertices {child_vertex} at {child_vertex_coords}")
    parent_vertex_coords = msh_c.geometry.x[parent_vertex]
    print(f"Parent {parent_cells[i]} has vertices {parent_vertex} at {parent_vertex_coords}")

i.e. this computes the relationship between the mesh topology and geometry before accessing the coordinates.

3

Seems like I’m getting the same output as before with your code (I only left both nx and ny of create_unit_square(comm, nx, ny) at 2). The problem persists.

4

Agreed. Illustrated with:

from mpi4py import MPI
import numpy as np
import dolfinx

msh_c = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
msh_c.topology.create_entities(1)
msh_f, parent_cells, _ = dolfinx.mesh.refine_plaza(msh_c, redistribute=False, option=dolfinx.mesh.RefinementOption.parent_cell)

child_vertices = dolfinx.mesh.entities_to_geometry(msh_f, 2, np.arange(len(parent_cells), dtype=np.int32))
parent_vertices = dolfinx.mesh.entities_to_geometry(msh_c, 2, parent_cells)


for i, (child_vertex, parent_vertex) in enumerate(zip(child_vertices, parent_vertices)):
    print(f"Child {i} has parent {parent_cells[i]}")
    child_vertex_coords = msh_f.geometry.x[child_vertex]
    print(f"Child {i} has vertices {child_vertex} at {child_vertex_coords}")
    parent_vertex_coords = msh_c.geometry.x[parent_vertex]
    print(f"Parent {parent_cells[i]} has vertices {parent_vertex} at {parent_vertex_coords}")



num_cells = msh_c.topology.index_map(msh_c.topology.dim).size_local
marker = np.arange(num_cells, dtype=np.int32)
ct = dolfinx.mesh.meshtags(msh_c, msh_c.topology.dim, marker, marker)

ct_refined = dolfinx.mesh.meshtags(msh_f, msh_f.topology.dim, np.arange(len(parent_cells), dtype=np.int32), parent_cells)


with dolfinx.io.XDMFFile(msh_c.comm, "ct.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh_c)
    xdmf.write_meshtags(ct, msh_c.geometry)


with dolfinx.io.XDMFFile(msh_c.comm, "ct_refined.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh_f)
    xdmf.write_meshtags(ct_refined, msh_f.geometry)

yielding:

image
image1397×597 25.3 KB

Issue made at: Parent cells not correct in refine_plaza · Issue #3303 · FEniCS/dolfinx · GitHub

1 Like
5

Turns out you can get the mapping with the following:

from mpi4py import MPI
import numpy as np
import dolfinx

msh_c = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 5, ghost_mode=dolfinx.mesh.GhostMode.none)
msh_c.topology.create_entities(1)
msh_f, parent_cells, _ = dolfinx.mesh.refine_plaza(msh_c, redistribute=False, option=dolfinx.mesh.RefinementOption.parent_cell)

child_vertices = dolfinx.mesh.entities_to_geometry(msh_f, 2, np.arange(len(parent_cells), dtype=np.int32))
parent_vertices = dolfinx.mesh.entities_to_geometry(msh_c, 2, parent_cells)

num_cells = msh_c.topology.index_map(msh_c.topology.dim).size_local
marker = np.arange(num_cells, dtype=np.int32)
ct = dolfinx.mesh.meshtags(msh_c, msh_c.topology.dim, marker, marker)


ct_refined = dolfinx.mesh.MeshTags(dolfinx.cpp.refinement.transfer_cell_meshtag(ct._cpp_object, msh_f.topology, parent_cells))




with dolfinx.io.XDMFFile(msh_c.comm, "ct.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh_c)
    xdmf.write_meshtags(ct, msh_c.geometry)


with dolfinx.io.XDMFFile(msh_f.comm, "ct_refined.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh_f)
    xdmf.write_meshtags(ct_refined, msh_f.geometry)

where ct_refined.values is the map between cells in the refined mesh and cells in the original mesh.

1 Like