From a5705393cd21ed12c54482c3166ed7e46ee4c6b9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 12 May 2020 22:19:04 -0500 Subject: [PATCH 1/6] 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 From 0ea3200e0bd95358b606c70ff96299af19bd7792 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 13 May 2020 11:27:37 -0500 Subject: [PATCH 2/6] improve is_affine check for general meshes * check all elements, not just the first one * check all second derivatives (including cross ones) * add tests for all of that --- meshmode/mesh/__init__.py | 32 +++++++++++++++++++++++++------- meshmode/mesh/generation.py | 1 + test/test_meshmode.py | 33 +++++++++++++++++++++++++++------ 3 files changed, 53 insertions(+), 13 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 61ad5852..1d2fbfcf 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -225,6 +225,7 @@ class MeshElementGroup(Record): def nunit_nodes(self): return self.unit_nodes.shape[-1] + @property @memoize_method def is_affine(self): return is_affine_group(self) @@ -1360,19 +1361,36 @@ def is_boundary_tag_empty(mesh, boundary_tag): # {{{ -def is_affine_group(group): +def is_affine_group(group, abs_tol=None): + if abs_tol is None: + abs_tol = 1.0e-13 + + # get matrices 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)) + 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]]))) - if not isinstance(diff_matrix, tuple): - diff_matrix = (diff_matrix,) - mat = tuple(vinv.dot(d.dot(d)) for d in diff_matrix) + # check if any element has a non-affine parametrization + 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 - ddx_coeffs = np.einsum("aij,bj->abi", mat, group.nodes[:, 0, :]) - return np.max(np.abs(ddx_coeffs)) < 1.0e-13 + ddx_coeffs = np.einsum("aij,bcj->abci", mats, group.nodes) + norm_inf = np.max(np.abs(ddx_coeffs)) + return norm_inf < abs_tol # }}} diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 92c49b7d..31221b57 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -843,6 +843,7 @@ def generate_warped_rect_mesh(dim, order, n): from meshmode.mesh.processing import map_mesh return map_mesh(mesh, m) + # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 1ae0b634..79c462b9 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1205,8 +1205,25 @@ 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", [ - "box", "warped_box", + "box2d", "box3d", + "warped_box2d", "warped_box3d", "cross_warped_box", "circle", "ellipse", "sphere", "torus" ]) @@ -1219,16 +1236,20 @@ def test_is_affine_group_check(mesh_name): order = 4 nelements = 32 - if mesh_name == "box": - dim = 2 + 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 == "warped_box": + elif mesh_name.startswith("warped_box"): + dim = int(mesh_name[-2]) is_affine = False - dim = 2 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( @@ -1248,7 +1269,7 @@ def test_is_affine_group_check(mesh_name): else: raise ValueError("unknown mesh name: {}".format(mesh_name)) - assert all(grp.is_affine() for grp in mesh.groups) == is_affine + assert all(grp.is_affine for grp in mesh.groups) == is_affine if __name__ == "__main__": -- GitLab From ca94f9fc4e855b7e4a0929c5a981fdd79910fc06 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 17 May 2020 15:49:19 -0500 Subject: [PATCH 3/6] specialize is_affine_check for simplex meshes for now --- meshmode/mesh/__init__.py | 14 +++++++++++--- test/test_meshmode.py | 2 +- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 1d2fbfcf..3bbf2c2c 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -226,9 +226,8 @@ class MeshElementGroup(Record): return self.unit_nodes.shape[-1] @property - @memoize_method def is_affine(self): - return is_affine_group(self) + raise NotImplementedError() def face_vertex_indices(self): """Return a tuple of tuples indicating which vertices @@ -310,6 +309,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 ( @@ -1361,10 +1365,14 @@ def is_boundary_tag_empty(mesh, boundary_tag): # {{{ -def is_affine_group(group, abs_tol=None): +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) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 79c462b9..77a25cd5 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1234,7 +1234,7 @@ def test_is_affine_group_check(mesh_name): generate_icosphere, generate_torus) order = 4 - nelements = 32 + nelements = 16 if mesh_name.startswith("box"): dim = int(mesh_name[-2]) -- GitLab From 0ed195a599c28c0e380ec7dfd7df53c6e6c20e6d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 18 May 2020 02:17:36 +0200 Subject: [PATCH 4/6] Apply suggestion to meshmode/mesh/__init__.py --- meshmode/mesh/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 3bbf2c2c..d0f6df3b 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1396,6 +1396,7 @@ def is_affine_simplex_group(group, abs_tol=None): 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 -- GitLab From 81d43197204d428d8e263d364fa3f2c16220dfd6 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 18 May 2020 02:17:58 +0200 Subject: [PATCH 5/6] Apply suggestion to meshmode/mesh/__init__.py --- meshmode/mesh/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index d0f6df3b..a1009564 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1390,7 +1390,7 @@ def is_affine_simplex_group(group, abs_tol=None): continue mats.append(vinv.dot(diff[n[0]].dot(diff[n[1]]))) - # check if any element has a non-affine parametrization + # 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: -- GitLab From ac0adfdaa4f9762a676b3e5f8ff98f3bb167532f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 17 May 2020 20:58:00 -0500 Subject: [PATCH 6/6] add docs for is_affine flag --- meshmode/discretization/__init__.py | 6 ++++++ meshmode/mesh/__init__.py | 6 +++++- meshmode/mesh/generation.py | 1 - 3 files changed, 11 insertions(+), 2 deletions(-) diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py index ea282587..a3847f2c 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): diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index a1009564..0af14583 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -159,7 +159,11 @@ class MeshElementGroup(Record): *Not* the ambient dimension, see :attr:`Mesh.ambient_dim` for that. - .. automethod:: is_affine + .. 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 diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 31221b57..92c49b7d 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -843,7 +843,6 @@ def generate_warped_rect_mesh(dim, order, n): from meshmode.mesh.processing import map_mesh return map_mesh(mesh, m) - # }}} -- GitLab