Extend submesh to higher order geometries #2233
Conversation
…nx into dokken/higher-order-submesh
|
My assembly tests all pass with these changes |
| { | ||
| const std::int32_t idx = submesh_to_mesh_entity_map[i]; | ||
| assert(!e_to_c->links(idx).empty()); | ||
| // Always pick the last cell to be consistent with the e_to_v connectivity |
There was a problem hiding this comment.
Is this guaranteed? Or is it exploiting implicit behaviour that could change in the future?
There was a problem hiding this comment.
This is how the implementation of entity_to_vertex is currently done. If we change how we create those connectivities, this would change how the local ordering of vertices is computed.
For exterior facet integrals, this has no effect as there is only one facet.
For interior facet integrals, one facet needs a permutation, to align the vertices of the child facet with the vertices of the parent cell.
The old implementation gives the correct vertices, but in a permuted order. This implementation makes sure we get the entities in the consistent order defined in basix.
|
|
||
| // Submesh geometry input_global_indices | ||
| // TODO Check this | ||
| const std::vector<std::int64_t>& mesh_igi | ||
| = mesh.geometry().input_global_indices(); | ||
| const std::vector<std::int64_t>& mesh_igi = geometry.input_global_indices(); |
There was a problem hiding this comment.
Maybe a more instructive name than mesh_igi?
There was a problem hiding this comment.
Could we do that in a separate PR, as I haven't actually changed anything wrt. mesh_igi, except of some formatting.
| @@ -470,7 +470,8 @@ mesh::entities_to_geometry(const Mesh& mesh, int dim, | |||
| for (std::size_t i = 0; i < entities.size(); ++i) | |||
| { | |||
| const std::int32_t idx = entities[i]; | |||
| const std::int32_t cell = e_to_c->links(idx).front(); | |||
| // Always pick the second cell to be consistent with the e_to_v connectivity | |||
There was a problem hiding this comment.
See earlier comment on ordering.
| points += [[i / order, j / order, 0] for j in range(order + 1) | ||
| for i in range(order + 1 - j)] | ||
| for k in range(1, order): | ||
| points += [[i / order, j / order + 0.1, k / order] for j in range(order + 1 - k) |
There was a problem hiding this comment.
Removing line break would improve readability.
| assert np.isclose(area, volume) | ||
|
|
||
|
|
||
| @pytest.mark.skip_in_parallel | ||
| @pytest.mark.parametrize('order', range(1, 5)) | ||
| def test_submesh(order): |
There was a problem hiding this comment.
This test is hard to follow. Explanation on what is being done would be helpful.
There was a problem hiding this comment.
The mesh generator is the same as we use everywhere else in the test_higher_order_mesh.py.
Integration entities are unrolled in Python to be 1D (cell integrals), 2D (Exterior facet), 3D (Interior facet). Preparing code for higher order function spaces. Currently does not work due to a bug in higher order submeshes, see: FEniCS/dolfinx#2233
* Flatten integration entities in C++. Integration entities are unrolled in Python to be 1D (cell integrals), 2D (Exterior facet), 3D (Interior facet). Preparing code for higher order function spaces. Currently does not work due to a bug in higher order submeshes, see: FEniCS/dolfinx#2233 * Various copy improvements
…nx into dokken/higher-order-submesh
|
@garth-wells Can we merge this now? It has been thoroughly tested by @jpdean and adds important functionality to the class |
…nx into dokken/higher-order-submesh
Current code does not actually work for higher order geometries (and does not throw an error).
This PR extends submeshes to higher order geometries (edges,facets,cell based submeshes).