From 3d16443e10490ea43193b40941a5beb478ad34c3 Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Fri, 8 May 2020 16:02:02 -0500
Subject: [PATCH 01/10] Small changes to get SBP-DG cases running with Grudge

---
 meshmode/discretization/__init__.py | 3 +--
 meshmode/mesh/generation.py         | 9 +++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index 6a4fcda7..70a93798 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -356,9 +356,8 @@ class Discretization(object):
                     "stride:auto,stride:auto,stride:auto")
             return knl
 
-        result = self.empty(dtype=self.real_dtype, extra_dims=(self.ambient_dim,))
-
         with cl.CommandQueue(self.cl_context) as queue:
+            result = self.empty(queue=queue, dtype=self.real_dtype, extra_dims=(self.ambient_dim,))
             for grp in self.groups:
                 if grp.nelements == 0:
                     continue
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 3540dadd..de5d9581 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -658,7 +658,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_tags=[]):
     """Create a semi-structured mesh.
 
     :param axis_coords: a tuple with a number of entries corresponding
@@ -776,14 +776,15 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 
     from meshmode.mesh import Mesh
     return Mesh(vertices, [grp],
-            is_conforming=True)
+            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_tags=[]):
     """Create a semi-structured rectangular mesh.
 
     :param a: the lower left hand point of the rectangle
@@ -797,7 +798,7 @@ 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_tags=boundary_tags)
 
 # }}}
 
-- 
GitLab


From fc8cbe68b841c7c05673816adad421cd79170701 Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Fri, 8 May 2020 16:38:44 -0500
Subject: [PATCH 02/10] Placate Flake8

---
 meshmode/discretization/__init__.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index 70a93798..9340a982 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -357,7 +357,8 @@ class Discretization(object):
             return knl
 
         with cl.CommandQueue(self.cl_context) as queue:
-            result = self.empty(queue=queue, dtype=self.real_dtype, extra_dims=(self.ambient_dim,))
+            result = self.empty(queue=queue, dtype=self.real_dtype,
+                                extra_dims=(self.ambient_dim,))
             for grp in self.groups:
                 if grp.nelements == 0:
                     continue
-- 
GitLab


From 5d563ef7974c92d1060d22e4d86d1df6488e9c03 Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Tue, 12 May 2020 09:49:19 -0500
Subject: [PATCH 03/10] Reverts incorrect change to nodes result array

---
 meshmode/discretization/__init__.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index 9340a982..6a4fcda7 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -356,9 +356,9 @@ class Discretization(object):
                     "stride:auto,stride:auto,stride:auto")
             return knl
 
+        result = self.empty(dtype=self.real_dtype, extra_dims=(self.ambient_dim,))
+
         with cl.CommandQueue(self.cl_context) as queue:
-            result = self.empty(queue=queue, dtype=self.real_dtype,
-                                extra_dims=(self.ambient_dim,))
             for grp in self.groups:
                 if grp.nelements == 0:
                     continue
-- 
GitLab


From ce9852dd24e5ef61566af08aa48fe09b2acac1d0 Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Tue, 12 May 2020 16:42:23 -0500
Subject: [PATCH 04/10] Beefs up boundary changes for box mesh generation

---
 meshmode/mesh/__init__.py   |  3 +-
 meshmode/mesh/generation.py | 56 +++++++++++++++++++++++++++++++++++--
 2 files changed, 55 insertions(+), 4 deletions(-)

diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 6d6191aa..16983a06 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -1141,7 +1141,8 @@ def _compute_facial_adjacency_from_vertices(groups, boundary_tags,
                     tag_mask = 0
                     for tag in tags:
                         tag_mask |= boundary_tag_bit(tag)
-                    fagrp.neighbors[idx] = -(-(fagrp.neighbors[idx]) | tag_mask)
+                    fagrp.neighbors[idx] = -(tag_mask | boundary_tag_bit(
+                                             BTAG_REALLY_ALL))
 
         else:
             raise RuntimeError("unexpected number of adjacent faces")
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index de5d9581..04b5288f 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -658,7 +658,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, boundary_tags=[]):
+        group_factory=None, boundary_dict={}):
     """Create a semi-structured mesh.
 
     :param axis_coords: a tuple with a number of entries corresponding
@@ -666,6 +666,10 @@ 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_dict: an optional dictionary for boundary configuration.
+        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 cardinal
+        directions.
 
     .. versionchanged:: 2017.1
 
@@ -774,8 +778,53 @@ 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_dict.keys())
+    if boundary_tags:
+        for tag_idx, tag in enumerate(boundary_tags):
+            # Need to map the correct face vertices to the boundary tags
+            for face in boundary_dict[tag]:
+                if face == "east":
+                    face_id = 1
+                    dim_crit = 0
+                    node_crit = axis_coords[0][0]
+                elif face == "west":
+                    face_id = 1
+                    dim_crit = 0
+                    node_crit = axis_coords[0][-1]
+                elif face == "north":
+                    face_id = 0
+                    dim_crit = 1
+                    node_crit = axis_coords[1][-1]
+                elif face == "south":
+                    face_id = 0
+                    dim_crit = 1
+                    node_crit = axis_coords[1][0]
+                for ielem in range(0, grp.nelements):
+                    if grp.nodes[dim_crit][ielem][0] == node_crit:
+                        fvi = grp.vertex_indices[ielem,
+                                                 grp.face_vertex_indices()[
+                                                     face_id]]
+                        if frozenset(fvi) not in face_vertex_indices_to_tags:
+                            face_vertex_indices_to_tags[frozenset(fvi)] = []
+                        face_vertex_indices_to_tags[frozenset(fvi)].append(tag)
+        # Need to add BTAG ALL and BTAG_REALLY_ALL to list of tags
+        from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL
+        boundary_tags.append(BTAG_ALL)
+        boundary_tags.append(BTAG_REALLY_ALL)
+        from meshmode.mesh import _compute_facial_adjacency_from_vertices
+        facial_adjacency_groups = _compute_facial_adjacency_from_vertices(
+                [grp], boundary_tags, np.int32, np.int8,
+                face_vertex_indices_to_tags)
+    else:
+        facial_adjacency_groups = _compute_facial_adjacency_from_vertices(
+                [grp], boundary_tags, np.int32, np.int8)
+
     from meshmode.mesh import Mesh
     return Mesh(vertices, [grp],
+            facial_adjacency_groups=facial_adjacency_groups,
             is_conforming=True, boundary_tags=boundary_tags)
 
 # }}}
@@ -784,7 +833,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 # {{{ generate_regular_rect_mesh
 
 def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
-                               boundary_tags=[]):
+                               boundary_dict={}):
     """Create a semi-structured rectangular mesh.
 
     :param a: the lower left hand point of the rectangle
@@ -798,7 +847,8 @@ 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, boundary_tags=boundary_tags)
+    return generate_box_mesh(axis_coords, order=order,
+                             boundary_dict=boundary_dict)
 
 # }}}
 
-- 
GitLab


From fbca9d3033decf10e5c9ac17927335cc3e14b096 Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Wed, 13 May 2020 11:47:27 -0500
Subject: [PATCH 05/10] Get Meshmode tests to pass by reinstating BTAG_ALL
 everywhere

---
 meshmode/mesh/__init__.py   | 3 +--
 meshmode/mesh/generation.py | 5 +----
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 16983a06..6d6191aa 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -1141,8 +1141,7 @@ def _compute_facial_adjacency_from_vertices(groups, boundary_tags,
                     tag_mask = 0
                     for tag in tags:
                         tag_mask |= boundary_tag_bit(tag)
-                    fagrp.neighbors[idx] = -(tag_mask | boundary_tag_bit(
-                                             BTAG_REALLY_ALL))
+                    fagrp.neighbors[idx] = -(-(fagrp.neighbors[idx]) | tag_mask)
 
         else:
             raise RuntimeError("unexpected number of adjacent faces")
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 04b5288f..2376cb44 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -782,6 +782,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
     facial_adjacency_groups = None
     face_vertex_indices_to_tags = {}
     boundary_tags = list(boundary_dict.keys())
+    from meshmode.mesh import _compute_facial_adjacency_from_vertices
     if boundary_tags:
         for tag_idx, tag in enumerate(boundary_tags):
             # Need to map the correct face vertices to the boundary tags
@@ -814,13 +815,9 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
         from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL
         boundary_tags.append(BTAG_ALL)
         boundary_tags.append(BTAG_REALLY_ALL)
-        from meshmode.mesh import _compute_facial_adjacency_from_vertices
         facial_adjacency_groups = _compute_facial_adjacency_from_vertices(
                 [grp], boundary_tags, np.int32, np.int8,
                 face_vertex_indices_to_tags)
-    else:
-        facial_adjacency_groups = _compute_facial_adjacency_from_vertices(
-                [grp], boundary_tags, np.int32, np.int8)
 
     from meshmode.mesh import Mesh
     return Mesh(vertices, [grp],
-- 
GitLab


From 12c3cb67500c7799e1e8b9382bd97681cf6ed492 Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Thu, 14 May 2020 10:05:46 -0500
Subject: [PATCH 06/10] Changes naming convention(s), adds check for
 dimensionality

---
 meshmode/mesh/generation.py | 52 ++++++++++++++++++++++++++++---------
 1 file changed, 40 insertions(+), 12 deletions(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 2376cb44..6fc9c8f5 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -658,7 +658,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, boundary_dict={}):
+        group_factory=None, face_to_boundary_tag={}):
     """Create a semi-structured mesh.
 
     :param axis_coords: a tuple with a number of entries corresponding
@@ -666,10 +666,10 @@ 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_dict: an optional dictionary for boundary configuration.
+    :param face_to_boundary_tag: an optional dictionary for boundary configuration.
         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 cardinal
-        directions.
+        a list of the faces on which they should be applied in terms of coordinate
+        directions (+x, -x, +y, -y, +z, -z).
 
     .. versionchanged:: 2017.1
 
@@ -781,28 +781,52 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
     # compute facial adjacency for Mesh if there is tag information
     facial_adjacency_groups = None
     face_vertex_indices_to_tags = {}
-    boundary_tags = list(boundary_dict.keys())
+    boundary_tags = list(face_to_boundary_tag.keys())
     from meshmode.mesh import _compute_facial_adjacency_from_vertices
     if boundary_tags:
         for tag_idx, tag in enumerate(boundary_tags):
             # Need to map the correct face vertices to the boundary tags
-            for face in boundary_dict[tag]:
-                if face == "east":
+            for face in face_to_boundary_tag[tag]:
+                if face == "-x":
                     face_id = 1
                     dim_crit = 0
                     node_crit = axis_coords[0][0]
-                elif face == "west":
+                elif face == "+x":
                     face_id = 1
                     dim_crit = 0
                     node_crit = axis_coords[0][-1]
-                elif face == "north":
+                elif face == "-y":
+                    # Check to make sure we are at least 2D.
+                    if dim < 2:
+                        raise ValueError("Attempting to implement y- "
+                                "boundary condition on 1D box mesh.")
                     face_id = 0
                     dim_crit = 1
                     node_crit = axis_coords[1][-1]
-                elif face == "south":
+                elif face == "+y":
+                    # Check to make sure we are at least 2D.
+                    if dim < 2:
+                        raise ValueError("Attempting to implement y- "
+                                "boundary condition on 1D box mesh.")
                     face_id = 0
                     dim_crit = 1
                     node_crit = axis_coords[1][0]
+                elif face == "-z":
+                    # Check to make sure we are at least 3D.
+                    if dim < 3:
+                        raise ValueError("Attempting to implement z- "
+                                "boundary condition on 2D box mesh.")
+                    face_id = 0
+                    dim_crit = 2
+                    node_crit = axis_coords[2][-1]
+                elif face == "+z":
+                    # Check to make sure we are at least 3D.
+                    if dim < 3:
+                        raise ValueError("Attempting to implement z- "
+                                "boundary condition on 2D box mesh.")
+                    face_id = 0
+                    dim_crit = 2
+                    node_crit = axis_coords[2][0]
                 for ielem in range(0, grp.nelements):
                     if grp.nodes[dim_crit][ielem][0] == node_crit:
                         fvi = grp.vertex_indices[ielem,
@@ -830,13 +854,17 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 # {{{ generate_regular_rect_mesh
 
 def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
-                               boundary_dict={}):
+                               face_to_boundary_tag={}):
     """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 face_to_boundary_tag: an optional dictionary for boundary configuration.
+        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).
     """
     if min(n) < 2:
         raise ValueError("need at least two points in each direction")
@@ -845,7 +873,7 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
             for a_i, b_i, n_i in zip(a, b, n)]
 
     return generate_box_mesh(axis_coords, order=order,
-                             boundary_dict=boundary_dict)
+                             face_to_boundary_tag=face_to_boundary_tag)
 
 # }}}
 
-- 
GitLab


From 96494e51b80c471ef4ad7376af1fed1cd32a7a6f Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Fri, 15 May 2020 13:12:06 -0500
Subject: [PATCH 07/10] Cuts down new code by using axis list

---
 meshmode/mesh/generation.py | 57 +++++++++++--------------------------
 1 file changed, 17 insertions(+), 40 deletions(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 6fc9c8f5..335751e9 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -782,51 +782,28 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
     facial_adjacency_groups = None
     face_vertex_indices_to_tags = {}
     boundary_tags = list(face_to_boundary_tag.keys())
+    axes = ["x", "y", "z"]
+    face_ids = [1, 0, 0]
     from meshmode.mesh import _compute_facial_adjacency_from_vertices
     if boundary_tags:
         for tag_idx, tag in enumerate(boundary_tags):
             # Need to map the correct face vertices to the boundary tags
             for face in face_to_boundary_tag[tag]:
-                if face == "-x":
-                    face_id = 1
-                    dim_crit = 0
-                    node_crit = axis_coords[0][0]
-                elif face == "+x":
-                    face_id = 1
-                    dim_crit = 0
-                    node_crit = axis_coords[0][-1]
-                elif face == "-y":
-                    # Check to make sure we are at least 2D.
-                    if dim < 2:
-                        raise ValueError("Attempting to implement y- "
-                                "boundary condition on 1D box mesh.")
-                    face_id = 0
-                    dim_crit = 1
-                    node_crit = axis_coords[1][-1]
-                elif face == "+y":
-                    # Check to make sure we are at least 2D.
-                    if dim < 2:
-                        raise ValueError("Attempting to implement y- "
-                                "boundary condition on 1D box mesh.")
-                    face_id = 0
-                    dim_crit = 1
-                    node_crit = axis_coords[1][0]
-                elif face == "-z":
-                    # Check to make sure we are at least 3D.
-                    if dim < 3:
-                        raise ValueError("Attempting to implement z- "
-                                "boundary condition on 2D box mesh.")
-                    face_id = 0
-                    dim_crit = 2
-                    node_crit = axis_coords[2][-1]
-                elif face == "+z":
-                    # Check to make sure we are at least 3D.
-                    if dim < 3:
-                        raise ValueError("Attempting to implement z- "
-                                "boundary condition on 2D box mesh.")
-                    face_id = 0
-                    dim_crit = 2
-                    node_crit = axis_coords[2][0]
+                for i_ax, axis in enumerate(axes):
+                    if face == "-" + axis:
+                        if dim < i_ax + 1:
+                            raise ValueError("Boundary condition dimension "
+                            "mismatch")
+                        face_id = face_ids[i_ax]
+                        dim_crit = i_ax
+                        node_crit = axis_coords[i_ax][0]
+                    elif face == "+" + axis:
+                        if dim < i_ax + 1:
+                            raise ValueError("Boundary condition dimension "
+                            "mismatch")
+                        face_id = face_ids[i_ax]
+                        dim_crit = i_ax
+                        node_crit = axis_coords[i_ax][-1]
                 for ielem in range(0, grp.nelements):
                     if grp.nodes[dim_crit][ielem][0] == node_crit:
                         fvi = grp.vertex_indices[ielem,
-- 
GitLab


From d1647c9fac3556dbab53e5223c795fc3245d681a Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Fri, 15 May 2020 13:16:20 -0500
Subject: [PATCH 08/10] Adds example of new face_to_boundary_tag input to
 docstring

---
 meshmode/mesh/generation.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 335751e9..d5153ec6 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -669,7 +669,9 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
     :param face_to_boundary_tag: an optional dictionary for boundary configuration.
         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).
+        directions (+x, -x, +y, -y, +z, -z). One example of this would be:
+
+        face_to_boundary_tag={"bdry_1": ["+x", "+y"], "bdry_2": ["-x"]}
 
     .. versionchanged:: 2017.1
 
-- 
GitLab


From c9835aacd63732a99df6c71c2a615a811e0e267a Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Fri, 15 May 2020 13:33:59 -0500
Subject: [PATCH 09/10] Adds test for 2D box with custom boundary tags

---
 test/test_meshmode.py | 50 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 915b89e8..697e0bbd 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -123,6 +123,56 @@ def test_boundary_tags():
 # }}}
 
 
+# {{{ test custom boundary tags on box mesh
+
+def test_box_boundary_tags():
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    from meshmode.mesh import is_boundary_tag_empty
+    from meshmode.mesh import check_bc_coverage
+    # Test this capability with a 2D mesh.
+    mesh = generate_regular_rect_mesh(a=(0, -1), b=(1, 1),
+                                      n=(20, 20), order=3,
+                                      face_to_boundary_tag={
+                                          "btag_test_1": ["+x", "-y"],
+                                          "btag_test_2": ["+y", "-x"]})
+
+    # ensure boundary is covered
+    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'])
+
+    # correct answers
+    num_on_bdy = 38
+
+    # 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", [
-- 
GitLab


From f0a06f69662e00c5a6107d8ee978bf6644e6bb37 Mon Sep 17 00:00:00 2001
From: Cory Mikida <cmikida2@illinois.edu>
Date: Wed, 27 May 2020 15:33:06 -0500
Subject: [PATCH 10/10] Catches bad face tag specifications, renames some
 things, and implements a 3D test

---
 meshmode/mesh/generation.py | 86 ++++++++++++++++++++++++++-----------
 test/test_meshmode.py       | 36 ++++++++++------
 2 files changed, 85 insertions(+), 37 deletions(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index d5153ec6..5c48f729 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -658,7 +658,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, face_to_boundary_tag={}):
+        group_factory=None, boundary_tag_to_face={}):
     """Create a semi-structured mesh.
 
     :param axis_coords: a tuple with a number of entries corresponding
@@ -666,12 +666,12 @@ 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 face_to_boundary_tag: an optional dictionary for boundary configuration.
+    :param boundary_tag_to_face: an optional dictionary for boundary configuration.
         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). One example of this would be:
+        directions (+x, -x, +y, -y, +z, -z, +w, -w). One example of this would be:
 
-        face_to_boundary_tag={"bdry_1": ["+x", "+y"], "bdry_2": ["-x"]}
+        boundary_tag_to_face={"bdry_1": ["+x", "+y"], "bdry_2": ["-x"]}
 
     .. versionchanged:: 2017.1
 
@@ -713,6 +713,9 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 
     el_vertices = []
 
+    vertex_indices_left = {}
+    vertex_indices_right = {}
+
     if dim == 1:
         for i in range(shape[0]-1):
             # a--b
@@ -723,6 +726,11 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
             el_vertices.append((a, b,))
 
     elif dim == 2:
+        # For custom boundary tagging.
+        vertex_indices_left["x"] = vertex_indices[0, :].flatten()
+        vertex_indices_right["x"] = vertex_indices[-1, :].flatten()
+        vertex_indices_left["y"] = vertex_indices[:, 0].flatten()
+        vertex_indices_right["y"] = vertex_indices[:, -1].flatten()
         for i in range(shape[0]-1):
             for j in range(shape[1]-1):
 
