From a5705393cd21ed12c54482c3166ed7e46ee4c6b9 Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Tue, 12 May 2020 22:19:04 -0500
Subject: [PATCH] add check if mesh group is affine
---
meshmode/discretization/__init__.py | 4 +++
meshmode/mesh/__init__.py | 28 ++++++++++++++++--
test/test_meshmode.py | 46 +++++++++++++++++++++++++++++
3 files changed, 76 insertions(+), 2 deletions(-)
diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index 6a4fcda7..ea282587 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -78,6 +78,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 211e764c..61ad5852 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,7 @@ class MeshElementGroup(Record):
*Not* the ambient dimension, see :attr:`Mesh.ambient_dim`
for that.
+ .. automethod:: is_affine
.. automethod:: face_vertex_indices
.. automethod:: vertex_unit_coordinates
@@ -223,6 +225,10 @@ class MeshElementGroup(Record):
def nunit_nodes(self):
return self.unit_nodes.shape[-1]
+ @memoize_method
+ def is_affine(self):
+ return is_affine_group(self)
+
def face_vertex_indices(self):
"""Return a tuple of tuples indicating which vertices
(in mathematically positive ordering) make up each face
@@ -1352,4 +1358,22 @@ def is_boundary_tag_empty(mesh, boundary_tag):
# }}}
+# {{{
+
+def is_affine_group(group):
+ basis = mp.simplex_best_available_basis(group.dim, group.order)
+ grad_basis = mp.grad_simplex_best_available_basis(group.dim, group.order)
+
+ diff_matrix = mp.differentiation_matrices(basis, grad_basis, group.unit_nodes)
+ vinv = la.inv(mp.vandermonde(basis, group.unit_nodes))
+
+ if not isinstance(diff_matrix, tuple):
+ diff_matrix = (diff_matrix,)
+ mat = tuple(vinv.dot(d.dot(d)) for d in diff_matrix)
+
+ ddx_coeffs = np.einsum("aij,bj->abi", mat, group.nodes[:, 0, :])
+ return np.max(np.abs(ddx_coeffs)) < 1.0e-13
+
+# }}}
+
# vim: foldmethod=marker
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index e2ee1a81..1ae0b634 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -1205,6 +1205,52 @@ def test_open_curved_mesh(curve_name):
closed=closed)
+@pytest.mark.parametrize("mesh_name", [
+ "box", "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 = 32
+
+ if mesh_name == "box":
+ dim = 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 == "warped_box":
+ is_affine = False
+ dim = 2
+ mesh = generate_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:
--
GitLab