Skip to content

Commit

Permalink
Let entities be a list of tuples and add tests for on_boundary = True…
Browse files Browse the repository at this point in the history
… and False
  • Loading branch information
finsberg committed Sep 24, 2024
1 parent a5b5d39 commit 18e029b
Show file tree
Hide file tree
Showing 3 changed files with 78 additions and 37 deletions.
18 changes: 4 additions & 14 deletions examples/mark_entities.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,15 +38,12 @@ def inner(x):

from scifem import create_meshtags

cell_tag = create_meshtags(
mesh,
mesh.topology.dim,
[{"tag": 1, "locator": left}, {"tag": 3, "locator": right}, {"tag": 7, "locator": inner}],
)
cell_tag = create_meshtags(mesh, mesh.topology.dim, [(1, left), (3, right), (7, inner)])

# Next we can plot these marked entities

import pyvista

pyvista.start_xvfb()
vtk_grid = dolfinx.plot.vtk_mesh(mesh, cell_tag.dim, cell_tag.indices)
grid = pyvista.UnstructuredGrid(*vtk_grid)
Expand All @@ -72,9 +69,7 @@ def top(x):
return x[1] > 0.9


facet_tags = create_meshtags(
mesh, mesh.topology.dim - 1, [{"tag": 2, "locator": top}, {"tag": 7, "locator": circle}]
)
facet_tags = create_meshtags(mesh, mesh.topology.dim - 1, [(2, top), (7, circle)])


facet_grid = dolfinx.plot.vtk_mesh(mesh, facet_tags.dim, facet_tags.indices)
Expand All @@ -91,12 +86,7 @@ def top(x):
# We can also exclude interior facets by adding `on_boundary: True` (by default this is set to False).

boundary_facet_tags = create_meshtags(
mesh,
mesh.topology.dim - 1,
[
{"tag": 2, "locator": top, "on_boundary": True},
{"tag": 7, "locator": circle, "on_boundary": False},
],
mesh, mesh.topology.dim - 1, [(2, top, True), (7, circle, False)]
)


Expand Down
28 changes: 12 additions & 16 deletions src/scifem/mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,11 @@
__all__ = ["create_meshtags"]


# Workaround to make id and locator required but on_boundary optional
# see https://peps.python.org/pep-0655/
class _TaggedEntities(typing.TypedDict):
tag: int
locator: typing.Callable[[npt.NDArray[np.floating]], npt.NDArray[np.bool_]]


class TaggedEntities(_TaggedEntities, total=False):
on_boundary: bool
# (tag, locator, on_boundary) where on_boundary is optional
TaggedEntities = (
tuple[int, typing.Callable[[npt.NDArray[np.floating]], npt.NDArray[np.bool_]]]
| tuple[int, typing.Callable[[npt.NDArray[np.floating]], npt.NDArray[np.bool_]], bool]
)


def create_meshtags(
Expand All @@ -27,12 +23,12 @@ def create_meshtags(
Args:
domain: A ``dolfinx.mesh.Mesh`` object
dim: Dimension of the entities to mark
entities_list: A list of dictionaries with the following keys
entities_list: A list of tuples with the following elements:
- ``tag``: The tag to assign to the entities
- ``locator``: A function that takes a point and returns a boolean array
- ``index 0``: The tag to assign to the entities
- ``index 1``: A function that takes a point and returns a boolean array
indicating whether the point is inside the entity
- ``on_boundary``: Optional, if True, the entities will be marked on the boundary
- ``index 2``: Optional, if True, the entities will be marked on the boundary
Returns:
A ``dolfinx.mesh.MeshTags`` object with the corresponding entities marked.
Expand Down Expand Up @@ -60,9 +56,9 @@ def create_meshtags(

# Concatenate and sort the arrays based on indices
for tagged_entity in entities_list:
on_boundary = tagged_entity.get("on_boundary", False)
entities = locate_entities(on_boundary)(domain, dim, tagged_entity["locator"])
markers[entities] = tagged_entity["tag"]
on_boundary = False if len(tagged_entity) == 2 else tagged_entity[2]
entities = locate_entities(on_boundary)(domain, dim, tagged_entity[1])
markers[entities] = tagged_entity[0]

facets = np.flatnonzero(markers != -1).astype(np.int32)
return dolfinx.mesh.meshtags(domain, dim, facets, markers[facets])
69 changes: 62 additions & 7 deletions tests/test_mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,18 +9,18 @@
@pytest.mark.parametrize(
"entities_list, set_values",
[
([{"tag": 1, "locator": lambda x: x[0] <= 1.0}], {1}),
([{"tag": 2, "locator": lambda x: x[0] <= 1.0}], {2}),
([(1, lambda x: x[0] <= 1.0)], {1}),
([(2, lambda x: x[0] <= 1.0)], {2}),
(
[
{"tag": 1, "locator": lambda x: x[0] <= 1.0},
{"tag": 2, "locator": lambda x: x[0] >= 0.0},
(1, lambda x: x[0] <= 1.0),
(2, lambda x: x[0] >= 0.0),
],
{2},
),
],
)
def test_create_meshtags_celltags_all(entities_list, set_values):
def test_create_celltags_celltags_all(entities_list, set_values):
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 3, 3)
cell_tag = scifem.create_meshtags(mesh, mesh.topology.dim, entities_list=entities_list)

Expand All @@ -32,14 +32,69 @@ def test_create_meshtags_celltags_all(entities_list, set_values):
assert set(cell_tag.values) == set_values


@pytest.mark.parametrize(
"entities_list, set_values",
[
([(1, lambda x: x[0] <= 1.0, False)], {1}),
([(2, lambda x: x[0] <= 1.0, False)], {2}),
(
[
(1, lambda x: x[0] <= 1.0, False),
(2, lambda x: x[0] >= 0.0, False),
],
{2},
),
],
)
def test_create_facet_tags_all_on_boundary_False(entities_list, set_values):
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 3, 3)
facet_tag = scifem.create_meshtags(mesh, mesh.topology.dim - 1, entities_list=entities_list)
im = mesh.topology.index_map(mesh.topology.dim - 1)

assert facet_tag.dim == mesh.topology.dim - 1
assert facet_tag.indices.shape[0] == im.size_local + im.num_ghosts
assert facet_tag.values.shape[0] == im.size_local + im.num_ghosts
assert facet_tag.values.dtype == np.int32
assert set(facet_tag.values) == set_values


@pytest.mark.parametrize(
"entities_list, set_values",
[
([(1, lambda x: x[0] <= 1.0, True)], {1}),
([(2, lambda x: x[0] <= 1.0, True)], {2}),
(
[
(1, lambda x: x[0] <= 1.0, True),
(2, lambda x: x[0] >= 0.0, True),
],
{2},
),
],
)
def test_create_facet_tags_all_on_boundary_True(entities_list, set_values):
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 3, 3)
facet_tag = scifem.create_meshtags(mesh, mesh.topology.dim - 1, entities_list=entities_list)

mesh.topology.create_entities(1)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
facet_indices = dolfinx.mesh.exterior_facet_indices(mesh.topology)
assert facet_tag.dim == mesh.topology.dim - 1
assert facet_tag.indices.shape[0] == len(facet_indices)
assert facet_tag.values.shape[0] == len(facet_indices)
assert facet_tag.values.dtype == np.int32
assert set(facet_tag.values) == set_values


@pytest.mark.parametrize(
"entities_list",
[
([{"tag": 1, "locator": lambda x: x[0] < 0.0}]),
([(1, lambda x: x[0] < 0.0)]),
([(1, lambda x: x[0] < 0.0)]),
([]),
],
)
def test_create_meshtags_celltags_empty(entities_list):
def test_create_celltags_empty(entities_list):
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 3, 3)
cell_tag = scifem.create_meshtags(mesh, mesh.topology.dim, entities_list=entities_list)

Expand Down

0 comments on commit 18e029b

Please sign in to comment.