@@ -742,6 +750,13 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
                     el_vertices.append((d, c, b))
 
     elif dim == 3:
+        # For custom boundary tagging.
+        vertex_indices_left["x"] = vertex_indices[0, :, :].flatten()
+        vertex_indices_right["x"] = vertex_indices[-1, :, :].flatten()
+        vertex_indices_left["y"] = vertex_indices[:, 0, :].flatten()
+        vertex_indices_right["y"] = vertex_indices[:, -1, :].flatten()
+        vertex_indices_left["z"] = vertex_indices[:, :, 0].flatten()
+        vertex_indices_right["z"] = vertex_indices[:, :, -1].flatten()
         for i in range(shape[0]-1):
             for j in range(shape[1]-1):
                 for k in range(shape[2]-1):
@@ -783,37 +798,60 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
     # compute facial adjacency for Mesh if there is tag information
     facial_adjacency_groups = None
     face_vertex_indices_to_tags = {}
-    boundary_tags = list(face_to_boundary_tag.keys())
-    axes = ["x", "y", "z"]
-    face_ids = [1, 0, 0]
+    boundary_tags = list(boundary_tag_to_face.keys())
+    axes = ["x", "y", "z", "w"]
     from meshmode.mesh import _compute_facial_adjacency_from_vertices
     if boundary_tags:
         for tag_idx, tag in enumerate(boundary_tags):
             # Need to map the correct face vertices to the boundary tags
