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

Wrap entities_to_geometry in Python #3038

Merged
merged 15 commits into from
Feb 20, 2024
8 changes: 4 additions & 4 deletions cpp/dolfinx/mesh/utils.h
Original file line number Diff line number Diff line change
Expand Up @@ -618,12 +618,12 @@ std::vector<std::int32_t> locate_entities_boundary(const Mesh<T>& mesh, int dim,
/// @warning This function should not be used unless there is no
/// alternative. It may be removed in the future.
///
/// @param[in] mesh The mesh
/// @param[in] dim Topological dimension of the entities of interest
/// @param[in] mesh The mesh.
/// @param[in] dim Topological dimension of the entities of interest.
/// @param[in] entities Entity indices (local) to compute the vertex
/// geometry indices for
/// geometry indices for.
/// @param[in] orient If true, in 3D, reorients facets to have
/// consistent normal direction
/// consistent normal direction.
/// @return Indices in the geometry array for the entity vertices. The
/// shape is `(num_entities, num_vertices_per_entity)` and the storage
/// is row-major. The index `indices[i, j]` is the position in the
Expand Down
25 changes: 24 additions & 1 deletion python/dolfinx/mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
"exterior_facet_indices", "compute_incident_entities", "create_cell_partitioner",
"create_interval", "create_unit_interval", "create_rectangle", "create_unit_square",
"create_box", "create_unit_cube", "to_type", "to_string", "refine_plaza",
"transfer_meshtag"]
"transfer_meshtag", "entities_to_geometry"]


class Mesh:
Expand Down Expand Up @@ -645,3 +645,26 @@ def create_unit_cube(comm: _MPI.Comm, nx: int, ny: int, nz: int, cell_type=CellT
"""
return create_box(comm, [np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0])],
[nx, ny, nz], cell_type, dtype, ghost_mode, partitioner)


def entities_to_geometry(
mesh: Mesh, dim: int, entities: npt.NDArray[np.int32], orient: bool
mscroggs marked this conversation as resolved.
Show resolved Hide resolved
) -> npt.NDArray[np.int32]:
"""Indices in the geometry data for each vertex of the given mesh entities.

Warning:
This function should not be used unless there is no alternative.
It may be removed in the future.

Args:
mesh: The mesh.
dim: Topological dimension of the entities of interest.
entities: Entity indices (local to the process) to determine the
vertex geometry indices for.
orient: If ``True``, in 3D, re-orients facets to have consistent
Copy link
Member

Choose a reason for hiding this comment

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

What does this mean? (I know it straight from the C++ docs). @jorgensd?

Copy link
Sponsor Member

Choose a reason for hiding this comment

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

This is a flag only used by BEMPP (at least at its introduction) to ensure consistent orientation of each facet (to make all cell normals point outwards).

@mscroggs might be able to comment on the necessity of this at the current stage of BEMPP

Copy link
Member

Choose a reason for hiding this comment

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

It it's still needed it needs an explanation of what is required for a consistent orientation, i.e. document how we related the normal vector to a facet to the vertex ordering. And does it apply for all facet shapes?

Copy link
Member Author

Choose a reason for hiding this comment

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

This is still used by Bempp-cl, so I'm keen to keep it for now. I think we've only tested it for triangle facets, and I'd be keen to work towards removing this once Bempp-rs can do coupling without needing this.

normal direction.

Returns:
Indices in the geometry array for the entity vertices.
"""
return _cpp.mesh.entities_to_geometry(mesh._cpp_object, dim, entities, orient)
6 changes: 3 additions & 3 deletions python/test/unit/mesh/test_mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,11 @@
from dolfinx import cpp as _cpp
from dolfinx import graph
from dolfinx import mesh as _mesh
from dolfinx.cpp.mesh import create_cell_partitioner, entities_to_geometry, is_simplex
from dolfinx.cpp.mesh import create_cell_partitioner, is_simplex
from dolfinx.fem import assemble_scalar, coordinate_element, form
from dolfinx.mesh import (CellType, DiagonalType, GhostMode, create_box, create_interval, create_rectangle,
create_submesh, create_unit_cube, create_unit_interval, create_unit_square,
exterior_facet_indices, locate_entities, locate_entities_boundary)
entities_to_geometry, exterior_facet_indices, locate_entities, locate_entities_boundary)


def submesh_topology_test(mesh, submesh, entity_map, vertex_map, entity_dim):
Expand Down Expand Up @@ -61,7 +61,7 @@ def submesh_geometry_test(mesh, submesh, entity_map, geom_map, entity_dim):
if len(entity_map) > 0:
assert mesh.geometry.dim == submesh.geometry.dim

e_to_g = entities_to_geometry(mesh._cpp_object, entity_dim, np.array(entity_map), False)
e_to_g = entities_to_geometry(mesh, entity_dim, np.array(entity_map), False)
for submesh_entity in range(len(entity_map)):
submesh_x_dofs = submesh.geometry.dofmap[submesh_entity]
# e_to_g[i] gets the mesh x_dofs of entities[i], which should
Expand Down
Loading