C++ Geometry object not constructible from Python

[idoby]

Hello everyone!

I was wondering why the Geometry_float32/float64 objects are no longer accessible from Python. Specifically, I'm trying to write code to refine a P2 mesh down to an equivalent P1 mesh for drawing and other operations.

This used to work in previous versions, but now seems not to be possible. It seems that a constructor for GF32/64 is not exposed to the Python side.

I'd appreciate your help very much, or if there is a better way to accomplish this, I'm also open to that (although I would like the flexibility since I have a few more uses in mind)

Thank you!

[francesco-ballarin]

It seems that a constructor for GF32/64 is not exposed to the Python side.

Looking at https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/wrappers/mesh.cpp#L154, indeed there does not seem to be a constructor.

This used to work in previous versions, but now seems not to be possible.

Can you tell us in which versions it used to work?

I am looking at https://github.com/FEniCS/dolfinx/blame/main/python/dolfinx/wrappers/mesh.cpp at the same line as above: I have gone back years, but I don't seem to be ever see a constructor being wrapped.

If you have a small snippet of code of what used to work, that would also be great.

[idoby]

Hi, thanks for replying

I think the older version I was using was 0.6.0, after upgrading to 0.7.2, it no longer seems to be possible to create geometry objects from Python. The objects may have been moved around/refactored/renamed, but I did have working refinement code whose output I was able to draw with matplotlib.

[idoby]

I see now that the fem.Mesh object was able to take a numpy array, and now it requires a Geometry object, which can't be created.

[francesco-ballarin]

From what I see, the last commit at which dolfinx.mesh.Mesh was not expecting a Geometry object was https://github.com/FEniCS/dolfinx/blame/e2e0d1a967acaa5a953e8dbf096770b4418e10a6/python/dolfinx/wrappers/mesh.cpp#L179, which is circa 2020, which is definitely not v0.6.0.

Afterwards, https://github.com/FEniCS/dolfinx/blame/c30e6e2f4a85fe3395c05118f3adb9b9257dc196/python/dolfinx/wrappers/mesh.cpp#L163, a Geometry object was required.

It would helpful if you could provide us with a small snippet of code of what used to work (as small as possible, not the entire refinement code).

[idoby]

Something along the lines of the following worked before I upgraded dolfinx:

new_connectivity = dofs[:, new_tri_dofs].reshape((-1, 3))
new_connectivity = fcpp.graph.AdjacencyList_int32(np.cast['int32'](new_connectivity))

new_index_map = fcpp.common.IndexMap(mesh.comm, 4 * index_map.size_local)
new_topology = fcpp.mesh.Topology(mesh.comm, fcpp.mesh.CellType.triangle)
new_topology.set_connectivity(new_connectivity, 2, 0)
new_topology.set_index_map(0, mesh.topology.index_map(0))
new_topology.set_index_map(1, mesh.topology.index_map(1))
new_topology.set_index_map(2, new_index_map)

new_coords = function_space.tabulate_dof_coordinates()
new_domain = ufl.Mesh(basix.ufl_wrapper.create_vector_element('Lagrange', 'triangle', 1))

return fmesh.Mesh(mesh.comm, new_topology, new_coords, new_domain)

where function_space is the P2 function space being refined, new_tri_dofs is an index array to resample the dofs to break up the triangles, etc

[jorgensd]

Usually meshes are created with the dolfinx.mesh.create_mesh function, as for instance shown in: https://jsdokken.com/FEniCS23-tutorial/src/mesh_generation.html#mesh-generation This is usually easier to work with from Python, as it works on numpy arrays and not DOLFINx objects.

However, it seems like your code is done something more custom, based on an existing mesh. It would be nice if the snippet could reproduce a use case (i.e. a code that we could run in 0.6.0 and get an output of, so we could understand the workflow).

[idoby]

Usually meshes are created with the dolfinx.mesh.create_mesh function, as for instance shown in: https://jsdokken.com/FEniCS23-tutorial/src/mesh_generation.html#mesh-generation This is usually easier to work with from Python, as it works on numpy arrays and not DOLFINx objects.

