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 211e764c87b95310d5b68995677f0d958a8c9550..0af14583cd9a14bfae0ca0060d4c72f331508d90 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__ = """ @@ -158,6 +159,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 @@ -223,6 +229,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 @@ -303,6 +313,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 ( @@ -1352,4 +1367,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: