Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generalize computation of mesh::h to arbitrary order and entity #2162

Closed
wants to merge 15 commits into from

Conversation

jorgensd
Copy link
Sponsor Member

@jorgensd jorgensd commented Apr 28, 2022

Fix computation of h for arbitrary order geometries and sub entities (maximum distance between any vertex in the entity).

The current implementation throws a memory error on higher order geometries.

The new implementation works for all cell types, and any sub-dimension (topological dimension 0 returns 0).

The new implementation removes all copying and avoids xtensor

@garth-wells
Copy link
Member

I think the PR title and description could be improved. The PR title will appear in the ChangeLog.

PR description should explain what 'generalise' means. What is computed for higher order element?

@jorgensd jorgensd changed the title Generalize h Generalize computation of CellSize to arbitrary order Apr 30, 2022
@jorgensd
Copy link
Sponsor Member Author

I think the PR title and description could be improved. The PR title will appear in the ChangeLog.

PR description should explain what 'generalise' means. What is computed for higher order element?

Generalize means that it takes the maximum distance between any node in the cell, i.e. For higher order cells, the distance between any of the nodes on the facets and a vertex can be larger than the distance between any set of vertices.

@jorgensd
Copy link
Sponsor Member Author

I think the PR title and description could be improved. The PR title will appear in the ChangeLog.

PR description should explain what 'generalise' means. What is computed for higher order element?

@garth-wells I've now updated the description of the Pull request, and extended it.
Now dolfinx::mesh::h can compute the size of any mesh entity (vertex, edge, facet, cell) for any cell type.

@jorgensd jorgensd changed the title Generalize computation of CellSize to arbitrary order Generalize computation of mesh::h to arbitrary order and entity Apr 30, 2022
@garth-wells
Copy link
Member

I'm not sure about the definition of h. Is this a standard definition?

@jorgensd
Copy link
Sponsor Member Author

jorgensd commented May 3, 2022

I'm not sure about the definition of h. Is this a standard definition?

h is usually the CellDiameter, which would be the largest distance between any two points in a cell.
If the cell is sufficiently curved, simply computing the distance between the vertices is not going to be sufficient.

We currently do not support this through the form compiler (FFCx), and neither does TSFC, see: firedrakeproject/tsfc#224

@jorgensd
Copy link
Sponsor Member Author

jorgensd commented May 3, 2022

@garth-wells Code has been rewritten to only use the vertices of the entity.

cpp/dolfinx/mesh/utils.cpp Outdated Show resolved Hide resolved
cpp/dolfinx/mesh/utils.cpp Outdated Show resolved Hide resolved
{
auto p0 = xt::row(points, i);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't quite see the point of this change while e_vertices is still an xtensor type. It seems a but muddled being a mixture of xtensor std STL.

Copy link
Sponsor Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are no problems with taking row-views og xtensors, as they are contiguous in memory.
Eventually entities_to_geometry should return an adjacency-list, such that it can support prisms and pyramids, and that would remove the need for this view.

The main improvements in this PR is that the code now works for higher order geometries, sub entities of cells, and does not do any copying when computing h.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's confusing to be mixing concepts.

/// Compute greatest distance between any two vertices in the entity
/// @param[in] mesh The mesh
/// @param[in] entities List of entities (local to process)
/// @param[in] dim The topological dimension of the entity
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Needs a @return in doc.

/// Compute greatest distance between any two vertices
/// Compute greatest distance between any two vertices in the entity
/// @param[in] mesh The mesh
/// @param[in] entities List of entities (local to process)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

/// @param[in] entities List of entity indices to compute the size measure for.

@@ -64,7 +64,10 @@ graph::AdjacencyList<std::int64_t>
extract_topology(const CellType& cell_type, const fem::ElementDofLayout& layout,
const graph::AdjacencyList<std::int64_t>& cells);

/// Compute greatest distance between any two vertices
/// Compute greatest distance between any two vertices in the entity
/// @param[in] mesh The mesh
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

/// @param[in] mesh Mesh that entities belong to.

assert(np.allclose(h, np.sqrt(3 / (N**2))))


@pytest.mark.parametrize("ct", [CellType.hexahedron, CellType.tetrahedron])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ct is cryptic.

// Get number of cell vertices
const CellType type = cell_entity_type(mesh.topology().cell_type(), dim, 0);
const int num_vertices = num_cell_vertices(type);
// Get entity vertices in geometry dofs
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment doesn't make sense to me.

{
auto p0 = xt::row(points, i);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's confusing to be mixing concepts.

const xtl::span<const double> x_g = geometry.x();

// Lambda function to compute norm
auto l2_norm = [](auto p0, auto p1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The auto make things a bit unclear. The return type is always double, so it would be clearer if it was explicit.

I'm not sure what p0 and p1 are. There are index into, but from where the function is called it looks like an xtensor iterator.

@garth-wells
Copy link
Member

Merged into #2193.

@jorgensd jorgensd deleted the dokken/mesh_h_arbitrary_order branch November 10, 2023 15:53
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants