From 5ee016f285b86ee06121875869132be0c5dc7eb3 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 11 May 2016 15:40:18 -0500 Subject: [PATCH] Support, test 4D simplicial meshes --- .../discretization/connection/__init__.py | 2 +- meshmode/discretization/poly_element.py | 52 ++++++++++++++---- meshmode/mesh/__init__.py | 32 ++++------- meshmode/mesh/generation.py | 6 ++- test/test_meshmode.py | 53 ++++++++++++++++++- 5 files changed, 109 insertions(+), 36 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index c4bcf4a..269ed5d 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -212,7 +212,7 @@ class DiscretizationConnection(object): from_grp = self.from_discr.groups[ibatch.from_group_index] result = mp.resampling_matrix( - mp.simplex_onb(self.from_discr.dim, from_grp.order), + from_grp.basis(), ibatch.result_unit_nodes, from_grp.unit_nodes) with cl.CommandQueue(self.cl_context) as queue: diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index ac8629d..a5319bc 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -38,6 +38,7 @@ Group types .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup .. autoclass:: PolynomialWarpAndBlendElementGroup +.. autoclass:: PolynomialEquidistantElementGroup Group factories ^^^^^^^^^^^^^^^ @@ -45,6 +46,7 @@ Group factories .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory .. autoclass:: PolynomialWarpAndBlendGroupFactory +.. autoclass:: PolynomialEquidistantGroupFactory """ # FIXME Most of the loopy kernels will break as soon as we start using multiple @@ -63,16 +65,16 @@ from meshmode.discretization import ElementGroupBase class PolynomialSimplexElementGroupBase(ElementGroupBase): def basis(self): - return mp.simplex_onb(self.dim, self.order) + return mp.simplex_best_available_basis(self.dim, self.order) def grad_basis(self): - return mp.grad_simplex_onb(self.dim, self.order) + return mp.grad_simplex_best_available_basis(self.dim, self.order) @memoize_method def from_mesh_interp_matrix(self): meg = self.mesh_el_group return mp.resampling_matrix( - mp.simplex_onb(meg.dim, meg.order), + mp.simplex_best_available_basis(meg.dim, meg.order), self.unit_nodes, meg.unit_nodes) @@ -150,11 +152,20 @@ class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): return self._quadrature_rule().weights -class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase): +class _MassMatrixQuadratureElementGroup(PolynomialSimplexElementGroupBase): + @property + @memoize_method + def weights(self): + return np.dot( + mp.mass_matrix(self.basis(), self.unit_nodes), + np.ones(len(self.basis()))) + + +class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): """Elemental discretization with a number of nodes matching the number of polynomials in :math:`P^k`, hence usable for differentiation and - interpolation. Interpolation nodes are present on the boundary of the - simplex. + interpolation. Interpolation nodes edge-clustered for avoidance of Runge + phenomena. Nodes are present on the boundary of the simplex. """ @property @memoize_method @@ -166,12 +177,24 @@ class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase): assert dim2 == dim return result + +class PolynomialEquidistantElementGroup(_MassMatrixQuadratureElementGroup): + """Elemental discretization with a number of nodes matching the number of + polynomials in :math:`P^k`, hence usable for differentiation and + interpolation. Interpolation nodes are present on the boundary of the + simplex. + + .. versionadded:: 2016.1 + """ @property @memoize_method - def weights(self): - return np.dot( - mp.mass_matrix(self.basis(), self.unit_nodes), - np.ones(len(self.basis()))) + def unit_nodes(self): + dim = self.mesh_el_group.dim + result = mp.equidistant_nodes(dim, self.order) + + dim2, nunit_nodes = result.shape + assert dim2 == dim + return result # }}} @@ -208,6 +231,15 @@ class PolynomialWarpAndBlendGroupFactory(OrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = PolynomialWarpAndBlendElementGroup + +class PolynomialEquidistantGroupFactory(OrderBasedGroupFactory): + """ + .. versionadded:: 2016.1 + """ + + mesh_group_class = _MeshSimplexElementGroup + group_class = PolynomialEquidistantElementGroup + # }}} diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 53916c7..5c2b5c1 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -251,7 +251,10 @@ class SimplexElementGroup(MeshElementGroup): raise TypeError("'dim' must be passed " "if 'unit_nodes' is not passed") - unit_nodes = mp.warp_and_blend_nodes(dim, order) + if dim <= 3: + unit_nodes = mp.warp_and_blend_nodes(dim, order) + else: + unit_nodes = mp.equidistant_nodes(dim, order) dims = unit_nodes.shape[0] @@ -286,25 +289,8 @@ class SimplexElementGroup(MeshElementGroup): raise NotImplementedError("dim=%d" % self.dim) def vertex_unit_coordinates(self): - if self.dim == 1: - return np.array([ - [-1], [1] - ], dtype=np.float64) - elif self.dim == 2: - return np.array([ - [-1, -1], - [1, -1], - [-1, 1], - ], dtype=np.float64) - elif self.dim == 3: - return np.array([ - [-1, -1, -1], - [1, -1, -1], - [-1, 1, -1], - [-1, -1, 1], - ], dtype=np.float64) - else: - raise NotImplementedError("dim=%d" % self.dim) + from modepy.tools import unit_vertices + return unit_vertices(self.dim) # }}} @@ -673,10 +659,10 @@ def _test_node_vertex_consistency_simplex(mesh, mgrp): if mgrp.nelements == 0: return True - from modepy.tools import UNIT_VERTICES resampling_mat = mp.resampling_matrix( - mp.simplex_onb(mgrp.dim, mgrp.order), - UNIT_VERTICES[mgrp.dim].T.copy(), mgrp.unit_nodes) + mp.simplex_best_available_basis(mgrp.dim, mgrp.order), + mgrp.vertex_unit_coordinates().T, + mgrp.unit_nodes) # dim, nelments, nnvertices map_vertices = np.einsum( diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 1023b52..6576158 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -234,7 +234,11 @@ def make_group_from_vertices(vertices, vertex_indices, order): dim = nspan_vectors # dim, nunit_nodes - unit_nodes = mp.warp_and_blend_nodes(dim, order) + if dim <= 3: + unit_nodes = mp.warp_and_blend_nodes(dim, order) + else: + unit_nodes = mp.equidistant_nodes(dim, order) + unit_nodes_01 = 0.5 + 0.5*unit_nodes nodes = np.einsum( diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 98d5d3b..2c4a5b4 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -35,7 +35,8 @@ from pyopencl.tools import ( # noqa from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory + PolynomialWarpAndBlendGroupFactory, + PolynomialEquidistantGroupFactory, ) from meshmode.mesh import BTAG_ALL from meshmode.discretization.connection import \ @@ -523,6 +524,56 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False): # }}} +# {{{ sanity check: volume interpolation on scipy/qhull delaunay meshes in nD + +@pytest.mark.parametrize("dim", [2, 3, 4]) +@pytest.mark.parametrize("order", [3]) +def test_sanity_qhull_nd(ctx_getter, dim, order): + pytest.importorskip("scipy") + + logging.basicConfig(level=logging.INFO) + + ctx = ctx_getter() + queue = cl.CommandQueue(ctx) + + from scipy.spatial import Delaunay + verts = np.random.rand(1000, dim) + dtri = Delaunay(verts) + + from meshmode.mesh.io import from_vertices_and_simplices + mesh = from_vertices_and_simplices(dtri.points.T, dtri.simplices, + fix_orientation=True) + + from meshmode.discretization import Discretization + low_discr = Discretization(ctx, mesh, + PolynomialEquidistantGroupFactory(order)) + high_discr = Discretization(ctx, mesh, + PolynomialEquidistantGroupFactory(order+1)) + + from meshmode.discretization.connection import make_same_mesh_connection + cnx = make_same_mesh_connection(high_discr, low_discr) + + def f(x): + return 0.1*cl.clmath.sin(x) + + x_low = low_discr.nodes()[0].with_queue(queue) + f_low = f(x_low) + + x_high = high_discr.nodes()[0].with_queue(queue) + f_high_ref = f(x_high) + + f_high_num = cnx(queue, f_low) + + err = (f_high_ref-f_high_num).get() + + err = la.norm(err, np.inf)/la.norm(f_high_ref.get(), np.inf) + + print(err) + assert err < 1e-2 + +# }}} + + # {{{ sanity checks: ball meshes # python test_meshmode.py 'test_sanity_balls(cl._csc, "disk-radius-1.step", 2, 2, visualize=True)' # noqa -- GitLab