diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index 6a4fcda7565487d44ee49045871d9226762af2db..a3847f2c3ab6c262c390f3c805619d74c85ab401 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -67,6 +67,12 @@ class ElementGroupBase(object):
Returns an array of length :attr:`nunit_nodes` containing
quadrature weights.
+
+ .. attribute:: is_affine
+
+ A :class:`bool` flag that is *True* if the local-to-global
+ parametrization of all the elements in the group is affine. Based on
+ :attr:`meshmode.mesh.MeshElementGroup.is_affine`.
"""
def __init__(self, mesh_el_group, order, node_nr_base):
@@ -78,6 +84,10 @@ class ElementGroupBase(object):
self.order = order
self.node_nr_base = node_nr_base
+ @property
+ def is_affine(self):
+ return self.mesh_el_group.is_affine
+
@property
def nelements(self):
return self.mesh_el_group.nelements
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 19ef8e5ddaf46b283ea9ca3431d672f57deca4fb..9df70d387b5d40d435982024c1dd0b78429915d8 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -26,9 +26,10 @@ from six.moves import range
import six
import numpy as np
-import modepy as mp
import numpy.linalg as la
-from pytools import Record
+
+import modepy as mp
+from pytools import Record, memoize_method
__doc__ = """
@@ -162,6 +163,11 @@ class MeshElementGroup(Record):
*Not* the ambient dimension, see :attr:`Mesh.ambient_dim`
for that.
+ .. attribute:: is_affine
+
+ A :class:`bool` flag that is *True* if the local-to-global
+ parametrization of all the elements in the group is affine.
+
.. automethod:: face_vertex_indices
.. automethod:: vertex_unit_coordinates
@@ -227,6 +233,10 @@ class MeshElementGroup(Record):
def nunit_nodes(self):
return self.unit_nodes.shape[-1]
+ @property
+ def is_affine(self):
+ raise NotImplementedError()
+
def face_vertex_indices(self):
"""Return a tuple of tuples indicating which vertices
(in mathematically positive ordering) make up each face
@@ -307,6 +317,11 @@ class SimplexElementGroup(MeshElementGroup):
super(SimplexElementGroup, self).__init__(order, vertex_indices, nodes,
element_nr_base, node_nr_base, unit_nodes, dim)
+ @property
+ @memoize_method
+ def is_affine(self):
+ return is_affine_simplex_group(self)
+
def face_vertex_indices(self):
if self.dim == 1:
return (
@@ -1367,4 +1382,44 @@ def is_boundary_tag_empty(mesh, boundary_tag):
# }}}
+# {{{
+
+def is_affine_simplex_group(group, abs_tol=None):
+ if abs_tol is None:
+ abs_tol = 1.0e-13
+
+ if not isinstance(group, SimplexElementGroup):
+ raise TypeError("expected a 'SimplexElementGroup' not '%s'" %
+ type(group).__name__)
+
+ # get matrices
+ basis = mp.simplex_best_available_basis(group.dim, group.order)
+ grad_basis = mp.grad_simplex_best_available_basis(group.dim, group.order)
+
+ vinv = la.inv(mp.vandermonde(basis, group.unit_nodes))
+ diff = mp.differentiation_matrices(basis, grad_basis, group.unit_nodes)
+ if not isinstance(diff, tuple):
+ diff = (diff,)
+
+ # construct all second derivative matrices (including cross terms)
+ from itertools import product
+ mats = []
+ for n in product(range(group.dim), repeat=2):
+ if n[0] > n[1]:
+ continue
+ mats.append(vinv.dot(diff[n[0]].dot(diff[n[1]])))
+
+ # check just the first element for a non-affine local-to-global mapping
+ ddx_coeffs = np.einsum("aij,bj->abi", mats, group.nodes[:, 0, :])
+ norm_inf = np.max(np.abs(ddx_coeffs))
+ if norm_inf > abs_tol:
+ return False
+
+ # check all elements for a non-affine local-to-global mapping
+ ddx_coeffs = np.einsum("aij,bcj->abci", mats, group.nodes)
+ norm_inf = np.max(np.abs(ddx_coeffs))
+ return norm_inf < abs_tol
+
+# }}}
+
# vim: foldmethod=marker
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index e2ee1a813b27629aba4bf31ec250dbcf758f7b1c..77a25cd5dd4de07e1207ca137ec5b1510d268a84 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -1205,6 +1205,73 @@ def test_open_curved_mesh(curve_name):
closed=closed)
+def _generate_cross_warped_rect_mesh(dim, order, n):
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(0,)*dim, b=(1,)*dim,
+ n=(n,)*dim, order=order)
+
+ def m(x):
+ results = np.empty_like(x)
+ results[0] = 1 + 1.5 * (x[0] + 0.25) * (x[1] + 0.3)
+ results[1] = x[1]
+ return results
+
+ from meshmode.mesh.processing import map_mesh
+ return map_mesh(mesh, m)
+
+
+@pytest.mark.parametrize("mesh_name", [
+ "box2d", "box3d",
+ "warped_box2d", "warped_box3d", "cross_warped_box",
+ "circle", "ellipse",
+ "sphere", "torus"
+ ])
+def test_is_affine_group_check(mesh_name):
+ from meshmode.mesh.generation import (
+ generate_regular_rect_mesh, generate_warped_rect_mesh,
+ make_curve_mesh, ellipse,
+ generate_icosphere, generate_torus)
+
+ order = 4
+ nelements = 16
+
+ if mesh_name.startswith("box"):
+ dim = int(mesh_name[-2])
+ is_affine = True
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*dim, b=(0.5,)*dim,
+ n=(nelements,)*dim, order=order)
+ elif mesh_name.startswith("warped_box"):
+ dim = int(mesh_name[-2])
+ is_affine = False
+ mesh = generate_warped_rect_mesh(dim, order, nelements)
+ elif mesh_name == "cross_warped_box":
+ dim = 2
+ is_affine = False
+ mesh = _generate_cross_warped_rect_mesh(dim, order, nelements)
+ elif mesh_name == "circle":
+ is_affine = False
+ mesh = make_curve_mesh(
+ lambda t: ellipse(1.0, t),
+ np.linspace(0.0, 1.0, nelements + 1), order=order)
+ elif mesh_name == "ellipse":
+ is_affine = False
+ mesh = make_curve_mesh(
+ lambda t: ellipse(2.0, t),
+ np.linspace(0.0, 1.0, nelements + 1), order=order)
+ elif mesh_name == "sphere":
+ is_affine = False
+ mesh = generate_icosphere(r=1.0, order=order)
+ elif mesh_name == "torus":
+ is_affine = False
+ mesh = generate_torus(10.0, 2.0, order=order)
+ else:
+ raise ValueError("unknown mesh name: {}".format(mesh_name))
+
+ assert all(grp.is_affine for grp in mesh.groups) == is_affine
+
+
if __name__ == "__main__":
import sys
if len(sys.argv) > 1: