diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py
index c4881c4f366619b4484156ae740cca3728ddc590..e7929937d1fdf7d360759e477aee7a10faa19ce3 100644
--- a/meshmode/discretization/connection/__init__.py
+++ b/meshmode/discretization/connection/__init__.py
@@ -31,8 +31,9 @@ from pytools import memoize_method, memoize_in
from meshmode.discretization.connection.same_mesh import \
make_same_mesh_connection
-from meshmode.discretization.connection.face import \
- make_face_restriction
+from meshmode.discretization.connection.face import (
+ FRESTR_INTERIOR_FACES, FRESTR_ALL_FACES,
+ make_face_restriction)
from meshmode.discretization.connection.opposite_face import \
make_opposite_face_connection
@@ -42,7 +43,9 @@ logger = logging.getLogger(__name__)
__all__ = [
"DiscretizationConnection",
- "make_same_mesh_connection", "make_face_restriction",
+ "make_same_mesh_connection",
+ "FRESTR_INTERIOR_FACES", "FRESTR_ALL_FACES",
+ "make_face_restriction",
"make_opposite_face_connection"
]
@@ -51,6 +54,8 @@ __doc__ = """
.. autofunction:: make_same_mesh_connection
+.. autofunction:: FRESTR_INTERIOR_FACES
+.. autofunction:: FRESTR_ALL_FACES
.. autofunction:: make_face_restriction
.. autofunction:: make_opposite_face_connection
diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py
index f03b5d8a744770d11192adbcf13b70e15ec56e62..14127b32a28396236e5d457bd079ac4d2f79ed43 100644
--- a/meshmode/discretization/connection/face.py
+++ b/meshmode/discretization/connection/face.py
@@ -36,49 +36,74 @@ import logging
logger = logging.getLogger(__name__)
+class FRESTR_INTERIOR_FACES:
+ """A special value to pass to
+ :func:`meshmode.discretization.connection.make_face_restriction`
+ to produce a discretization consisting of all interior faces
+ in a discretization.
+ """
+
+
+class FRESTR_ALL_FACES:
+ """A special value to pass to
+ :func:`meshmode.discretization.connection.make_face_restriction`
+ to produce a discretization consisting of all faces (interior and boundary)
+ faces in a discretization.
+ """
+
+
# {{{ boundary connection
class _ConnectionBatchData(Record):
pass
-def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data):
+def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data,
+ per_face_groups):
from meshmode.discretization.connection import (
InterpolationBatch, DiscretizationConnectionElementGroup,
DiscretizationConnection)
+ ibdry_grp = 0
+ batches = []
+
connection_groups = []
- for igrp, (vol_grp, bdry_grp) in enumerate(
- zip(vol_discr.groups, bdry_discr.groups)):
- connection_batches = []
+ for igrp, vol_grp in enumerate(vol_discr.groups):
mgrp = vol_grp.mesh_el_group
- for face_id in range(len(mgrp.face_vertex_indices())):
+ for face_id in range(mgrp.nfaces):
+ bdry_grp = bdry_discr.groups[ibdry_grp]
data = connection_data[igrp, face_id]
bdry_unit_nodes_01 = (bdry_grp.unit_nodes + 1)*0.5
result_unit_nodes = (np.dot(data.A, bdry_unit_nodes_01).T + data.b).T
- connection_batches.append(
- InterpolationBatch(
- from_group_index=igrp,
- from_element_indices=cl.array.to_device(
- queue,
- vol_grp.mesh_el_group.element_nr_base
- + data.group_source_element_indices)
- .with_queue(None),
- to_element_indices=cl.array.to_device(
- queue,
- bdry_grp.mesh_el_group.element_nr_base
- + data.group_target_element_indices)
- .with_queue(None),
- result_unit_nodes=result_unit_nodes,
- to_element_face=face_id
- ))
-
- connection_groups.append(
- DiscretizationConnectionElementGroup(
- connection_batches))
+ batches.append(
+ InterpolationBatch(
+ from_group_index=igrp,
+ from_element_indices=cl.array.to_device(
+ queue,
+ vol_grp.mesh_el_group.element_nr_base
+ + data.group_source_element_indices)
+ .with_queue(None),
+ to_element_indices=cl.array.to_device(
+ queue,
+ data.group_target_element_indices)
+ .with_queue(None),
+ result_unit_nodes=result_unit_nodes,
+ to_element_face=face_id
+ ))
+
+ is_last_face = face_id + 1 == mgrp.nfaces
+
+ if per_face_groups or is_last_face:
+ connection_groups.append(
+ DiscretizationConnectionElementGroup(batches))
+ batches = []
+
+ ibdry_grp += 1
+
+ assert ibdry_grp == len(bdry_discr.groups)
return DiscretizationConnection(
vol_discr, bdry_discr, connection_groups,
@@ -93,7 +118,7 @@ def _get_face_vertices(mesh, boundary_tag):
# a set of volume vertex numbers
bdry_vertex_vol_nrs = set()
- if boundary_tag is not None:
+ if boundary_tag not in [FRESTR_INTERIOR_FACES, FRESTR_ALL_FACES]:
# {{{ boundary faces
btag_bit = mesh.boundary_tag_bit(boundary_tag)
@@ -121,15 +146,17 @@ def _get_face_vertices(mesh, boundary_tag):
# }}}
else:
- # For interior faces, this is likely every vertex in the book.
+ # For INTERIOR_FACES, this is likely every vertex in the book.
# Don't ever bother trying to cut the list down.
+ # For ALL_FACES, it literally is every single vertex.
return np.arange(mesh.nvertices, dtype=np.intp)
# }}}
-def make_face_restriction(discr, group_factory, boundary_tag):
+def make_face_restriction(discr, group_factory, boundary_tag,
+ per_face_groups=False):
"""Create a mesh, a discretization and a connection to restrict
a function on *discr* to its values on the edges of element faces
denoted by *boundary_tag*.
@@ -142,6 +169,17 @@ def make_face_restriction(discr, group_factory, boundary_tag):
to make a discretization consisting of all (interior and
boundary) faces.
+ :arg per_face_groups: If *True*, the resulting discretization is
+ guaranteed to have groups organized as::
+
+ (grp0, face0), (grp0, face1), ... (grp0, faceN),
+ (grp1, face0), (grp1, face1), ... (grp1, faceN), ...
+
+ each with the elements in the same order as the originating
+ group. If *False*, volume and boundary groups correspond with
+ each other one-to-one, and an interpolation batch is created
+ per face.
+
:return: a :class:`meshmode.discretization.connection.DiscretizationConnection`
representing the new connection. The new boundary discretization can be
obtained from the
@@ -149,16 +187,15 @@ def make_face_restriction(discr, group_factory, boundary_tag):
attribute of the return value, and the corresponding new boundary mesh
from that.
- The resulting discretization is guaranteed to have groups
- organized as::
-
- (grp0, face0), (grp0, face1), ... (grp0, faceN),
- (grp1, face0), (grp1, face1), ... (grp1, faceN), ...
-
- each with the elements in the same order as the originating
- group.
"""
+ if boundary_tag is None:
+ boundary_tag = FRESTR_INTERIOR_FACES
+ from warnings import warn
+ warn("passing *None* for boundary_tag is deprecated--pass "
+ "INTERIOR_FACES instead",
+ DeprecationWarning, stacklevel=2)
+
logger.info("building face restriction: start")
# {{{ gather boundary vertices
@@ -195,18 +232,7 @@ def make_face_restriction(discr, group_factory, boundary_tag):
group_boundary_faces = []
- if boundary_tag is not None:
- bdry_grp = fagrp_map.get(None)
- if bdry_grp is not None:
- nb_el_bits = -bdry_grp.neighbors
- face_relevant_flags = (nb_el_bits & btag_bit) != 0
-
- group_boundary_faces.extend(
- zip(
- bdry_grp.elements[face_relevant_flags],
- bdry_grp.element_faces[face_relevant_flags]))
-
- else:
+ if boundary_tag is FRESTR_INTERIOR_FACES:
for fagrp in six.itervalues(fagrp_map):
if fagrp.ineighbor_group is None:
# boundary faces -> not looking for those
@@ -215,34 +241,34 @@ def make_face_restriction(discr, group_factory, boundary_tag):
group_boundary_faces.extend(
zip(fagrp.elements, fagrp.element_faces))
- # }}}
-
- # {{{ Preallocate arrays for mesh group
-
- ngroup_bdry_elements = len(group_boundary_faces)
- vertex_indices = np.empty(
- (ngroup_bdry_elements, mgrp.dim+1-1),
- mgrp.vertex_indices.dtype)
+ elif boundary_tag is FRESTR_ALL_FACES:
+ group_boundary_faces.extend(
+ (iel, iface)
+ for iface in range(grp.mesh_el_group.nfaces)
+ for iel in range(grp.nelements)
+ )
- bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order)
- bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5
+ else:
+ bdry_grp = fagrp_map.get(None)
+ if bdry_grp is not None:
+ nb_el_bits = -bdry_grp.neighbors
+ face_relevant_flags = (nb_el_bits & btag_bit) != 0
- vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order)
- nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1]
- nodes = np.empty(
- (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes),
- dtype=np.float64)
+ group_boundary_faces.extend(
+ zip(
+ bdry_grp.elements[face_relevant_flags],
+ bdry_grp.element_faces[face_relevant_flags]))
# }}}
grp_face_vertex_indices = mgrp.face_vertex_indices()
grp_vertex_unit_coordinates = mgrp.vertex_unit_coordinates()
- # batch by face_id
-
batch_base = 0
- for face_id in range(len(grp_face_vertex_indices)):
+ # group by face_id
+
+ for face_id in range(mgrp.nfaces):
batch_boundary_el_numbers_in_grp = np.array(
[
ibface_el
@@ -250,9 +276,34 @@ def make_face_restriction(discr, group_factory, boundary_tag):
if ibface_face == face_id],
dtype=np.intp)
- new_el_numbers = np.arange(
- batch_base,
- batch_base + len(batch_boundary_el_numbers_in_grp))
+ # {{{ Preallocate arrays for mesh group
+
+ nbatch_elements = len(batch_boundary_el_numbers_in_grp)
+
+ if per_face_groups or face_id == 0:
+ if per_face_groups:
+ ngroup_bdry_elements = nbatch_elements
+ else:
+ ngroup_bdry_elements = len(group_boundary_faces)
+
+ vertex_indices = np.empty(
+ (ngroup_bdry_elements, mgrp.dim+1-1),
+ mgrp.vertex_indices.dtype)
+
+ bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order)
+ bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5
+
+ vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order)
+ nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1]
+ nodes = np.empty(
+ (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes),
+ dtype=np.float64)
+
+ # }}}
+
+ new_el_numbers = batch_base + np.arange(nbatch_elements)
+ if not per_face_groups:
+ batch_base += nbatch_elements
# {{{ no per-element axes in these computations
@@ -303,11 +354,13 @@ def make_face_restriction(discr, group_factory, boundary_tag):
b=b,
)
- batch_base += len(batch_boundary_el_numbers_in_grp)
+ is_last_face = face_id + 1 == mgrp.nfaces
- bdry_mesh_group = SimplexElementGroup(
- mgrp.order, vertex_indices, nodes, unit_nodes=bdry_unit_nodes)
- bdry_mesh_groups.append(bdry_mesh_group)
+ if per_face_groups or is_last_face:
+ bdry_mesh_group = SimplexElementGroup(
+ mgrp.order, vertex_indices, nodes,
+ unit_nodes=bdry_unit_nodes)
+ bdry_mesh_groups.append(bdry_mesh_group)
bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups)
@@ -317,7 +370,8 @@ def make_face_restriction(discr, group_factory, boundary_tag):
with cl.CommandQueue(discr.cl_context) as queue:
connection = _build_boundary_connection(
- queue, discr, bdry_discr, connection_data)
+ queue, discr, bdry_discr, connection_data,
+ per_face_groups)
logger.info("building face restriction: done")
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 03a6661d3d032d4aeefbd2b5e6eaea7b055bcd24..122b0b9666f958afa14ac0fc1b6e6c8cd159e636 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -38,6 +38,8 @@ from meshmode.discretization.poly_element import (
PolynomialWarpAndBlendGroupFactory
)
from meshmode.mesh import BTAG_ALL
+from meshmode.discretization.connection import \
+ FRESTR_ALL_FACES, FRESTR_INTERIOR_FACES
import pytest
@@ -80,16 +82,17 @@ def test_circle_mesh(do_plot=False):
])
@pytest.mark.parametrize("boundary_tag", [
BTAG_ALL,
- # all internal faces:
- None
+ FRESTR_ALL_FACES,
+ FRESTR_INTERIOR_FACES,
])
@pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
("blob", 2, [1e-1, 8e-2, 5e-2]),
("warp", 2, [10, 20, 30]),
("warp", 3, [10, 20, 30]),
])
+@pytest.mark.parametrize("per_face_groups", [False, True])
def test_boundary_interpolation(ctx_getter, group_factory, boundary_tag,
- mesh_name, dim, mesh_pars):
+ mesh_name, dim, mesh_pars, per_face_groups):
cl_ctx = ctx_getter()
queue = cl.CommandQueue(cl_ctx)
@@ -141,7 +144,7 @@ def test_boundary_interpolation(ctx_getter, group_factory, boundary_tag,
bdry_connection = make_face_restriction(
vol_discr, group_factory(order),
- boundary_tag)
+ boundary_tag, per_face_groups=per_face_groups)
bdry_discr = bdry_connection.to_discr
bdry_x = bdry_discr.nodes()[0].with_queue(queue)
@@ -152,7 +155,8 @@ def test_boundary_interpolation(ctx_getter, group_factory, boundary_tag,
mat = bdry_connection.full_resample_matrix(queue).get(queue)
bdry_f_2_by_mat = mat.dot(vol_f.get())
- assert(la.norm(bdry_f_2.get(queue=queue) - bdry_f_2_by_mat)) < 1e-14
+ mat_error = la.norm(bdry_f_2.get(queue=queue) - bdry_f_2_by_mat)
+ assert mat_error < 1e-14, mat_error
err = la.norm((bdry_f-bdry_f_2).get(), np.inf)
eoc_rec.add_data_point(h, err)