-            for face in face_to_boundary_tag[tag]:
+            for face in boundary_tag_to_face[tag]:
+                # Check if face to boundary input is formatted properly
+                if face.startswith("-"):
+                    pass
+                elif face.startswith("+"):
+                    pass
+                else:
+                    raise ValueError("Face ", face, " given for boundary tag ",
+                            tag, " is not formatted properly. Use +x, -x, +y, "
+                            "+z, -z, +w, or -w")
+                if face.endswith("x"):
+                    pass
+                elif face.endswith("y"):
+                    pass
+                elif face.endswith("z"):
+                    pass
+                elif face.endswith("w"):
+                    pass
+                else:
+                    raise ValueError("Face ", face, " given for boundary tag ",
+                            tag, " is not formatted properly. Use +x, -x, +y, "
+                            "+z, -z, +w, or -w")
                 for i_ax, axis in enumerate(axes):
                     if face == "-" + axis:
                         if dim < i_ax + 1:
                             raise ValueError("Boundary condition dimension "
-                            "mismatch")
-                        face_id = face_ids[i_ax]
-                        dim_crit = i_ax
-                        node_crit = axis_coords[i_ax][0]
+                            "mismatch: facial dimension ", face, " for tag ",
+                            tag, " is higher than the dimension of the "
+                            "problem")
+                        vert_crit = vertex_indices_left[axis]
                     elif face == "+" + axis:
                         if dim < i_ax + 1:
                             raise ValueError("Boundary condition dimension "
-                            "mismatch")
-                        face_id = face_ids[i_ax]
-                        dim_crit = i_ax
-                        node_crit = axis_coords[i_ax][-1]
+                            "mismatch: facial dimension ", face, " for tag ",
+                            tag, " is higher than the dimension of the "
+                            "problem")
+                        vert_crit = vertex_indices_right[axis]
                 for ielem in range(0, grp.nelements):
-                    if grp.nodes[dim_crit][ielem][0] == node_crit:
+                    # Loop through potential boundary-facing faces.
+                    for i in range(0, dim+1):
+                        # Element vertices:
                         fvi = grp.vertex_indices[ielem,
                                                  grp.face_vertex_indices()[
-                                                     face_id]]
-                        if frozenset(fvi) not in face_vertex_indices_to_tags:
-                            face_vertex_indices_to_tags[frozenset(fvi)] = []
-                        face_vertex_indices_to_tags[frozenset(fvi)].append(tag)
+                                                     i]]
+                        if all(ind in vert_crit for ind in fvi):
+                            if frozenset(fvi) not in face_vertex_indices_to_tags:
+                                face_vertex_indices_to_tags[frozenset(fvi)] = []
+                            face_vertex_indices_to_tags[frozenset(fvi)].append(tag)
         # Need to add BTAG ALL and BTAG_REALLY_ALL to list of tags
         from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL
         boundary_tags.append(BTAG_ALL)
@@ -833,14 +871,14 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 # {{{ generate_regular_rect_mesh
 
 def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
-                               face_to_boundary_tag={}):
+                               boundary_tag_to_face={}):
     """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 face_to_boundary_tag: an optional dictionary for boundary configuration.
+    :param boundary_tag_to_face: an optional dictionary for boundary configuration.
         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).
@@ -852,7 +890,7 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
             for a_i, b_i, n_i in zip(a, b, n)]
 
     return generate_box_mesh(axis_coords, order=order,
-                             face_to_boundary_tag=face_to_boundary_tag)
+                             boundary_tag_to_face=boundary_tag_to_face)
 
 # }}}
 
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 697e0bbd..7af223d5 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -124,26 +124,36 @@ def test_boundary_tags():
 
 
 # {{{ test custom boundary tags on box mesh
-
-def test_box_boundary_tags():
+@pytest.mark.parametrize(("dim", "nelem"), [
+    (2, 20),
+    (3, 20),
+    ])
+def test_box_boundary_tags(dim, nelem):
     from meshmode.mesh.generation import generate_regular_rect_mesh
     from meshmode.mesh import is_boundary_tag_empty
     from meshmode.mesh import check_bc_coverage
-    # Test this capability with a 2D mesh.
-    mesh = generate_regular_rect_mesh(a=(0, -1), b=(1, 1),
-                                      n=(20, 20), order=3,
-                                      face_to_boundary_tag={
-                                          "btag_test_1": ["+x", "-y"],
-                                          "btag_test_2": ["+y", "-x"]})
+    if 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)
+    # correct answer
+    num_on_bdy = dim*(dim-1)*(nelem-1)**(dim-1)
 
-    # ensure boundary is covered
     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'])
 
-    # correct answers
-    num_on_bdy = 38
-
     # check how many elements are marked on each boundary
     num_marked_bdy_1 = 0
     num_marked_bdy_2 = 0
@@ -165,11 +175,11 @@ def test_box_boundary_tags():
     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))
 
+
 # }}}
 
 
-- 
GitLab