diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 8e5288bb8e05b7d27a6ed68d0cb364a546ecf11c..ca5ddfb4f5d46f1b1a92c66983ca383129e572e1 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -673,7 +673,7 @@ def generate_urchin(order, m, n, est_rel_interp_tolerance, min_rad=0.2):
# {{{ generate_box_mesh
def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
- group_factory=None):
+ group_factory=None, boundary_tag_to_face=None):
"""Create a semi-structured mesh.
:param axis_coords: a tuple with a number of entries corresponding
@@ -681,12 +681,27 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
specifying the coordinates to be used along that axis.
:param group_factory: One of :class:`meshmode.mesh.SimplexElementGroup`
or :class:`meshmode.mesh.TensorProductElementGroup`.
+ :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
+ The keys correspond to custom boundary tags, with the values giving
+ a list of the faces on which they should be applied in terms of coordinate
+ directions (``+x``, ``-x``, ``+y``, ``-y``, ``+z``, ``-z``, ``+w``, ``-w``).
+
+ For example::
+
+ boundary_tag_to_face={"bdry_1": ["+x", "+y"], "bdry_2": ["-x"]}
.. versionchanged:: 2017.1
*group_factory* parameter added.
+
+ .. versionchanged:: 2020.1
+
+ *boundary_tag_to_face* parameter added.
"""
+ if boundary_tag_to_face is None:
+ boundary_tag_to_face = {}
+
for iaxis, axc in enumerate(axis_coords):
if len(axc) < 2:
raise ValueError("need at least two points along axis %d"
@@ -780,8 +795,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
el_vertices.append((a011, a111, a101, a110))
else:
- raise NotImplementedError("box meshes of dimension %d"
- % dim)
+ raise NotImplementedError("box meshes of dimension %d" % dim)
el_vertices = np.array(el_vertices, dtype=np.int32)
@@ -789,22 +803,85 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
vertices.reshape(dim, -1), el_vertices, order,
group_factory=group_factory)
+ # {{{ compute facial adjacency for mesh if there is tag information
+
+ facial_adjacency_groups = None
+ face_vertex_indices_to_tags = {}
+ boundary_tags = list(boundary_tag_to_face.keys())
+ axes = ["x", "y", "z", "w"]
+
+ if boundary_tags:
+ from pytools import indices_in_shape
+ vert_index_to_tuple = {
+ vertex_indices[itup]: itup
+ for itup in indices_in_shape(shape)}
+
+ for tag_idx, tag in enumerate(boundary_tags):
+ # Need to map the correct face vertices to the boundary tags
+ for face in boundary_tag_to_face[tag]:
+ if len(face) != 2:
+ raise ValueError("face identifier '%s' does not "
+ "consist of exactly two characters" % face)
+
+ side, axis = face
+ try:
+ axis = axes.index(axis)
+ except ValueError:
+ raise ValueError("unrecognized axis in face identifier '%s'" % face)
+ if axis >= dim:
+ raise ValueError("axis in face identifier '%s' does not exist in %dD"
+ % (face, dim))
+
+ if side == "-":
+ vert_crit = 0
+ elif side == "+":
+ vert_crit = shape[axis] - 1
+ else:
+ raise ValueError("first character of face identifier '%s' is not"
+ "'+' or '-'" % face)
+
+ for ielem in range(0, grp.nelements):
+ for ref_fvi in grp.face_vertex_indices():
+ fvi = grp.vertex_indices[ielem, ref_fvi]
+ fvi_tuples = [vert_index_to_tuple[i] for i in fvi]
+
+ if all(fvi_tuple[axis] == vert_crit for fvi_tuple in fvi_tuples):
+ key = frozenset(fvi)
+ face_vertex_indices_to_tags.setdefault(key, []).append(tag)
+
+ if boundary_tags:
+ from meshmode.mesh import (
+ _compute_facial_adjacency_from_vertices, BTAG_ALL, BTAG_REALLY_ALL)
+ boundary_tags.extend([BTAG_ALL, BTAG_REALLY_ALL])
+ facial_adjacency_groups = _compute_facial_adjacency_from_vertices(
+ [grp], boundary_tags, np.int32, np.int8,
+ face_vertex_indices_to_tags)
+ else:
+ facial_adjacency_groups = None
+
+ # }}}
+
from meshmode.mesh import Mesh
return Mesh(vertices, [grp],
- is_conforming=True)
+ facial_adjacency_groups=facial_adjacency_groups,
+ is_conforming=True, boundary_tags=boundary_tags)
# }}}
# {{{ generate_regular_rect_mesh
-def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
+def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
+ boundary_tag_to_face=None,
+ group_factory=None):
"""Create a semi-structured rectangular mesh.
:param a: the lower left hand point of the rectangle
:param b: the upper right hand point of the rectangle
:param n: a tuple of integers indicating the total number of points
on [a,b].
+ :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
+ See :func:`generate_box_mesh`.
"""
if min(n) < 2:
raise ValueError("need at least two points in each direction")
@@ -812,7 +889,9 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
axis_coords = [np.linspace(a_i, b_i, n_i)
for a_i, b_i, n_i in zip(a, b, n)]
- return generate_box_mesh(axis_coords, order=order)
+ return generate_box_mesh(axis_coords, order=order,
+ boundary_tag_to_face=boundary_tag_to_face,
+ group_factory=group_factory)
# }}}
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 77a25cd5dd4de07e1207ca137ec5b1510d268a84..ded91fe297a14c1d2241fb67486fbc0843fa237a 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -33,12 +33,13 @@ from pyopencl.tools import ( # noqa
pytest_generate_tests_for_pyopencl
as pytest_generate_tests)
+from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup
from meshmode.discretization.poly_element import (
InterpolatoryQuadratureSimplexGroupFactory,
PolynomialWarpAndBlendGroupFactory,
PolynomialEquidistantSimplexGroupFactory,
)
-from meshmode.mesh import BTAG_ALL
+from meshmode.mesh import Mesh, BTAG_ALL
from meshmode.discretization.connection import \
FACE_RESTR_ALL, FACE_RESTR_INTERIOR
import meshmode.mesh.generation as mgen
@@ -123,6 +124,84 @@ def test_boundary_tags():
# }}}
+# {{{ test custom boundary tags on box mesh
+
+@pytest.mark.parametrize(("dim", "nelem"), [
+ (1, 20),
+ (2, 20),
+ (3, 10),
+ ])
+@pytest.mark.parametrize("group_factory", [
+ SimplexElementGroup,
+
+ # FIXME: Not implemented: TPE.face_vertex_indices
+ # TensorProductElementGroup
+ ])
+def test_box_boundary_tags(dim, nelem, group_factory):
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ from meshmode.mesh import is_boundary_tag_empty
+ from meshmode.mesh import check_bc_coverage
+ if dim == 1:
+ a = (0,)
+ b = (1,)
+ n = (nelem,)
+ btag_to_face = {"btag_test_1": ["+x"],
+ "btag_test_2": ["-x"]}
+ elif dim == 2:
+ a = (0, -1)
+ b = (1, 1)
+ n = (nelem, nelem)
+ btag_to_face = {"btag_test_1": ["+x", "-y"],
+ "btag_test_2": ["+y", "-x"]}
+ elif dim == 3:
+ a = (0, -1, -1)
+ b = (1, 1, 1)
+ n = (nelem, nelem, nelem)
+ btag_to_face = {"btag_test_1": ["+x", "-y", "-z"],
+ "btag_test_2": ["+y", "-x", "+z"]}
+ mesh = generate_regular_rect_mesh(a=a, b=b,
+ n=n, order=3,
+ boundary_tag_to_face=btag_to_face,
+ group_factory=group_factory)
+ # correct answer
+ if dim == 1:
+ num_on_bdy = 1
+ else:
+ num_on_bdy = dim*(dim-1)*(nelem-1)**(dim-1)
+
+ assert not is_boundary_tag_empty(mesh, "btag_test_1")
+ assert not is_boundary_tag_empty(mesh, "btag_test_2")
+ check_bc_coverage(mesh, ['btag_test_1', 'btag_test_2'])
+
+ # check how many elements are marked on each boundary
+ num_marked_bdy_1 = 0
+ num_marked_bdy_2 = 0
+ btag_1_bit = mesh.boundary_tag_bit("btag_test_1")
+ btag_2_bit = mesh.boundary_tag_bit("btag_test_2")
+ for igrp in range(len(mesh.groups)):
+ bdry_fagrp = mesh.facial_adjacency_groups[igrp].get(None, None)
+
+ if bdry_fagrp is None:
+ continue
+
+ for i, nbrs in enumerate(bdry_fagrp.neighbors):
+ if (-nbrs) & btag_1_bit:
+ num_marked_bdy_1 += 1
+ if (-nbrs) & btag_2_bit:
+ num_marked_bdy_2 += 1
+
+ # raise errors if wrong number of elements marked
+ if num_marked_bdy_1 != num_on_bdy:
+ raise ValueError("%i marked on custom boundary 1, should be %i" %
+ (num_marked_bdy_1, num_on_bdy))
+ if num_marked_bdy_2 != num_on_bdy:
+ raise ValueError("%i marked on custom boundary 2, should be %i" %
+ (num_marked_bdy_2, num_on_bdy))
+
+
+# }}}
+
+
# {{{ convergence of boundary interpolation
@pytest.mark.parametrize("group_factory", [
@@ -527,7 +606,6 @@ def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False):
def test_merge_and_map(ctx_factory, visualize=False):
from meshmode.mesh.io import generate_gmsh, FileSource
from meshmode.mesh.generation import generate_box_mesh
- from meshmode.mesh import TensorProductElementGroup
from meshmode.discretization.poly_element import (
PolynomialWarpAndBlendGroupFactory,
LegendreGaussLobattoTensorProductGroupFactory)
@@ -1009,7 +1087,6 @@ def no_test_quad_mesh_3d():
def test_quad_single_element():
from meshmode.mesh.generation import make_group_from_vertices
- from meshmode.mesh import Mesh, TensorProductElementGroup
vertices = np.array([
[0.91, 1.10],
@@ -1037,7 +1114,6 @@ def test_quad_single_element():
def test_quad_multi_element():
from meshmode.mesh.generation import generate_box_mesh
- from meshmode.mesh import TensorProductElementGroup
mesh = generate_box_mesh(
(
np.linspace(3, 8, 4),