Dolfinx: Bug in mesh refinement and question about parent-children relationship

mesh
1

Hello,
I want to refine a 2D mesh on dolfinx. But after the refinement, the geometry is always 3D.

points = np.array([[0, 0],  [1.2, 0],  [1.5, 1],  [0.2, 1]])
cells = np.array([[0, 1, 2],  [0, 2, 3]])
cell = ufl.Cell("triangle", geometric_dimension=2)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, 1))
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, points, domain)
print(mesh.geometry.x.shape)  #[4,**3**]
print(mesh.topology.dim, mesh.geometry.dim)  # 2   2

mesh2=refine(mesh)
print(mesh2.topology.dim, mesh2.geometry.dim)  # 2 **3**

cells2 = mesh2.geometry.dofmap.array().reshape((-1, mesh.topology.dim + 1))
xy=np.delete(mesh2.geometry.x,2,1)  #remove z-axis
mesh3=dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells2, xy, domain)
print(mesh3.topology.dim, mesh3.geometry.dim) # 2 2

#dolfinx.plotting.plot(mesh3)

A trick is to create a mesh3 to make the dimension right. Is there a Pythonic way to do this?

Question2: I am trying to build the relationship between the original mesh and the refined mesh.

Create a matrix such that
each row contains the sub_cell_id’s of the refined mesh that are refined from a cell of the original mesh.

Thank you.

2

To address your first question, this is a bug, and thanks for reporting it.
I have created a suggested fix in dolfinx, which hopefully will be merged soon.

For question 2: I am not aware of any built-in functionality for this, but I might be mistaken there.
Also, you need to be more specific about you want this map to do, should it work for only CG 1 function spaces or arbitrary function spaces?
You also need to be very specific if this map should work in parallel or not, and if you want to redistribute the refined mesh or not.

3

Thanks, dokken.

For Q2, the map is simple one, for example
mesh refmesh

The map matrix is
| 0 | 1 | 2 | 3 |
| 4 | 5 | 6 | 7 |

Thank you.

4

select_colliding_cells is good enough. The middle point of each sub_cell collides with a cell of the original mesh. Perhaps this is not efficient but it works.

from dolfinx.cpp.geometry import select_colliding_cells
NE=cells.shape[0] # cell count of mesh
cells3 = mesh3.geometry.dofmap.array().reshape((-1, mesh.topology.dim + 1))
xy3=mesh3.geometry.x
mx3=np.sum(xy3[cells3,:], axis=1)/3 #middle points of refined mesh 
inx=np.squeeze(np.array([select_colliding_cells(mesh, range(NE), mx3[mi,:],1) 
                                          for mi in range(mx3.shape[0])]))
refine_map =np.array([np.where(inx==i)[0] for i in range(NE)])
print(refine_map)
5

The fix has been merged, so if you update your version of dolfinx, you should be able to get the correct gdim when refining.
Note that you can use

dolfinx.cpp.mesh.midpoints

to get the midpoints of your cells. See the following link for syntax.