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),