Create mixed topology mesh

[chrisrichardson]

Add a new experimental create_mesh() interface, which accepts mixed topology input, calls the partitioner and creates a mesh, in the same way as for single cell type.

[jorgensd]

I can't get this branch to work with mixing P1 and P2 triangles in the same mesh (even if they are completely disconnected.

Here is an image of my minimal mesh: If i change the order of the second triangle to 1, the code runs smoothly. I do not have minimal working example. But I've attached the mesh I've used here, and the code I've used for reading the mesh is available at: https://gist.github.com/jorgensd/fcbcee61f773e8889eccacf2f75158a2 (needs mesh_converter installed: https://github.com/jorgensd/mesh_converter/ test.vtkhdf.zip

[jorgensd]

Error message is:

python3: /root/shared/dolfinx/cpp/dolfinx/fem/dofmapbuilder.cpp:299: std::tuple<std::vector<{anonymous}::dofmap_t, std::allocator<{anonymous}::dofmap_t> >, std::vector<long int, std::allocator<long int> >, std::vector<std::pair<signed char, int>, std::allocator<std::pair<signed char, int> > >, std::vector<std::shared_ptr<const dolfinx::common::IndexMap>, std::allocator<std::shared_ptr<const dolfinx::common::IndexMap> > >, long int> {anonymous}::build_basic_dofmaps(const dolfinx::mesh::Topology&, const std::vector<dolfinx::fem::ElementDofLayout>&): Assertion `(int)num_entity_dofs == num_entity_dofs_et[k]' failed.
Aborted (core dumped)

[jorgensd]

Another thing I realized when looking through this is that we need to think about how to wrap this as a Python object, as the C++ object is not really usable in Python.

[chrisrichardson]

It will not currently work with mixed order meshes - the geometry dofmap would have hanging nodes... and I guess this causes a problem, even if the meshes are disconnected. I am sure this can be fixed in a future PR (e.g. by using DG geometry).

[chrisrichardson]

Another thing I realized when looking through this is that we need to think about how to wrap this as a Python object, as the C++ object is not really usable in Python.

It's too experimental at the moment. Accessing through dolfinx.cpp should be fine for developers. We still have other issues to resolve: e.g. FunctionSpace and assembly... https://github.com/orgs/FEniCS/projects/7

[jorgensd]

It will not currently work with mixed order meshes - the geometry dofmap would have hanging nodes... and I guess this causes a problem, even if the meshes are disconnected. I am sure this can be fixed in a future PR (e.g. by using DG geometry).

I don't understand why the dofmap would have hanging nodes? Shouldn't there be a dofmap per coordinate element type? Hanging nodes would only happen if a higher order cell shared two vertices with a lower order (linear) cell?

[chrisrichardson]

The dofmap construction needs the same number of dofs en each entity type (edges, facets etc.) or it will currently fail, even if the mesh is disjoint. Whilst each element has its own dofmap, they have to match if they are sharing degrees of freedom between cells. I think we can fix it - either by using DG geometry, or by updating the dofmap constructor. It's a future challenge.

[jorgensd]

The dofmap construction needs the same number of dofs en each entity type (edges, facets etc.) or it will currently fail, even if the mesh is disjoint. Whilst each element has its own dofmap, they have to match if they are sharing degrees of freedom between cells. I think we can fix it - either by using DG geometry, or by updating the dofmap constructor. It's a future challenge.

So that means that we are limited to linear geometries for now? As if we did some second order tetrahedra, prism and hexahedra the sub entities would have different dofs on them (i.e. a facet of a tetrahedra would have less dofs than a facet of a hexahedra?

Now it only works since the only place with dofs associated with them are the vertices, and each vertex has one dof regardless of the shape?

[chrisrichardson]

So that means that we are limited to linear geometries for now?

No, just compatible ones. If everything is 3rd order, in 3D for example, the dofs will agree on edges and facets (2 dofs on each edge, 1 dof on a triangle, 4 on a quad). This would actually be a good one to test, since the permutations have not been tested yet for mixed topology.