However, it seems like your code is done something more custom, based on an existing mesh. It would be nice if the snippet could reproduce a use case (i.e. a code that we could run in 0.6.0 and get an output of, so we could understand the workflow).

That did seem to work, thanks! Now the refined mesh seems to be correct, but the dof ordering seems to have changed. How is the dof ordering represented and how can I make sure it stays the same in the refined mesh? Basically, I'm trying to create a P1 mesh that is as close as possible to the original P2 mesh so I can draw it (I'm not really happy with pyvista/VTK's ability to generate publication-quality output).

Thank you for the willingness to help, I'll try to get some runnable code after my submission deadline. Would you be interested in accepting matplotlib-based plotting code into dolfinx?

[jorgensd]

Usually meshes are created with the dolfinx.mesh.create_mesh function, as for instance shown in: https://jsdokken.com/FEniCS23-tutorial/src/mesh_generation.html#mesh-generation This is usually easier to work with from Python, as it works on numpy arrays and not DOLFINx objects. However, it seems like your code is done something more custom, based on an existing mesh. It would be nice if the snippet could reproduce a use case (i.e. a code that we could run in 0.6.0 and get an output of, so we could understand the workflow).

That did seem to work, thanks! Now the refined mesh seems to be correct, but the dof ordering seems to have changed. How is the dof ordering represented and how can I make sure it stays the same in the refined mesh? Basically, I'm trying to create a P1 mesh that is as close as possible to the original P2 mesh so I can draw it (I'm not really happy with pyvista/VTK's ability to generate publication-quality output).

Thank you for the willingness to help, I'll try to get some runnable code after my submission deadline. Would you be interested in accepting matplotlib-based plotting code into dolfinx?

The main reason for us moving away from matplotlib plotting is that it is only suitable for 2D affine simplicies (triangles or tetrahedra with linear edges). It does not work for non-affine geometries, and in particular for unstructured quads/hexes.

We are happy for people to provide extensions (as a stand-alone library) that does plotting with matplotlib, for instance by posting to the FEniCS slack or discourse.

I guess in general, we should let the user to create a Topology and Geometry and pass it down into C++, or what do you think @garth-wells ?

[idoby]

Usually meshes are created with the dolfinx.mesh.create_mesh function, as for instance shown in: https://jsdokken.com/FEniCS23-tutorial/src/mesh_generation.html#mesh-generation This is usually easier to work with from Python, as it works on numpy arrays and not DOLFINx objects. However, it seems like your code is done something more custom, based on an existing mesh. It would be nice if the snippet could reproduce a use case (i.e. a code that we could run in 0.6.0 and get an output of, so we could understand the workflow).

That did seem to work, thanks! Now the refined mesh seems to be correct, but the dof ordering seems to have changed. How is the dof ordering represented and how can I make sure it stays the same in the refined mesh? Basically, I'm trying to create a P1 mesh that is as close as possible to the original P2 mesh so I can draw it (I'm not really happy with pyvista/VTK's ability to generate publication-quality output). Thank you for the willingness to help, I'll try to get some runnable code after my submission deadline. Would you be interested in accepting matplotlib-based plotting code into dolfinx?

The main reason for us moving away from matplotlib plotting is that it is only suitable for 2D affine simplicies (triangles or tetrahedra with linear edges). It does not work for non-affine geometries, and in particular for unstructured quads/hexes.

We are happy for people to provide extensions (as a stand-alone library) that does plotting with matplotlib, for instance by posting to the FEniCS slack or discourse.

I guess in general, we should let the user to create a Topology and Geometry and pass it down into C++, or what do you think @garth-wells ?

I can definitely see the problems with using MPL, and have evidently run into them myself, but currently dolfinx offers no easy way to create publication-quality vector imagery as PyVista, VTK and ParaView all target interactive rendering and will output only rasterized PDFs and SVGs.

Conversely, I believe most mesh types (even with higher-order interpolation and higher-order edges) can eventually be rendered into vector formats using MPL and a some coding work. For example, I have written code to render quadrilateral meshes. Has anyone, to your knowledge, published papers using dolfinx and has created vector-based figures for their papers? I'd like to ask them how they did it, if possible.

For example, I have created this figure for my paper. I don't know how else I could have accomplished it. I'm posting a rasterized image here, but MPL outputs true vector PDF or SVG figures just as easily.

Either way, I'd be happy if dolfinx could expose as much as possible to the Python side so I could manipulate the objects to create these figures.

[idoby]

Hi there!

Using @jorgensd 's advice I was able to make progress on this. However, it's still not working perfectly. Maybe you guys could offer some advice, and I would appreciate it.

Basically, I'm trying to take a P2 mesh, for example, this one with the P2 DoF locations marked in red:

I successfully refine it to an equivalent P1 mesh:

However, DOLFINx seems to be doing some magic behind the scenes and reordering the DoFs, resulting in a garbled plot:

This was supposed to plot the function lambda x: x[0], i.e., the x coordinate of the node. As you can see, it did quite a lot of reordering. Either that or my refinement code is missing some assumption about the ordering of the DoFs within an element or globally.

Here's the refinement code:

_P2_to_P1_TRI_REFINEMENT = ((0, 5, 4), (5, 3, 4),
                            (5, 1, 3), (4, 3, 2))

# Break up all triangles into 4 smaller triangles each
new_connectivity = dofs[:, _P2_to_P1_TRI_REFINEMENT].reshape((-1, 3))
new_connectivity = fcpp.graph.AdjacencyList_int64(np.cast['int64'](new_connectivity))
new_geometry = function_space.tabulate_dof_coordinates()[:, :2]

new_domain = ufl.Mesh(basix.ufl.element('Lagrange', 'triangle', 1, shape=(2,)))
return fmesh.create_mesh(mesh.comm, new_connectivity, new_geometry, new_domain)

Thanks for any advice!

[jorgensd]

Hi @idoby,

Sorry for the late reply, this fell out of my inbox. Please note that we have interpolation over non-matching meshes that you can use, see for instance: https://github.com/FEniCS/dolfinx/blob/7bb01adaa736a325ca048e37aec0b51cf9ff1f2e/python/test/unit/fem/test_interpolation.py#L913-L934 which I believe would resolve your problem.

There is currently some work on geometric multigrid by @schnellerhase that might help you as well.

Anyhow, do you see this issue as resolved or do you still want the C++ geometry?

Note that I've also covered how to map input indices (in mesh constructor) to dofs in: https://fenicsproject.discourse.group/t/how-to-apply-a-boundary-condition-to-a-range-of-dofs/15758/7

[schnellerhase]

The reordering of the vertices (and thus DOF's) is part of create_mesh. We faced the same problem in https://github.com/FEniCS/dolfinx/pull/3359 where we are trying to check the generation process of meshes, basically just like you trying to figure out which index goes where - we are working on a fix for this. Will let you know once this is resolved.

Additionally if you are executing this in parallel you will also see reordering happening due to the redistribution/-partitioning of the mesh in create_rectangle. The fenicsx own refinement structure has a custom repartitioner implemented to maintain the former ownership structure on the refined mesh, see https://github.com/FEniCS/dolfinx/blob/7bb01adaa736a325ca048e37aec0b51cf9ff1f2e/cpp/dolfinx/refinement/utils.h#L293.

[idoby]

Hi @jorgensd, @schnellerhase, thanks for replying.

I think the entire approach here is backwards. Trying to solve this specific problem by creating another big mechanism to do something specific, is the exact approach that led to both this issue and @schnellerhase's issue. This is because create_mesh does too much and assumes too much, and I can only assume the new function takes the same approach.

In my humble opinion, the approach should be to provide as many tools into the toolbox as possible. It's great to also provide high-level mechanisms to do common things but this should not be done in place of exposing all the underlying tools, i.e., whatever create_mesh does - I should be able to copy its code and do whatever customizations I need on top of the original code.

In summary: yes, unless there's a really good reason not to allow creation of a Geometry object and subsequently a Mesh object from Python, please expose the constructor and anything else that is not yet exposed to Python.

[idoby]

@jorgensd Thanks a lot!