diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 679d98b34184cc2b596b39640527491796f75c60..170b0a006fb8ea8ac612a7759f9e2fb0abe195f5 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -33,110 +33,161 @@ __doc__ = """
"""
-def find_volume_mesh_element_orientations(mesh):
- result = np.empty(mesh.nelements, dtype=np.float64)
+def find_volume_mesh_element_group_orientation(mesh, grp):
+ """Return a positive floating point number for each positively
+ oriented element, and a negative floating point number for
+ each negatively oriented element.
+ """
from meshmode.mesh import SimplexElementGroup
- for grp in mesh.groups:
- if not isinstance(grp, SimplexElementGroup):
- raise NotImplementedError(
- "finding element orientations "
- "only supported on "
- "exclusively SimplexElementGroup-based meshes")
+ if not isinstance(grp, SimplexElementGroup):
+ raise NotImplementedError(
+ "finding element orientations "
+ "only supported on "
+ "exclusively SimplexElementGroup-based meshes")
+
+ # (ambient_dim, nelements, nvertices)
+ vertices = mesh.vertices[:, grp.vertex_indices]
+
+ # (ambient_dim, nelements, nspan_vectors)
+ spanning_vectors = (
+ vertices[:, :, 1:] - vertices[:, :, 0][:, :, np.newaxis])
- # (ambient_dim, nelements, nvertices)
- vertices = mesh.vertices[:, grp.vertex_indices]
+ ambient_dim = spanning_vectors.shape[0]
+ nspan_vectors = spanning_vectors.shape[-1]
- # (ambient_dim, nelements, nspan_vectors)
- spanning_vectors = (
- vertices[:, :, 1:] - vertices[:, :, 0][:, :, np.newaxis])
+ spanning_object_array = np.empty(
+ (nspan_vectors, ambient_dim),
+ dtype=np.object)
- ambient_dim = spanning_vectors.shape[0]
- nspan_vectors = spanning_vectors.shape[-1]
+ for ispan in xrange(nspan_vectors):
+ for idim in xrange(ambient_dim):
+ spanning_object_array[ispan, idim] = \
+ spanning_vectors[idim, :, ispan]
- spanning_object_array = np.empty(
- (nspan_vectors, ambient_dim),
- dtype=np.object)
+ from pymbolic.geometric_algebra import MultiVector
- for ispan in xrange(nspan_vectors):
- for idim in xrange(ambient_dim):
- spanning_object_array[ispan, idim] = \
- spanning_vectors[idim, :, ispan]
+ mvs = [MultiVector(vec) for vec in spanning_object_array]
- from pymbolic.geometric_algebra import MultiVector
+ from operator import xor
+ outer_prod = -reduce(xor, mvs)
- mvs = [MultiVector(vec) for vec in spanning_object_array]
+ return (outer_prod.I | outer_prod).as_scalar()
- from operator import xor
- outer_prod = reduce(xor, mvs)
- signed_area_element = (outer_prod.I | outer_prod).as_scalar()
+def find_volume_mesh_element_orientations(mesh, tolerate_unimplemented_checks=False):
+ """Return a positive floating point number for each positively
+ oriented element, and a negative floating point number for
+ each negatively oriented element.
- result[grp.element_nr_base:grp.element_nr_base + grp.nelements] \
- = signed_area_element
+ :arg tolerate_unimplemented_checks: If *True*, elements for which no
+ check is available will return NaN.
+ """
+
+ result = np.empty(mesh.nelements, dtype=np.float64)
+
+ for grp in mesh.groups:
+ result_grp_view = result[
+ grp.element_nr_base:grp.element_nr_base + grp.nelements]
+
+ if tolerate_unimplemented_checks:
+ try:
+ signed_area_elements = \
+ find_volume_mesh_element_group_orientation(mesh, grp)
+ except NotImplementedError:
+ result_grp_view[:] = float("nan")
+ else:
+ assert not np.isnan(signed_area_elements).any()
+ result_grp_view[:] = signed_area_elements
+ else:
+ signed_area_elements = \
+ find_volume_mesh_element_group_orientation(mesh, grp)
+ assert not np.isnan(signed_area_elements).any()
+ result_grp_view[:] = signed_area_elements
return result
+def test_volume_mesh_element_orientations(mesh):
+ area_elements = find_volume_mesh_element_orientations(
+ mesh, tolerate_unimplemented_checks=True)
+
+ valid = ~np.isnan(area_elements)
+
+ return (area_elements[valid] > 0).all()
+
+
# {{{ flips
-def perform_flips(mesh, flip_flags):
+def flip_simplex_element_group(vertices, grp, grp_flip_flags):
from modepy.tools import barycentric_to_unit, unit_to_barycentric
- flip_flags = flip_flags.astype(np.bool)
+ from meshmode.mesh import SimplexElementGroup
- from meshmode.mesh import Mesh, SimplexElementGroup
+ if not isinstance(grp, SimplexElementGroup):
+ raise NotImplementedError("flips only supported on "
+ "exclusively SimplexElementGroup-based meshes")
- new_groups = []
- for grp in mesh.groups:
- if not isinstance(grp, SimplexElementGroup):
- raise NotImplementedError("flips only supported on "
- "exclusively SimplexElementGroup-based meshes")
+ # Swap the first two vertices on elements to be flipped.
- grp_flip_flags = flip_flags[
- grp.element_nr_base:grp.element_nr_base+grp.nelements]
+ new_vertex_indices = grp.vertex_indices.copy()
+ new_vertex_indices[grp_flip_flags, 0] \
+ = grp.vertex_indices[grp_flip_flags, 1]
+ new_vertex_indices[grp_flip_flags, 1] \
+ = grp.vertex_indices[grp_flip_flags, 0]
- # Swap the first two vertices on elements to be flipped.
+ # Generate a resampling matrix that corresponds to the
+ # first two barycentric coordinates being swapped.
- new_vertex_indices = grp.vertex_indices.copy()
- new_vertex_indices[grp_flip_flags, 0] \
- = grp.vertex_indices[grp_flip_flags, 1]
- new_vertex_indices[grp_flip_flags, 1] \
- = grp.vertex_indices[grp_flip_flags, 0]
+ bary_unit_nodes = unit_to_barycentric(grp.unit_nodes)
- # Generate a resampling matrix that corresponds to the
- # first two barycentric coordinates being swapped.
+ flipped_bary_unit_nodes = bary_unit_nodes.copy()
+ flipped_bary_unit_nodes[0, :] = bary_unit_nodes[1, :]
+ flipped_bary_unit_nodes[1, :] = bary_unit_nodes[0, :]
+ flipped_unit_nodes = barycentric_to_unit(flipped_bary_unit_nodes)
- bary_unit_nodes = unit_to_barycentric(grp.unit_nodes)
+ flip_matrix = mp.resampling_matrix(
+ mp.simplex_onb(grp.dim, grp.order),
+ flipped_unit_nodes, grp.unit_nodes)
- flipped_bary_unit_nodes = bary_unit_nodes.copy()
- flipped_bary_unit_nodes[0, :] = bary_unit_nodes[1, :]
- flipped_bary_unit_nodes[1, :] = bary_unit_nodes[0, :]
- flipped_unit_nodes = barycentric_to_unit(flipped_bary_unit_nodes)
+ flip_matrix[np.abs(flip_matrix) < 1e-15] = 0
- flip_matrix = mp.resampling_matrix(
- mp.simplex_onb(grp.dim, grp.order),
- flipped_unit_nodes, grp.unit_nodes)
+ # Flipping twice should be the identity
+ assert la.norm(
+ np.dot(flip_matrix, flip_matrix)
+ - np.eye(len(flip_matrix))) < 1e-13
- flip_matrix[np.abs(flip_matrix) < 1e-15] = 0
+ # Apply the flip matrix to the nodes.
+ new_nodes = grp.nodes.copy()
+ new_nodes[:, grp_flip_flags] = np.einsum(
+ "ij,dej->dei",
+ flip_matrix, grp.nodes[:, grp_flip_flags])
- # Flipping twice should be the identity
- assert la.norm(
- np.dot(flip_matrix, flip_matrix)
- - np.eye(len(flip_matrix))) < 1e-13
+ return SimplexElementGroup(
+ grp.order, new_vertex_indices, new_nodes,
+ unit_nodes=grp.unit_nodes)
- # Apply the flip matrix to the nodes.
- new_nodes = grp.nodes.copy()
- new_nodes[:, grp_flip_flags] = np.einsum(
- "ij,dej->dei",
- flip_matrix, grp.nodes[:, grp_flip_flags])
- new_groups.append(SimplexElementGroup(
- grp.order, new_vertex_indices, new_nodes,
- unit_nodes=grp.unit_nodes))
+def perform_flips(mesh, flip_flags, skip_tests=False):
+ flip_flags = flip_flags.astype(np.bool)
+
+ from meshmode.mesh import Mesh
+
+ new_groups = []
+ for grp in mesh.groups:
+ grp_flip_flags = flip_flags[
+ grp.element_nr_base:grp.element_nr_base+grp.nelements]
+
+ if grp_flip_flags.any():
+ new_grp = flip_simplex_element_group(
+ mesh.vertices, grp, grp_flip_flags)
+ else:
+ new_grp = grp.copy()
+
+ new_groups.append(new_grp)
- return Mesh(mesh.vertices, new_groups)
+ return Mesh(mesh.vertices, new_groups, skip_tests=skip_tests)
# }}}