From 7dc2d82321731ea558361f62c19f7da0e4700814 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 24 Oct 2020 14:47:27 -0500 Subject: [PATCH 01/51] refactor tensor products a bit --- meshmode/discretization/poly_element.py | 67 ++++++++++++++++++------- 1 file changed, 49 insertions(+), 18 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 22cd60cd4..64a8034f3 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -346,28 +346,37 @@ def unit_nodes(self): # {{{ concrete element groups for tensor product elements -class LegendreGaussLobattoTensorProductElementGroup(PolynomialElementGroupBase): +class TensorProductElementGroupBase(PolynomialElementGroupBase): def is_orthogonal_basis(self): return True - def basis(self): - from modepy.modes import tensor_product_basis, jacobi + def basis_1d(self): + from modepy.modes import jacobi + from functools import partial + return tuple(partial(jacobi, 0, 0, i) for i in range(self.order + 1)) + + def grad_basis_1d(self): + from modepy.modes import grad_jacobi from functools import partial - return tensor_product_basis( - self.dim, tuple( - partial(jacobi, 0, 0, i) - for i in range(self.order+1))) + return tuple(partial(grad_jacobi, 0, 0, i) for i in range(self.order + 1)) + + def unit_nodes_1d(self): + raise NotImplementedError + + def basis(self): + from modepy.modes import tensor_product_basis + return tensor_product_basis(self.dim, self.basis_1d()) def grad_basis(self): - raise NotImplementedError() + from modepy.modes import grad_tensor_product_basis + return grad_tensor_product_basis(self.dim, + self.basis_1d(), self.grad_basis_1d()) @property @memoize_method def unit_nodes(self): from modepy.nodes import tensor_product_nodes - from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes - return tensor_product_nodes( - self.dim, legendre_gauss_lobatto_nodes(self.order)) + return tensor_product_nodes(self.dim, self.unit_nodes_1d) @memoize_method def from_mesh_interp_matrix(self): @@ -375,17 +384,39 @@ def from_mesh_interp_matrix(self): return mp.resampling_matrix( self.basis(), self.unit_nodes, - meg.unit_nodes) + meg.unit_nodes, + least_squares_ok=False) + + +class InterpolatoryQuadratureTensorProductElementGroup(TensorProductElementGroup): + @memoize_method + def _quadrature_rule_1d(self): + from modepy.quadrature import LegendreGaussQuadrature + return LegendreGaussQuadrature(self.order) + @property + def unit_nodes_1d(self): + return self._quadrature_rule().nodes -class EquidistantTensorProductElementGroup( - LegendreGaussLobattoTensorProductElementGroup): @property @memoize_method - def unit_nodes(self): - from modepy.nodes import tensor_product_nodes, equidistant_nodes - return tensor_product_nodes( - self.dim, equidistant_nodes(1, self.order)[0]) + def weights(self): + from modepy.nodes import tensor_product_nodes + return tensor_product_nodes(self.dim, self._quadrature_rule().weights) + + +class LegendreGaussLobattoTensorProductElementGroup(TensorProductElementGroupBase): + @property + def unit_nodes_1d(self): + from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes + return legendre_gauss_lobatto_nodes(self.order) + + +class EquidistantTensorProductElementGroup(TensorProductElementGroupBase): + @property + def unit_nodes_1d(self): + from modepy.nodes import equidistant_nodes + return equidistant_nodes(1, self.order)[0] # }}} From 759497143037cd78b37cb7c7b38eec23b6c7743b Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 30 Oct 2020 21:15:19 -0500 Subject: [PATCH 02/51] implement face_vertex_indices for tensor products --- meshmode/discretization/poly_element.py | 71 ++++++++++++++++++++----- meshmode/mesh/__init__.py | 44 ++++++++++++--- meshmode/mesh/generation.py | 11 ++-- test/test_meshmode.py | 47 ++++++++++------ 4 files changed, 130 insertions(+), 43 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 64a8034f3..b2b3ef3be 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -33,7 +33,7 @@ __doc__ = """ Group types -^^^^^^^^^^^ +^^^^^^^^^^^^^^^^^^^ .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup @@ -41,6 +41,8 @@ .. autoclass:: PolynomialRecursiveNodesElementGroup .. autoclass:: PolynomialEquidistantSimplexElementGroup .. autoclass:: PolynomialGivenNodesElementGroup + +.. autoclass:: InterpolatoryQuadratureTensorProductElementGroup .. autoclass:: LegendreGaussLobattoTensorProductElementGroup .. autoclass:: EquidistantTensorProductElementGroup @@ -48,12 +50,16 @@ ^^^^^^^^^^^^^^^ .. autoclass:: ElementGroupFactory +.. autoclass:: OrderAndTypeBasedGroupFactory + .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory .. autoclass:: PolynomialWarpAndBlendGroupFactory .. autoclass:: PolynomialRecursiveNodesGroupFactory .. autoclass:: PolynomialEquidistantSimplexGroupFactory .. autoclass:: PolynomialGivenNodesGroupFactory + +.. autoclass:: InterpolatoryQuadratureTensorProductGroupFactory .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ @@ -346,9 +352,8 @@ def unit_nodes(self): # {{{ concrete element groups for tensor product elements -class TensorProductElementGroupBase(PolynomialElementGroupBase): - def is_orthogonal_basis(self): - return True +class _TensorProductElementGroupBase(PolynomialElementGroupBase): + # {{{ def basis_1d(self): from modepy.modes import jacobi @@ -363,6 +368,11 @@ def grad_basis_1d(self): def unit_nodes_1d(self): raise NotImplementedError + # }}} + + def is_orthogonal_basis(self): + return True + def basis(self): from modepy.modes import tensor_product_basis return tensor_product_basis(self.dim, self.basis_1d()) @@ -376,25 +386,32 @@ def grad_basis(self): @memoize_method def unit_nodes(self): from modepy.nodes import tensor_product_nodes - return tensor_product_nodes(self.dim, self.unit_nodes_1d) + return tensor_product_nodes(self.dim, self.unit_nodes_1d()) @memoize_method def from_mesh_interp_matrix(self): + from modepy.modes import legendre_tensor_product_basis meg = self.mesh_el_group return mp.resampling_matrix( - self.basis(), + legendre_tensor_product_basis(self.dim, self.order), self.unit_nodes, - meg.unit_nodes, - least_squares_ok=False) + meg.unit_nodes) + + +class InterpolatoryQuadratureTensorProductElementGroup( + _TensorProductElementGroupBase): + """Elemental discretization supplying a high-order quadrature rule + with a number of nodes matching the number of polynomials in the tensor + product basis, hence usable for differentiation and interpolation. + No interpolation nodes are present on the boundary of the hypercube. + """ -class InterpolatoryQuadratureTensorProductElementGroup(TensorProductElementGroup): @memoize_method def _quadrature_rule_1d(self): from modepy.quadrature import LegendreGaussQuadrature return LegendreGaussQuadrature(self.order) - @property def unit_nodes_1d(self): return self._quadrature_rule().nodes @@ -405,15 +422,30 @@ def weights(self): return tensor_product_nodes(self.dim, self._quadrature_rule().weights) -class LegendreGaussLobattoTensorProductElementGroup(TensorProductElementGroupBase): - @property +class LegendreGaussLobattoTensorProductElementGroup( + _TensorProductElementGroupBase): + """Elemental discretization supplying a high-order quadrature rule + with a number of nodes matching the number of polynomials in the tensor + product basis, hence usable for differentiation and interpolation. + Nodes are present on the boundary of the hypercube. + + Uses :func:`~modepy.quadrature.legendre_gauss_lobatto_nodes`. + """ + def unit_nodes_1d(self): from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes return legendre_gauss_lobatto_nodes(self.order) -class EquidistantTensorProductElementGroup(TensorProductElementGroupBase): - @property +class EquidistantTensorProductElementGroup(_TensorProductElementGroupBase): + """Elemental discretization supplying a high-order quadrature rule + with a number of nodes matching the number of polynomials in the tensor + product basis, hence usable for differentiation and interpolation. + Nodes are present on the boundary of the hypercube. + + Uses :func:`~modepy.equidistant_nodes`. + """ + def unit_nodes_1d(self): from modepy.nodes import equidistant_nodes return equidistant_nodes(1, self.order)[0] @@ -514,9 +546,18 @@ def __call__(self, mesh_el_group, index): return PolynomialGivenNodesElementGroup( mesh_el_group, self.order, self.unit_nodes, index) + # }}} +# {{{ group factories for tensor products + +class InterpolatoryQuadratureTensorProductGroupFactory( + HomogeneousOrderBasedGroupFactory): + mesh_group_class = _MeshTensorProductElementGroup + group_class = InterpolatoryQuadratureTensorProductElementGroup + + class LegendreGaussLobattoTensorProductGroupFactory( HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshTensorProductElementGroup @@ -524,5 +565,7 @@ class LegendreGaussLobattoTensorProductGroupFactory( # }}} +# }}} + # vim: fdm=marker diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 494d03c25..6fd7ddc93 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -366,7 +366,7 @@ def face_vertex_indices(self): (1, 2, 3) ) else: - raise NotImplementedError("dim=%d" % self.dim) + raise NotImplementedError(f"dimension {self.dim}") def vertex_unit_coordinates(self): from modepy.tools import unit_vertices @@ -412,10 +412,33 @@ def __init__(self, order, vertex_indices, nodes, element_nr_base, node_nr_base, unit_nodes) def face_vertex_indices(self): - raise NotImplementedError() + if self.dim == 1: + return ( + (0,), + (1,), + ) + elif self.dim == 2: + return ( + (0, 1), + (3, 2), + (2, 0), + (1, 3), + ) + elif self.dim == 3: + return ( + (0, 2, 1, 3), + (4, 6, 5, 7), + (0, 4, 5, 1), + (2, 6, 7, 3), + (0, 4, 6, 2), + (1, 5, 7, 3), + ) + else: + raise NotImplementedError(f"dimension {self.dim}") def vertex_unit_coordinates(self): - raise NotImplementedError() + from pytools import generate_nonnegative_integer_tuples_below as gnitb + return np.array(list(gnitb(2, self.dim)), dtype=np.float64)*2.0 - 1.0 # }}} @@ -971,15 +994,22 @@ def __ne__(self, other): # {{{ node-vertex consistency test -def _test_node_vertex_consistency_simplex(mesh, mgrp, tol): +def _test_node_vertex_consistency_resampling(mesh, mgrp, tol): if mesh.vertices is None: return True if mgrp.nelements == 0: return True + if isinstance(mgrp, SimplexElementGroup): + basis = mp.simplex_best_available_basis(mgrp.dim, mgrp.order) + elif isinstance(mgrp, TensorProductElementGroup): + basis = mp.legendre_tensor_product_basis(mgrp.dim, mgrp.order) + else: + raise TypeError(f"unsupported group type: {type(mgrp).__name__}") + resampling_mat = mp.resampling_matrix( - mp.simplex_best_available_basis(mgrp.dim, mgrp.order), + basis, mgrp.vertex_unit_coordinates().T, mgrp.unit_nodes) @@ -1013,8 +1043,8 @@ def _test_node_vertex_consistency(mesh, tol): """ for mgrp in mesh.groups: - if isinstance(mgrp, SimplexElementGroup): - assert _test_node_vertex_consistency_simplex(mesh, mgrp, tol) + if isinstance(mgrp, (SimplexElementGroup, TensorProductElementGroup)): + assert _test_node_vertex_consistency_resampling(mesh, mgrp, tol) else: from warnings import warn warn("not implemented: node-vertex consistency check for '%s'" diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index d283c0031..b1ca6ca93 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -767,7 +767,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64, if dim == 1: if mesh_type is not None: - raise ValueError("unsupported mesh_type") + raise ValueError(f"unsupported mesh type: '{mesh_type}'") for i in range(shape[0]-1): # a--b @@ -801,7 +801,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64, pass else: - raise ValueError("unsupported mesh_type") + raise ValueError(f"unsupported mesh type: '{mesh_type}'") for i in range(shape[0]-1): for j in range(shape[1]-1): @@ -957,7 +957,8 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1, on [a,b]. :param boundary_tag_to_face: an optional dictionary for tagging boundaries. See :func:`generate_box_mesh`. - :param mesh_type: See :func:`generate_box_mesh`. + :param group_cls: see :func:`generate_box_mesh`. + :param mesh_type: see :func:`generate_box_mesh`. """ if min(n) < 2: raise ValueError("need at least two points in each direction") @@ -975,7 +976,7 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1, # {{{ generate_warped_rect_mesh -def generate_warped_rect_mesh(dim, order, n): +def generate_warped_rect_mesh(dim, order, n, *, group_cls=None): """Generate a mesh of a warped line/square/cube. Mainly useful for testing functionality with curvilinear meshes. """ @@ -983,7 +984,7 @@ def generate_warped_rect_mesh(dim, order, n): assert dim in [1, 2, 3] mesh = generate_regular_rect_mesh( a=(-0.5,)*dim, b=(0.5,)*dim, - n=(n,)*dim, order=order) + n=(n,)*dim, order=order, group_cls=group_cls) def m(x): result = np.empty_like(x) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cf90befb7..869065bed 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -285,14 +285,16 @@ def test_boundary_tags(): ]) @pytest.mark.parametrize("group_cls", [ SimplexElementGroup, - - # FIXME: Not implemented: TPE.face_vertex_indices - # TensorProductElementGroup + TensorProductElementGroup ]) def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): + if group_cls is TensorProductElementGroup and mesh_type is not None: + pytest.skip("mesh type not supported on tensor product elements") + from meshmode.mesh.generation import generate_regular_rect_mesh from meshmode.mesh import is_boundary_tag_empty from meshmode.mesh import check_bc_coverage + if dim == 1: a = (0,) b = (1,) @@ -317,7 +319,7 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): group_cls=group_cls, mesh_type=mesh_type) - if visualize: + if visualize and dim == 2: from meshmode.mesh.visualization import draw_2d_mesh draw_2d_mesh(mesh, draw_element_numbers=False, draw_vertex_numbers=False) import matplotlib.pyplot as plt @@ -326,8 +328,10 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): # correct answer if dim == 1: num_on_bdy = 1 - else: - num_on_bdy = dim*(dim-1)*(nelem-1)**(dim-1) + elif group_cls is TensorProductElementGroup: + num_on_bdy = dim * (nelem-1)**(dim-1) + elif group_cls is SimplexElementGroup: + num_on_bdy = dim * (dim-1) * (nelem-1)**(dim-1) assert not is_boundary_tag_empty(mesh, "btag_test_1") assert not is_boundary_tag_empty(mesh, "btag_test_2") @@ -358,7 +362,6 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): raise ValueError("%i marked on custom boundary 2, should be %i" % (num_marked_bdy_2, num_on_bdy)) - # }}} @@ -368,7 +371,10 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), - #partial(PolynomialRecursiveNodesGroupFactory, family="gc"), + # partial(PolynomialRecursiveNodesGroupFactory, family="gc"), + + # FIXME: need support in make_face_restriction + # LegendreGaussLobattoTensorProductGroupFactory, ]) @pytest.mark.parametrize("boundary_tag", [ BTAG_ALL, @@ -383,8 +389,17 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): @pytest.mark.parametrize("per_face_groups", [False, True]) def test_boundary_interpolation(actx_factory, group_factory, boundary_tag, mesh_name, dim, mesh_pars, per_face_groups): + if (issubclass(group_factory, LegendreGaussLobattoTensorProductGroupFactory) + and mesh_name == "blob"): + pytest.skip("tensor products not implemented on blobs") + actx = actx_factory() + if issubclass(group_factory, LegendreGaussLobattoTensorProductGroupFactory): + group_cls = TensorProductElementGroup + else: + group_cls = SimplexElementGroup + from meshmode.discretization import Discretization from meshmode.discretization.connection import ( make_face_restriction, check_connection) @@ -420,7 +435,8 @@ def f(x): force_ambient_dim=2) elif mesh_name == "warp": from meshmode.mesh.generation import generate_warped_rect_mesh - mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par) + mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par, + group_cls=group_cls) h = 1/mesh_par else: @@ -586,7 +602,7 @@ def f(x): ("segment", 1, [8, 16, 32]), ("blob", 2, [1e-1, 8e-2, 5e-2]), ("warp", 2, [3, 5, 7]), - ("warp", 3, [3, 5]), + ("warp", 3, [3, 5]) ]) def test_opposite_face_interpolation(actx_factory, group_factory, mesh_name, dim, mesh_pars): @@ -642,8 +658,7 @@ def f(x): # }}} - vol_discr = Discretization(actx, mesh, - group_factory(order)) + vol_discr = Discretization(actx, mesh, group_factory(order)) print("h=%s -> %d elements" % ( h, sum(mgrp.nelements for mgrp in mesh.groups))) @@ -767,9 +782,7 @@ def test_orientation_3d(actx_factory, what, mesh_gen_func, visualize=False): @pytest.mark.parametrize("group_cls", [ SimplexElementGroup, - - # NOTE: needs TensorProductElementGroup.face_vertex_indices - # TensorProductElementGroup + TensorProductElementGroup ]) def test_merge_and_map(actx_factory, group_cls, visualize=False): from meshmode.mesh.io import generate_gmsh, FileSource @@ -1242,7 +1255,7 @@ def no_test_quad_mesh_3d(): # {{{ test_quad_single_element -def test_quad_single_element(): +def test_quad_single_element(visualize=False): from meshmode.mesh.generation import make_group_from_vertices vertices = np.array([ @@ -1257,7 +1270,7 @@ def test_quad_single_element(): 30, group_cls=TensorProductElementGroup) Mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None) - if 0: + if visualize: import matplotlib.pyplot as plt plt.plot( mg.nodes[0].reshape(-1), From 0164616aeac7bb09e802d26408ed5dd53f0f02c8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 30 Oct 2020 21:21:01 -0500 Subject: [PATCH 03/51] fix doc header --- meshmode/discretization/poly_element.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index b2b3ef3be..b1fccdca5 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -33,7 +33,7 @@ __doc__ = """ Group types -^^^^^^^^^^^^^^^^^^^ +^^^^^^^^^^^ .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup @@ -429,7 +429,7 @@ class LegendreGaussLobattoTensorProductElementGroup( product basis, hence usable for differentiation and interpolation. Nodes are present on the boundary of the hypercube. - Uses :func:`~modepy.quadrature.legendre_gauss_lobatto_nodes`. + Uses :func:`~modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes`. """ def unit_nodes_1d(self): From fa414dbfa16bc1ea5986c99afbf7048005e32b96 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 31 Oct 2020 14:43:32 -0500 Subject: [PATCH 04/51] remove inappropriate uses of issubclass --- test/test_meshmode.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 869065bed..cad77eb3d 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -389,13 +389,13 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): @pytest.mark.parametrize("per_face_groups", [False, True]) def test_boundary_interpolation(actx_factory, group_factory, boundary_tag, mesh_name, dim, mesh_pars, per_face_groups): - if (issubclass(group_factory, LegendreGaussLobattoTensorProductGroupFactory) + if (group_factory is LegendreGaussLobattoTensorProductGroupFactory and mesh_name == "blob"): pytest.skip("tensor products not implemented on blobs") actx = actx_factory() - if issubclass(group_factory, LegendreGaussLobattoTensorProductGroupFactory): + if group_factory is LegendreGaussLobattoTensorProductGroupFactory: group_cls = TensorProductElementGroup else: group_cls = SimplexElementGroup From 554214b7e0256e89d3d60d85fc380819411b6ce7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 2 Nov 2020 14:34:59 -0600 Subject: [PATCH 05/51] Implement tensor product face restrictions (with voodoo element basis reversal, FIXME) --- meshmode/discretization/connection/face.py | 86 ++++++++++++++++++---- meshmode/mesh/__init__.py | 38 ++++++++-- test/test_meshmode.py | 20 ++++- 3 files changed, 118 insertions(+), 26 deletions(-) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 9eddc07c7..a2cfe9226 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -75,7 +75,9 @@ def _build_boundary_connection(actx, vol_discr, bdry_discr, connection_data, data = connection_data[igrp, face_id] bdry_unit_nodes_01 = (bdry_grp.unit_nodes + 1)*0.5 - result_unit_nodes = (np.dot(data.A, bdry_unit_nodes_01).T + data.b).T + result_unit_nodes = ( + np.dot(data.face_basis.T, bdry_unit_nodes_01).T + + data.face_offset).T batches.append( InterpolationBatch( @@ -208,7 +210,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # }}} - from meshmode.mesh import Mesh, SimplexElementGroup + from meshmode.mesh import Mesh, SimplexElementGroup, TensorProductElementGroup bdry_mesh_groups = [] connection_data = {} @@ -220,9 +222,10 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, mgrp = grp.mesh_el_group - if not isinstance(mgrp, SimplexElementGroup): + if not isinstance(mgrp, (SimplexElementGroup, TensorProductElementGroup)): raise NotImplementedError("can only take boundary of " - "SimplexElementGroup-based meshes") + "meshes based on SimplexElementGroup and " + "TensorProductElementGroup") # {{{ pull together per-group face lists @@ -282,14 +285,31 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, else: ngroup_bdry_elements = len(group_boundary_faces) + # make up some not-terrible nodes for the boundary Mesh + if isinstance(mgrp, SimplexElementGroup): + bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order) + vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order) + + nbdry_el_vertices = (mgrp.dim-1) + 1 + + elif isinstance(mgrp, TensorProductElementGroup): + bdry_unit_nodes = \ + mp.legendre_gauss_lobatto_tensor_product_nodes( + mgrp.dim-1, mgrp.order) + vol_basis = mp.legendre_tensor_product_basis( + mgrp.dim, mgrp.order) + + nbdry_el_vertices = 2**(mgrp.dim-1) + + else: + assert False + vertex_indices = np.empty( - (ngroup_bdry_elements, mgrp.dim+1-1), + (ngroup_bdry_elements, nbdry_el_vertices), mgrp.vertex_indices.dtype) - bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order) bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5 - vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order) nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1] nodes = np.empty( (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes), @@ -310,16 +330,52 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, face_vertex_unit_coordinates = \ grp_vertex_unit_coordinates[loc_face_vertices] - # Find A, b such that A [e_1 e_2] + b = [r_1 r_2] + # Find face_basis, face_offset such that + # face_basis.T @ [e_1 e_2] + face_offset == [r_1 r_2]. + # # (Notation assumes that the volume is 3D and the face is 2D. # Code does not.) - b = face_vertex_unit_coordinates[0] - A = ( # noqa + face_offset = face_vertex_unit_coordinates[0] + face_basis = ( face_vertex_unit_coordinates[1:] - - face_vertex_unit_coordinates[0]).T + - face_offset) + + if isinstance(mgrp, SimplexElementGroup): + # no further action required + pass + + elif isinstance(mgrp, TensorProductElementGroup): + if mgrp.dim <= 2: + # faces are simplices, no action required + pass + + elif mgrp.dim == 3: + # Drop the trailing basis vector, so that the convex + # combination of the remaining two is a subset of the face. + reduced_face_basis = face_basis[:-1] + + # FIXME: This makes the tests pass, but I do not understand + # why this reversal is necessary. + reduced_face_basis = reduced_face_basis[::-1] + + # assert that all basis vectors are axis-aligned + for bvec in reduced_face_basis: + assert sum(abs(entry) > 1e-13 for entry in bvec) == 1 + + face_basis = reduced_face_basis + del reduced_face_basis + + else: + raise NotImplementedError( + f"face restriction of a {mgrp.dim}D hypercube") + + else: + assert False - face_unit_nodes = (np.dot(A, bdry_unit_nodes_01).T + b).T + face_unit_nodes = ( + np.dot(face_basis.T, bdry_unit_nodes_01).T + + face_offset).T resampling_mat = mp.resampling_matrix( vol_basis, @@ -346,14 +402,14 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, connection_data[igrp, face_id] = _ConnectionBatchData( group_source_element_indices=batch_boundary_el_numbers_in_grp, group_target_element_indices=new_el_numbers, - A=A, - b=b, + face_basis=face_basis, + face_offset=face_offset, ) is_last_face = face_id + 1 == mgrp.nfaces if per_face_groups or is_last_face: - bdry_mesh_group = SimplexElementGroup( + bdry_mesh_group = type(mgrp)( mgrp.order, vertex_indices, nodes, unit_nodes=bdry_unit_nodes) bdry_mesh_groups.append(bdry_mesh_group) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 6fd7ddc93..9e266dfbd 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -385,7 +385,7 @@ class TensorProductElementGroup(MeshElementGroup): def __init__(self, order, vertex_indices, nodes, element_nr_base=None, node_nr_base=None, unit_nodes=None): - """ + r""" :arg order: the maximum total degree used for interpolation. :arg nodes: ``[ambient_dim, nelements, nunit_nodes]`` The nodes are assumed to be mapped versions of *unit_nodes*. @@ -393,6 +393,23 @@ def __init__(self, order, vertex_indices, nodes, The unit nodes of which *nodes* is a mapped version. + The vertex order in 3D is as follows :: + + + ^ s + | 6-----7 + | /| /| + |/ | / | + 2-----3 | + | | | | + | 4--|--5 + | / | / + |/ |/ + 0-----1---->r + + For general dimensions, follow binary upward counting: + 000, 001, 010, 011, ... + Do not supply *element_nr_base* and *node_nr_base*, they will be automatically assigned. """ @@ -426,12 +443,19 @@ def face_vertex_indices(self): ) elif self.dim == 3: return ( - (0, 2, 1, 3), - (4, 6, 5, 7), - (0, 4, 5, 1), - (2, 6, 7, 3), - (0, 4, 6, 2), - (1, 5, 7, 3), + # binary is tsr + + # rs-aligned + (0b000, 0b001, 0b010, 0b011,), + (0b100, 0b101, 0b110, 0b111,), + + # st-aligned + (0b000, 0b010, 0b100, 0b110,), + (0b001, 0b011, 0b101, 0b111,), + + # rt-aligned + (0b000, 0b001, 0b100, 0b101,), + (0b010, 0b011, 0b110, 0b111,), ) else: raise NotImplementedError(f"dimension {self.dim}") diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cad77eb3d..144626743 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -371,10 +371,11 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), + + # Redundant, no information gain. # partial(PolynomialRecursiveNodesGroupFactory, family="gc"), - # FIXME: need support in make_face_restriction - # LegendreGaussLobattoTensorProductGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory, ]) @pytest.mark.parametrize("boundary_tag", [ BTAG_ALL, @@ -383,10 +384,14 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ ("blob", 2, ["8e-2", "6e-2", "4e-2"]), + + # If "warp" works, "rect" likely does, too. + # ("rect", 3, [10, 20, 30]), + ("warp", 2, [10, 20, 30]), ("warp", 3, [10, 20, 30]), ]) -@pytest.mark.parametrize("per_face_groups", [False, True]) +@pytest.mark.parametrize("per_face_groups", [True, False]) def test_boundary_interpolation(actx_factory, group_factory, boundary_tag, mesh_name, dim, mesh_pars, per_face_groups): if (group_factory is LegendreGaussLobattoTensorProductGroupFactory @@ -439,6 +444,13 @@ def f(x): group_cls=group_cls) h = 1/mesh_par + + elif mesh_name == "rect": + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh(a=(0,)*dim, b=(1,)*dim, + order=4, n=(mesh_par,)*dim, group_cls=group_cls) + + h = 1/mesh_par else: raise ValueError("mesh_name not recognized") @@ -478,7 +490,7 @@ def f(x): print(eoc_rec) assert ( eoc_rec.order_estimate() >= order-order_slack - or eoc_rec.max_error() < 3e-14) + or eoc_rec.max_error() < 1e-13) # }}} From 8dac59e1cd2bdc87f0fd96f68a6a0eb6694199fe Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 3 Nov 2020 14:28:06 -0600 Subject: [PATCH 06/51] improve docs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/discretization/poly_element.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index b1fccdca5..234c01a29 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -353,7 +353,7 @@ def unit_nodes(self): # {{{ concrete element groups for tensor product elements class _TensorProductElementGroupBase(PolynomialElementGroupBase): - # {{{ + # {{{ 1D basis def basis_1d(self): from modepy.modes import jacobi @@ -427,7 +427,7 @@ class LegendreGaussLobattoTensorProductElementGroup( """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. - Nodes are present on the boundary of the hypercube. + Nodes sufficient for unisolvency are present on the boundary of the hypercube. Uses :func:`~modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes`. """ @@ -441,7 +441,7 @@ class EquidistantTensorProductElementGroup(_TensorProductElementGroupBase): """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. - Nodes are present on the boundary of the hypercube. + Nodes sufficient for unisolvency are present on the boundary of the hypercube. Uses :func:`~modepy.equidistant_nodes`. """ From fa013344ff4a70ba059dfcb62198d827eb908970 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 3 Nov 2020 14:51:49 -0600 Subject: [PATCH 07/51] underscore 1d tensor product methods --- meshmode/discretization/poly_element.py | 33 ++++++++++++++----------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 234c01a29..ee9949798 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -355,17 +355,17 @@ def unit_nodes(self): class _TensorProductElementGroupBase(PolynomialElementGroupBase): # {{{ 1D basis - def basis_1d(self): + def _basis_1d(self): from modepy.modes import jacobi from functools import partial return tuple(partial(jacobi, 0, 0, i) for i in range(self.order + 1)) - def grad_basis_1d(self): + def _grad_basis_1d(self): from modepy.modes import grad_jacobi from functools import partial return tuple(partial(grad_jacobi, 0, 0, i) for i in range(self.order + 1)) - def unit_nodes_1d(self): + def _unit_nodes_1d(self): raise NotImplementedError # }}} @@ -375,18 +375,18 @@ def is_orthogonal_basis(self): def basis(self): from modepy.modes import tensor_product_basis - return tensor_product_basis(self.dim, self.basis_1d()) + return tensor_product_basis(self.dim, self._basis_1d()) def grad_basis(self): from modepy.modes import grad_tensor_product_basis return grad_tensor_product_basis(self.dim, - self.basis_1d(), self.grad_basis_1d()) + self._basis_1d(), self._grad_basis_1d()) @property @memoize_method def unit_nodes(self): from modepy.nodes import tensor_product_nodes - return tensor_product_nodes(self.dim, self.unit_nodes_1d()) + return tensor_product_nodes(self.dim, self._unit_nodes_1d()) @memoize_method def from_mesh_interp_matrix(self): @@ -409,17 +409,20 @@ class InterpolatoryQuadratureTensorProductElementGroup( @memoize_method def _quadrature_rule_1d(self): - from modepy.quadrature import LegendreGaussQuadrature - return LegendreGaussQuadrature(self.order) + import modepy as mp + return mp.LegendreGaussQuadrature(self.order) - def unit_nodes_1d(self): - return self._quadrature_rule().nodes + def _unit_nodes_1d(self): + return self._quadrature_rule_1d().nodes @property - @memoize_method def weights(self): - from modepy.nodes import tensor_product_nodes - return tensor_product_nodes(self.dim, self._quadrature_rule().weights) + from itertools import product + weights = self._quadrature_rule_1d().weights + return np.fromiter( + (np.prod(w) for w in product(weights, repeat=self.dim)), + dtype=np.float, + count=self.order**self.dim) class LegendreGaussLobattoTensorProductElementGroup( @@ -432,7 +435,7 @@ class LegendreGaussLobattoTensorProductElementGroup( Uses :func:`~modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes`. """ - def unit_nodes_1d(self): + def _unit_nodes_1d(self): from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes return legendre_gauss_lobatto_nodes(self.order) @@ -446,7 +449,7 @@ class EquidistantTensorProductElementGroup(_TensorProductElementGroupBase): Uses :func:`~modepy.equidistant_nodes`. """ - def unit_nodes_1d(self): + def _unit_nodes_1d(self): from modepy.nodes import equidistant_nodes return equidistant_nodes(1, self.order)[0] From 07f315dbfc487f73d7ebb77359e9ac6037b58106 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 3 Nov 2020 15:30:58 -0600 Subject: [PATCH 08/51] expand more tests to tensor products --- test/test_meshmode.py | 50 ++++++++++++++++++++++++++++++++----------- 1 file changed, 38 insertions(+), 12 deletions(-) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 144626743..975ed6263 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -497,16 +497,30 @@ def f(x): # {{{ boundary-to-all-faces connecttion +@pytest.mark.parametrize("group_factory", [ + InterpolatoryQuadratureSimplexGroupFactory, + PolynomialWarpAndBlendGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory, + ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ ("blob", 2, [1e-1, 8e-2, 5e-2]), ("warp", 2, [10, 20, 30]), ("warp", 3, [10, 20, 30]), ]) @pytest.mark.parametrize("per_face_groups", [False, True]) -def test_all_faces_interpolation(actx_factory, mesh_name, dim, mesh_pars, - per_face_groups): +def test_all_faces_interpolation(actx_factory, group_factory, + mesh_name, dim, mesh_pars, per_face_groups): + if (group_factory is LegendreGaussLobattoTensorProductGroupFactory + and mesh_name == "blob"): + pytest.skip("tensor products not implemented on blobs") + actx = actx_factory() + if group_factory is LegendreGaussLobattoTensorProductGroupFactory: + group_cls = TensorProductElementGroup + else: + group_cls = SimplexElementGroup + from meshmode.discretization import Discretization from meshmode.discretization.connection import ( make_face_restriction, make_face_to_all_faces_embedding, @@ -540,7 +554,8 @@ def f(x): print("END GEN") elif mesh_name == "warp": from meshmode.mesh.generation import generate_warped_rect_mesh - mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par) + mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par, + group_cls=group_cls) h = 1/mesh_par else: @@ -548,13 +563,12 @@ def f(x): # }}} - vol_discr = Discretization(actx, mesh, - PolynomialWarpAndBlendGroupFactory(order)) + vol_discr = Discretization(actx, mesh, group_factory(order)) print("h=%s -> %d elements" % ( h, sum(mgrp.nelements for mgrp in mesh.groups))) all_face_bdry_connection = make_face_restriction( - actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order), + actx, vol_discr, group_factory(order), FACE_RESTR_ALL, per_face_groups=per_face_groups) all_face_bdry_discr = all_face_bdry_connection.to_discr @@ -579,7 +593,7 @@ def f(x): FACE_RESTR_INTERIOR, ]: bdry_connection = make_face_restriction( - actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order), + actx, vol_discr, group_factory(order), boundary_tag, per_face_groups=per_face_groups) bdry_discr = bdry_connection.to_discr @@ -608,19 +622,29 @@ def f(x): @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory + PolynomialWarpAndBlendGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory, ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ ("segment", 1, [8, 16, 32]), ("blob", 2, [1e-1, 8e-2, 5e-2]), ("warp", 2, [3, 5, 7]), - ("warp", 3, [3, 5]) + ("warp", 3, [5, 7]) ]) def test_opposite_face_interpolation(actx_factory, group_factory, mesh_name, dim, mesh_pars): + if (group_factory is LegendreGaussLobattoTensorProductGroupFactory + and mesh_name in ["segment", "blob"]): + pytest.skip("tensor products not implemented on blobs") + logging.basicConfig(level=logging.INFO) actx = actx_factory() + if group_factory is LegendreGaussLobattoTensorProductGroupFactory: + group_cls = TensorProductElementGroup + else: + group_cls = SimplexElementGroup + from meshmode.discretization import Discretization from meshmode.discretization.connection import ( make_face_restriction, make_opposite_face_connection, @@ -643,7 +667,8 @@ def f(x): from meshmode.mesh.generation import generate_box_mesh mesh = generate_box_mesh( [np.linspace(-0.5, 0.5, mesh_par)], - order=order) + order=order, + group_cls=group_cls) h = 1.0 / mesh_par elif mesh_name == "blob": assert dim == 2 @@ -662,7 +687,8 @@ def f(x): print("END GEN") elif mesh_name == "warp": from meshmode.mesh.generation import generate_warped_rect_mesh - mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par) + mesh = generate_warped_rect_mesh(dim, order=order, n=mesh_par, + group_cls=group_cls) h = 1/mesh_par else: @@ -692,7 +718,7 @@ def f(x): print(eoc_rec) assert ( eoc_rec.order_estimate() >= order-0.5 - or eoc_rec.max_error() < 1e-13) + or eoc_rec.max_error() < 1.5e-13) # }}} From 4dcb4c1ff600251c31db94eeba8c01a7abddb74f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 3 Nov 2020 15:40:10 -0600 Subject: [PATCH 09/51] make ascii cube slightly larger --- meshmode/mesh/__init__.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 9e266dfbd..cc151bdf0 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -397,15 +397,17 @@ def __init__(self, order, vertex_indices, nodes, ^ s - | 6-----7 - | /| /| - |/ | / | - 2-----3 | - | | | | - | 4--|--5 - | / | / - |/ |/ - 0-----1---->r + | 6-----------7 + | /| /| + | / | ^ / | + |/ | / t / | + 2-----------3 | + | |/ | | + | 4-------|---5 + | / | / + | / | / + |/ |/ + 0-----------1---->r For general dimensions, follow binary upward counting: 000, 001, 010, 011, ... From 6bc483dd07f00e9020290d120b392535d7998e0c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 4 Nov 2020 19:21:39 -0600 Subject: [PATCH 10/51] update comment --- meshmode/discretization/connection/face.py | 17 +- test/cubed-cube.msh | 864 --------------------- test/test_meshmode.py | 2 +- 3 files changed, 6 insertions(+), 877 deletions(-) delete mode 100644 test/cubed-cube.msh diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index a2cfe9226..7cb7a33ea 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -142,9 +142,9 @@ def _get_face_vertices(mesh, boundary_tag): # }}} else: - # For FRESTR_INTERIOR_FACES, this is likely every vertex in the book. + # For FACE_RESTR_INTERIOR, this is likely every vertex in the book. # Don't ever bother trying to cut the list down. - # For FRESTR_ALL_FACES, it literally is every single vertex. + # For FACE_RESTR_ALL, it literally is every single vertex. return np.arange(mesh.nvertices, dtype=np.intp) @@ -338,7 +338,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, face_offset = face_vertex_unit_coordinates[0] face_basis = ( - face_vertex_unit_coordinates[1:] + face_vertex_unit_coordinates[1:mgrp.dim] - face_offset) if isinstance(mgrp, SimplexElementGroup): @@ -351,21 +351,14 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, pass elif mgrp.dim == 3: - # Drop the trailing basis vector, so that the convex - # combination of the remaining two is a subset of the face. - reduced_face_basis = face_basis[:-1] - # FIXME: This makes the tests pass, but I do not understand # why this reversal is necessary. - reduced_face_basis = reduced_face_basis[::-1] + face_basis = face_basis[::-1] # assert that all basis vectors are axis-aligned - for bvec in reduced_face_basis: + for bvec in face_basis: assert sum(abs(entry) > 1e-13 for entry in bvec) == 1 - face_basis = reduced_face_basis - del reduced_face_basis - else: raise NotImplementedError( f"face restriction of a {mgrp.dim}D hypercube") diff --git a/test/cubed-cube.msh b/test/cubed-cube.msh deleted file mode 100644 index dde325a0a..000000000 --- a/test/cubed-cube.msh +++ /dev/null @@ -1,864 +0,0 @@ -$MeshFormat -2.2 0 8 -$EndMeshFormat -$Nodes -343 -1 0 0 0 -2 0 0 10 -3 0 10 0 -4 0 10 10 -5 10 0 0 -6 10 0 10 -7 10 10 0 -8 10 10 10 -9 0 0 1.666666666666672 -10 0 0 3.333333333333344 -11 0 0 5.000000000000009 -12 0 0 6.666666666666672 -13 0 0 8.333333333333336 -14 0 1.666666666666672 0 -15 0 3.333333333333344 0 -16 0 5.000000000000009 0 -17 0 6.666666666666672 0 -18 0 8.333333333333336 0 -19 0 10 1.666666666666672 -20 0 10 3.333333333333344 -21 0 10 5.000000000000009 -22 0 10 6.666666666666672 -23 0 10 8.333333333333336 -24 0 1.666666666666672 10 -25 0 3.333333333333344 10 -26 0 5.000000000000009 10 -27 0 6.666666666666672 10 -28 0 8.333333333333336 10 -29 10 0 1.666666666666672 -30 10 0 3.333333333333344 -31 10 0 5.000000000000009 -32 10 0 6.666666666666672 -33 10 0 8.333333333333336 -34 10 1.666666666666672 0 -35 10 3.333333333333344 0 -36 10 5.000000000000009 0 -37 10 6.666666666666672 0 -38 10 8.333333333333336 0 -39 10 10 1.666666666666672 -40 10 10 3.333333333333344 -41 10 10 5.000000000000009 -42 10 10 6.666666666666672 -43 10 10 8.333333333333336 -44 10 1.666666666666672 10 -45 10 3.333333333333344 10 -46 10 5.000000000000009 10 -47 10 6.666666666666672 10 -48 10 8.333333333333336 10 -49 1.666666666666672 0 0 -50 3.333333333333344 0 0 -51 5.000000000000009 0 0 -52 6.666666666666672 0 0 -53 8.333333333333336 0 0 -54 1.666666666666672 0 10 -55 3.333333333333344 0 10 -56 5.000000000000009 0 10 -57 6.666666666666672 0 10 -58 8.333333333333336 0 10 -59 1.666666666666672 10 0 -60 3.333333333333344 10 0 -61 5.000000000000009 10 0 -62 6.666666666666672 10 0 -63 8.333333333333336 10 0 -64 1.666666666666672 10 10 -65 3.333333333333344 10 10 -66 5.000000000000009 10 10 -67 6.666666666666672 10 10 -68 8.333333333333336 10 10 -69 0 1.666666666666668 3.333333333333338 -70 0 8.333333333333336 5.000000000000004 -71 0 5.000000000000005 8.333333333333336 -72 0 5.000000000000006 1.666666666666667 -73 0 8.333333333333334 3.333333333333338 -74 0 3.333333333333338 8.333333333333334 -75 0 8.333333333333336 8.333333333333336 -76 0 1.666666666666668 5.000000000000005 -77 0 6.666666666666671 1.666666666666668 -78 0 1.666666666666668 6.666666666666671 -79 0 1.66666666666667 1.666666666666669 -80 0 1.66666666666667 8.333333333333334 -81 0 8.333333333333334 1.66666666666667 -82 0 5.000000000000003 6.666666666666669 -83 0 3.333333333333337 3.333333333333335 -84 0 3.333333333333335 6.666666666666669 -85 0 6.666666666666669 5.000000000000002 -86 0 5.000000000000003 3.333333333333334 -87 0 3.333333333333339 1.666666666666668 -88 0 6.666666666666671 8.333333333333336 -89 0 8.333333333333336 6.666666666666671 -90 0 6.666666666666668 3.333333333333335 -91 0 5.000000000000002 5.000000000000002 -92 0 3.333333333333335 5.000000000000002 -93 0 6.66666666666667 6.666666666666671 -94 10 1.666666666666668 3.333333333333338 -95 10 8.333333333333336 5.000000000000004 -96 10 5.000000000000005 8.333333333333336 -97 10 5.000000000000004 1.666666666666668 -98 10 8.333333333333337 8.333333333333336 -99 10 8.333333333333334 3.333333333333338 -100 10 3.333333333333338 8.333333333333334 -101 10 1.666666666666667 5.000000000000006 -102 10 6.666666666666671 1.666666666666668 -103 10 1.666666666666668 6.666666666666671 -104 10 1.666666666666669 1.66666666666667 -105 10 1.666666666666669 8.333333333333334 -106 10 8.333333333333334 1.666666666666669 -107 10 6.666666666666669 5.000000000000004 -108 10 3.333333333333333 5.000000000000004 -109 10 6.666666666666668 3.333333333333337 -110 10 3.333333333333334 3.333333333333337 -111 10 3.333333333333337 1.666666666666668 -112 10 6.666666666666671 8.333333333333336 -113 10 8.333333333333337 6.666666666666671 -114 10 5.000000000000003 6.666666666666671 -115 10 3.333333333333335 6.66666666666667 -116 10 5.000000000000001 5.000000000000003 -117 10 5.000000000000002 3.333333333333337 -118 10 6.666666666666671 6.666666666666671 -119 4.963214356275155 0 1.668784957944797 -120 3.317582237776889 0 8.331760196778665 -121 8.315733273001639 0 5.000155748895926 -122 1.623165746041504 0 4.99391383873798 -123 1.636937952053827 0 1.662777951906097 -124 6.654194529315052 0 8.332831488754294 -125 8.323518108830998 0 1.667524092257246 -126 1.656504071342595 0 8.331974911699989 -127 8.323288233378031 0 6.666461769940527 -128 8.327361873953064 0 8.333156029608485 -129 1.625841868965969 0 3.328919181967938 -130 8.316557036413666 0 3.334049509706089 -131 6.638537756742497 0 6.665974027330344 -132 1.642552721132448 0 6.663695878577564 -133 4.982918784821877 0 8.332255412505511 -134 6.64528602688278 0 1.670374112938904 -135 3.278505243836852 0 3.33055000927853 -136 6.62594906800776 0 4.999916920658317 -137 4.937771380378987 0 3.33355951443616 -138 6.62749052424298 0 3.334574749423636 -139 3.286658535168235 0 6.66196704619375 -140 4.958530187990556 0 6.664632353924778 -141 4.943783761023866 0 4.998321579025965 -142 3.296459439269126 0 1.665682739796669 -143 3.270150200748045 0 4.995167900443951 -144 8.319477494274667 10 3.333169972868349 -145 1.642401289952334 10 4.995994122053824 -146 4.990499273962526 10 8.332390464609082 -147 4.974718901986735 10 1.665089100516185 -148 1.650806613158973 10 1.663879441264121 -149 1.666666666666671 10 8.333333333333332 -150 8.326507895108154 10 6.666387147756299 -151 8.324217698423221 10 1.666577446189837 -152 8.330421346341454 10 8.33317246986971 -153 1.637778017620791 10 3.327652077469392 -154 3.324767314036656 10 8.332279894336549 -155 6.64477579136364 10 1.666338457530614 -156 6.660165660811743 10 8.332843092640404 -157 8.321592194848483 10 4.999653680812338 -158 1.651821825244905 10 6.664257878478076 -159 3.28326020484872 10 3.326681050109922 -160 6.634104856166324 10 3.332502948194829 -161 4.97579818627555 10 6.664430588732456 -162 3.284967251501693 10 4.993482825700283 -163 6.641946928853725 10 4.998695732656651 -164 3.308721581528395 10 6.663488840038644 -165 6.651358102186887 10 6.665729135182278 -166 4.956018879142214 10 3.330340837991536 -167 4.961068569266943 10 4.996722349190394 -168 3.306257787745393 10 1.663550909007059 -169 8.332031501668784 3.327655810609839 0 -170 1.661583967397951 4.985925016393396 0 -171 4.99888149235128 8.328781672181401 0 -172 4.993595831504386 1.645761942295224 0 -173 1.660497912810677 1.651563432065564 0 -174 8.332958056923808 8.3315936822691 0 -175 1.665868128287409 8.330677281028967 0 -176 8.333333333333337 1.666666666666669 0 -177 8.332074246359143 4.994289507172165 0 -178 6.662704487145006 3.317607957790754 0 -179 3.323175186310005 1.639797188626136 0 -180 3.332149742463294 8.329006867792472 0 -181 8.332511595730896 6.662870524665738 0 -182 6.664258683158848 1.657777476463054 0 -183 6.665924415024502 8.330016659525276 0 -184 1.664258409514333 6.659472931888667 0 -185 4.991559204424686 3.303632966319623 0 -186 4.996765277527716 6.653675864145596 0 -187 3.32335131306381 3.303903840008758 0 -188 3.325134218047016 4.973125396957162 0 -189 3.329655718808052 6.653762669731537 0 -190 6.663143942894656 4.984579210873332 0 -191 6.66479451352513 6.658229402591791 0 -192 4.992487513742513 4.968476001600282 0 -193 1.660403493758395 3.316764122600103 0 -194 3.332833974604121 8.330588915852831 10 -195 4.990396923629646 1.639192438518225 10 -196 8.331700204977253 4.99348613953155 10 -197 1.66359950713175 4.986809017439535 10 -198 1.659476797213571 1.641166790796788 10 -199 8.331602596644176 1.661100149132812 10 -200 6.666262451736792 8.331036472916939 10 -201 1.666360512036942 8.331775733188479 10 -202 8.333176767431947 8.332461623855972 10 -203 3.324752363659667 1.640654676620429 10 -204 8.331143154549407 3.325601742294733 10 -205 8.332585110528434 6.663185298000827 10 -206 1.665690323527357 6.661884085335791 10 -207 4.999238771603533 8.329411344880642 10 -208 6.660931836157782 1.649309184852187 10 -209 3.322620058061121 3.296169448190562 10 -210 3.331386262994759 6.656799061971727 10 -211 6.662741535306179 4.984497939484297 10 -212 4.988687405622154 3.294823441600734 10 -213 6.660538599871673 3.312736539238734 10 -214 4.997593572036546 6.655043997091015 10 -215 6.66508040097065 6.65903418390889 10 -216 3.327865240662655 4.977063523995539 10 -217 4.993810674308104 4.973188365824303 10 -218 1.660298963080085 3.305559657733506 10 -219 3.323175186310005 1.639797188626136 1.666666665894576 -220 8.333333333333337 1.666666666666669 1.666666666159922 -221 8.333333332883761 8.333333333333336 5.000000000000004 -222 1.666666667399848 6.666666666666669 5.000000000000002 -223 8.333333335154675 3.333333333333338 8.333333333333334 -224 1.666666665276616 5.000000000000006 1.666666666666667 -225 8.333333332534451 8.333333333333337 8.333333333333336 -226 1.66666666735759 5.000000000000002 5.000000000000002 -227 3.325134218047016 4.973125396957162 1.666666665215716 -228 6.66479451352513 6.658229402591791 1.666666668556142 -229 4.992487513742513 4.968476001600282 1.666666662890887 -230 1.666666667636574 6.66666666666667 6.666666666666671 -231 8.333333332812531 8.333333333333334 1.666666666666669 -232 8.333333334638025 5.000000000000003 6.666666666666671 -233 6.666262451736792 8.331036472916939 8.33333333154998 -234 3.317582237776889 1.666666666712437 8.331760196778665 -235 4.963214356275155 1.666666665965598 1.668784957944797 -236 1.666666666034267 3.333333333333335 6.666666666666669 -237 6.66508040097065 6.65903418390889 8.333333333387067 -238 6.634104856166324 8.333333334715782 3.332502948194829 -239 1.666666665693501 6.666666666666671 8.333333333333336 -240 8.332074246359143 4.994289507172165 1.666666666314359 -241 1.666666666304344 6.666666666666668 3.333333333333335 -242 3.284967251501693 8.333333334781612 4.993482825700283 -243 8.332511595730896 6.662870524665738 1.666666666937509 -244 1.659476797213571 1.641166790796788 8.333333332432243 -245 6.651358102186887 8.333333331446589 6.665729135182278 -246 4.996765277527716 6.653675864145596 1.666666669892875 -247 6.660931836157782 1.649309184852187 8.333333333926896 -248 6.638537756742497 1.666666666445091 6.665974027330344 -249 3.331386262994759 6.656799061971727 8.333333334876567 -250 4.988687405622154 3.294823441600734 8.333333333444344 -251 6.660538599871673 3.312736539238734 8.333333333244161 -252 8.332031501668784 3.327655810609839 1.666666663483084 -253 1.660497912810677 1.651563432065564 1.666666666532634 -254 1.666666666666671 8.333333332815677 8.333333333333332 -255 1.623165746041504 1.6666666666034 4.99391383873798 -256 8.333333333273636 3.333333333333333 5.000000000000004 -257 3.332833974604121 8.330588915852831 8.33333333503149 -258 4.97579818627555 8.333333331156744 6.664430588732456 -259 1.660298963080085 3.305559657733506 8.33333333203548 -260 3.28326020484872 8.333333336724735 3.326681050109922 -261 8.331602596644176 1.661100149132812 8.333333335117558 -262 1.666666666603764 8.333333333333336 5.000000000000004 -263 6.665924415024502 8.330016659525276 1.666666667557128 -264 3.308721581528395 8.333333331465084 6.663488840038644 -265 3.286658535168235 1.666666666599032 6.66196704619375 -266 4.937771380378987 1.666666665951131 3.33355951443616 -267 8.326507895108154 8.333333331451792 6.666387147756299 -268 8.332585110528434 6.663185298000827 8.333333332584449 -269 4.982918784821877 1.6666666659068 8.332255412505511 -270 1.664258409514333 6.659472931888667 1.666666667233425 -271 3.322620058061121 3.296169448190562 8.33333333167775 -272 8.33333333683429 5.000000000000005 8.333333333333336 -273 8.333333333060155 6.666666666666669 5.000000000000004 -274 1.666666666327729 1.666666666666668 3.333333333333338 -275 4.99888149235128 8.328781672181401 1.666666669469232 -276 8.333333333095132 1.666666666666668 3.333333333333338 -277 8.333333334111209 6.666666666666671 6.666666666666671 -278 1.666666666231004 1.666666666666668 6.666666666666671 -279 8.319477494274667 8.333333334546891 3.333169972868349 -280 8.333333333713604 5.000000000000002 3.333333333333337 -281 8.333333335124941 3.333333333333335 6.66666666666667 -282 6.62749052424298 1.66666666623669 3.334574749423636 -283 1.66666666508557 3.333333333333335 5.000000000000002 -284 8.333333332051213 6.666666666666668 3.333333333333337 -285 1.637778017620791 8.333333342488848 3.327652077469392 -286 1.666666666089987 5.000000000000003 3.333333333333334 -287 6.62594906800776 1.666666666783335 4.999916920658317 -288 4.991559204424686 3.303632966319623 1.666666666035724 -289 1.666666665979094 3.333333333333337 3.333333333333335 -290 1.666666665657806 5.000000000000005 8.333333333333336 -291 1.666666666853695 5.000000000000003 6.666666666666669 -292 1.651821825244905 8.333333330438013 6.664257878478076 -293 3.270150200748045 1.66666666673903 4.995167900443951 -294 1.666666665802565 8.333333333333334 1.66666666666667 -295 3.332149742463294 8.329006867792472 1.666666666152793 -296 4.956018879142214 8.333333342722435 3.330340837991536 -297 4.961068569266943 8.333333337999449 4.996722349190394 -298 4.997593572036546 6.655043997091015 8.33333333485422 -299 3.32335131306381 3.303903840008758 1.666666665489136 -300 3.327865240662655 4.977063523995539 8.333333333931193 -301 4.993810674308104 4.973188365824303 8.333333333793638 -302 4.943783761023866 1.666666666230662 4.998321579025965 -303 4.990499273962526 8.333333328923789 8.332390464609082 -304 6.662704487145006 3.317607957790754 1.666666665319158 -305 6.64528602688278 1.666666667201381 1.670374112938904 -306 8.315733273001639 1.666666666751797 5.000155748895926 -307 6.663143942894656 4.984579210873332 1.666666668138241 -308 6.662741535306179 4.984497939484297 8.333333333366063 -309 8.333333332710026 3.333333333333334 3.333333333333337 -310 1.660403493758395 3.316764122600103 1.666666666364711 -311 6.641946928853725 8.33333332998944 4.998695732656651 -312 3.329655718808052 6.653762669731537 1.666666670056743 -313 8.333333334210733 1.666666666666668 6.666666666666671 -314 3.278505243836852 1.666666665982843 3.33055000927853 -315 8.333333333196052 5.000000000000001 5.000000000000003 -316 4.958530187990556 1.666666666701269 6.664632353924778 -317 3.333333335751559 6.666666666666669 5.000000000000002 -318 3.33333333945139 5.000000000000002 5.000000000000002 -319 3.325134218047016 4.973125396957162 3.333333331397721 -320 6.66479451352513 6.658229402591791 3.333333336627152 -321 4.992487513742513 4.968476001600282 3.333333328810805 -322 3.333333335727594 6.66666666666667 6.666666666666671 -323 6.66666666911012 5.000000000000003 6.666666666666671 -324 3.333333333488395 3.333333333333335 6.666666666666669 -325 6.66508040097065 6.65903418390889 6.666666667932162 -326 3.333333334384611 6.666666666666668 3.333333333333335 -327 4.996765277527716 6.653675864145596 3.333333337116327 -328 6.638537756742497 3.333333332536744 6.665974027330344 -329 4.988687405622154 3.294823441600734 6.666666666679443 -330 6.666666667093498 3.333333333333333 5.000000000000004 -331 4.97579818627555 6.666666665298601 6.664430588732456 -332 4.937771380378987 3.333333331698917 3.33355951443616 -333 6.666666667487235 6.666666666666669 5.000000000000004 -334 6.666666667036193 5.000000000000002 3.333333333333337 -335 6.62749052424298 3.333333332299917 3.334574749423636 -336 3.33333333143754 3.333333333333335 5.000000000000002 -337 3.333333331313443 3.333333333333337 3.333333333333335 -338 3.333333334724964 5.000000000000003 6.666666666666669 -339 4.961068569266943 6.666666676876886 4.996722349190394 -340 4.993810674308104 4.973188365824303 6.666666667605718 -341 4.943783761023866 3.33333333244255 4.998321579025965 -342 6.666666666438698 5.000000000000001 5.000000000000003 -343 5.000000007071594 5.000000000000002 5.000000000000002 -$EndNodes -$Elements -512 -1 15 2 0 1 1 -2 15 2 0 2 2 -3 15 2 0 3 3 -4 15 2 0 4 4 -5 15 2 0 5 5 -6 15 2 0 6 6 -7 15 2 0 7 7 -8 15 2 0 8 8 -9 1 2 0 1 1 9 -10 1 2 0 1 9 10 -11 1 2 0 1 10 11 -12 1 2 0 1 11 12 -13 1 2 0 1 12 13 -14 1 2 0 1 13 2 -15 1 2 0 2 1 14 -16 1 2 0 2 14 15 -17 1 2 0 2 15 16 -18 1 2 0 2 16 17 -19 1 2 0 2 17 18 -20 1 2 0 2 18 3 -21 1 2 0 3 3 19 -22 1 2 0 3 19 20 -23 1 2 0 3 20 21 -24 1 2 0 3 21 22 -25 1 2 0 3 22 23 -26 1 2 0 3 23 4 -27 1 2 0 4 2 24 -28 1 2 0 4 24 25 -29 1 2 0 4 25 26 -30 1 2 0 4 26 27 -31 1 2 0 4 27 28 -32 1 2 0 4 28 4 -33 1 2 0 5 5 29 -34 1 2 0 5 29 30 -35 1 2 0 5 30 31 -36 1 2 0 5 31 32 -37 1 2 0 5 32 33 -38 1 2 0 5 33 6 -39 1 2 0 6 5 34 -40 1 2 0 6 34 35 -41 1 2 0 6 35 36 -42 1 2 0 6 36 37 -43 1 2 0 6 37 38 -44 1 2 0 6 38 7 -45 1 2 0 7 7 39 -46 1 2 0 7 39 40 -47 1 2 0 7 40 41 -48 1 2 0 7 41 42 -49 1 2 0 7 42 43 -50 1 2 0 7 43 8 -51 1 2 0 8 6 44 -52 1 2 0 8 44 45 -53 1 2 0 8 45 46 -54 1 2 0 8 46 47 -55 1 2 0 8 47 48 -56 1 2 0 8 48 8 -57 1 2 0 9 1 49 -58 1 2 0 9 49 50 -59 1 2 0 9 50 51 -60 1 2 0 9 51 52 -61 1 2 0 9 52 53 -62 1 2 0 9 53 5 -63 1 2 0 10 2 54 -64 1 2 0 10 54 55 -65 1 2 0 10 55 56 -66 1 2 0 10 56 57 -67 1 2 0 10 57 58 -68 1 2 0 10 58 6 -69 1 2 0 11 3 59 -70 1 2 0 11 59 60 -71 1 2 0 11 60 61 -72 1 2 0 11 61 62 -73 1 2 0 11 62 63 -74 1 2 0 11 63 7 -75 1 2 0 12 4 64 -76 1 2 0 12 64 65 -77 1 2 0 12 65 66 -78 1 2 0 12 66 67 -79 1 2 0 12 67 68 -80 1 2 0 12 68 8 -81 3 2 0 1 82 71 88 93 -82 3 2 0 1 82 93 85 91 -83 3 2 0 1 82 84 92 91 -84 3 2 0 1 91 86 83 92 -85 3 2 0 1 86 91 85 90 -86 3 2 0 1 4 28 75 23 -87 3 2 0 1 75 28 27 88 -88 3 2 0 1 18 81 77 17 -89 3 2 0 1 77 90 86 72 -90 3 2 0 1 75 23 22 89 -91 3 2 0 1 14 15 87 79 -92 3 2 0 1 26 27 88 71 -93 3 2 0 1 17 77 72 16 -94 3 2 0 1 89 75 88 93 -95 3 2 0 1 15 87 72 16 -96 3 2 0 1 87 83 86 72 -97 3 2 0 1 89 70 85 93 -98 3 2 0 1 87 83 69 79 -99 3 2 0 1 3 18 81 19 -100 3 2 0 1 26 25 74 71 -101 3 2 0 1 25 74 80 24 -102 3 2 0 1 70 21 22 89 -103 3 2 0 1 80 13 2 24 -104 3 2 0 1 12 78 80 13 -105 3 2 0 1 19 20 73 81 -106 3 2 0 1 82 71 74 84 -107 3 2 0 1 85 70 73 90 -108 3 2 0 1 20 73 70 21 -109 3 2 0 1 74 80 78 84 -110 3 2 0 1 14 79 9 1 -111 3 2 0 1 73 90 77 81 -112 3 2 0 1 10 9 79 69 -113 3 2 0 1 78 12 11 76 -114 3 2 0 1 78 84 92 76 -115 3 2 0 1 69 10 11 76 -116 3 2 0 1 69 76 92 83 -117 3 2 0 2 114 116 107 118 -118 3 2 0 2 108 110 117 116 -119 3 2 0 2 115 114 116 108 -120 3 2 0 2 96 100 115 114 -121 3 2 0 2 117 109 107 116 -122 3 2 0 2 48 98 112 47 -123 3 2 0 2 8 43 98 48 -124 3 2 0 2 114 118 112 96 -125 3 2 0 2 38 37 102 106 -126 3 2 0 2 111 35 34 104 -127 3 2 0 2 35 111 97 36 -128 3 2 0 2 37 102 97 36 -129 3 2 0 2 110 94 104 111 -130 3 2 0 2 98 43 42 113 -131 3 2 0 2 110 111 97 117 -132 3 2 0 2 97 117 109 102 -133 3 2 0 2 100 105 103 115 -134 3 2 0 2 113 95 107 118 -135 3 2 0 2 6 33 105 44 -136 3 2 0 2 118 113 98 112 -137 3 2 0 2 96 112 47 46 -138 3 2 0 2 45 44 105 100 -139 3 2 0 2 32 33 105 103 -140 3 2 0 2 45 46 96 100 -141 3 2 0 2 7 39 106 38 -142 3 2 0 2 41 95 113 42 -143 3 2 0 2 34 5 29 104 -144 3 2 0 2 106 39 40 99 -145 3 2 0 2 95 99 40 41 -146 3 2 0 2 94 30 29 104 -147 3 2 0 2 102 109 99 106 -148 3 2 0 2 107 109 99 95 -149 3 2 0 2 31 32 103 101 -150 3 2 0 2 94 30 31 101 -151 3 2 0 2 115 108 101 103 -152 3 2 0 2 110 94 101 108 -153 3 2 0 3 136 141 137 138 -154 3 2 0 3 143 135 137 141 -155 3 2 0 3 6 58 128 33 -156 3 2 0 3 127 128 33 32 -157 3 2 0 3 53 5 29 125 -158 3 2 0 3 141 140 131 136 -159 3 2 0 3 30 29 125 130 -160 3 2 0 3 57 58 128 124 -161 3 2 0 3 143 141 140 139 -162 3 2 0 3 31 121 127 32 -163 3 2 0 3 136 121 127 131 -164 3 2 0 3 130 138 134 125 -165 3 2 0 3 125 134 52 53 -166 3 2 0 3 30 31 121 130 -167 3 2 0 3 137 138 134 119 -168 3 2 0 3 136 138 130 121 -169 3 2 0 3 133 56 55 120 -170 3 2 0 3 124 133 140 131 -171 3 2 0 3 126 13 2 54 -172 3 2 0 3 120 133 140 139 -173 3 2 0 3 124 133 56 57 -174 3 2 0 3 127 131 124 128 -175 3 2 0 3 12 13 126 132 -176 3 2 0 3 126 54 55 120 -177 3 2 0 3 52 51 119 134 -178 3 2 0 3 139 132 126 120 -179 3 2 0 3 1 49 123 9 -180 3 2 0 3 123 49 50 142 -181 3 2 0 3 143 122 132 139 -182 3 2 0 3 10 9 123 129 -183 3 2 0 3 132 12 11 122 -184 3 2 0 3 142 119 51 50 -185 3 2 0 3 123 142 135 129 -186 3 2 0 3 129 122 11 10 -187 3 2 0 3 129 135 143 122 -188 3 2 0 3 119 137 135 142 -189 3 2 0 4 4 64 149 23 -190 3 2 0 4 68 152 43 8 -191 3 2 0 4 152 43 42 150 -192 3 2 0 4 162 167 166 159 -193 3 2 0 4 67 156 152 68 -194 3 2 0 4 166 167 163 160 -195 3 2 0 4 67 66 146 156 -196 3 2 0 4 154 65 64 149 -197 3 2 0 4 154 65 66 146 -198 3 2 0 4 7 63 151 39 -199 3 2 0 4 157 150 165 163 -200 3 2 0 4 157 41 42 150 -201 3 2 0 4 23 22 158 149 -202 3 2 0 4 156 165 150 152 -203 3 2 0 4 39 151 144 40 -204 3 2 0 4 165 163 167 161 -205 3 2 0 4 146 156 165 161 -206 3 2 0 4 157 144 40 41 -207 3 2 0 4 160 144 157 163 -208 3 2 0 4 59 148 19 3 -209 3 2 0 4 146 161 164 154 -210 3 2 0 4 158 149 154 164 -211 3 2 0 4 161 164 162 167 -212 3 2 0 4 145 158 22 21 -213 3 2 0 4 151 63 62 155 -214 3 2 0 4 162 145 158 164 -215 3 2 0 4 145 162 159 153 -216 3 2 0 4 62 61 147 155 -217 3 2 0 4 19 148 153 20 -218 3 2 0 4 61 147 168 60 -219 3 2 0 4 144 151 155 160 -220 3 2 0 4 60 168 148 59 -221 3 2 0 4 153 145 21 20 -222 3 2 0 4 155 160 166 147 -223 3 2 0 4 147 166 159 168 -224 3 2 0 4 168 148 153 159 -225 3 2 0 5 189 180 171 186 -226 3 2 0 5 187 185 192 188 -227 3 2 0 5 189 184 175 180 -228 3 2 0 5 189 188 192 186 -229 3 2 0 5 38 37 181 174 -230 3 2 0 5 169 35 34 176 -231 3 2 0 5 169 177 36 35 -232 3 2 0 5 183 191 181 174 -233 3 2 0 5 37 181 177 36 -234 3 2 0 5 183 191 186 171 -235 3 2 0 5 18 17 184 175 -236 3 2 0 5 169 178 190 177 -237 3 2 0 5 177 190 191 181 -238 3 2 0 5 7 63 174 38 -239 3 2 0 5 187 193 173 179 -240 3 2 0 5 176 53 5 34 -241 3 2 0 5 14 15 193 173 -242 3 2 0 5 185 178 190 192 -243 3 2 0 5 17 184 170 16 -244 3 2 0 5 193 15 16 170 -245 3 2 0 5 192 190 191 186 -246 3 2 0 5 53 176 182 52 -247 3 2 0 5 59 175 18 3 -248 3 2 0 5 193 170 188 187 -249 3 2 0 5 182 178 169 176 -250 3 2 0 5 174 63 62 183 -251 3 2 0 5 170 188 189 184 -252 3 2 0 5 62 61 171 183 -253 3 2 0 5 61 171 180 60 -254 3 2 0 5 52 182 172 51 -255 3 2 0 5 60 180 175 59 -256 3 2 0 5 14 1 49 173 -257 3 2 0 5 49 50 179 173 -258 3 2 0 5 182 172 185 178 -259 3 2 0 5 50 51 172 179 -260 3 2 0 5 187 185 172 179 -261 3 2 0 6 210 216 217 214 -262 3 2 0 6 4 28 201 64 -263 3 2 0 6 216 217 212 209 -264 3 2 0 6 48 47 205 202 -265 3 2 0 6 68 8 48 202 -266 3 2 0 6 215 200 202 205 -267 3 2 0 6 201 28 27 206 -268 3 2 0 6 200 67 68 202 -269 3 2 0 6 6 58 199 44 -270 3 2 0 6 46 196 205 47 -271 3 2 0 6 207 194 210 214 -272 3 2 0 6 67 200 207 66 -273 3 2 0 6 45 204 199 44 -274 3 2 0 6 210 194 201 206 -275 3 2 0 6 65 64 201 194 -276 3 2 0 6 204 196 46 45 -277 3 2 0 6 65 194 207 66 -278 3 2 0 6 215 200 207 214 -279 3 2 0 6 214 217 211 215 -280 3 2 0 6 197 206 27 26 -281 3 2 0 6 205 215 211 196 -282 3 2 0 6 208 199 58 57 -283 3 2 0 6 209 216 197 218 -284 3 2 0 6 213 211 196 204 -285 3 2 0 6 213 212 217 211 -286 3 2 0 6 213 208 199 204 -287 3 2 0 6 197 206 210 216 -288 3 2 0 6 26 25 218 197 -289 3 2 0 6 25 218 198 24 -290 3 2 0 6 212 209 203 195 -291 3 2 0 6 195 208 213 212 -292 3 2 0 6 56 195 203 55 -293 3 2 0 6 198 24 2 54 -294 3 2 0 6 208 195 56 57 -295 3 2 0 6 54 198 203 55 -296 3 2 0 6 209 218 198 203 -297 5 2 0 1 114 118 107 116 232 277 273 315 -298 5 2 0 1 108 116 117 110 256 315 280 309 -299 5 2 0 1 281 232 114 115 256 315 116 108 -300 5 2 0 1 82 291 230 93 71 290 239 88 -301 5 2 0 1 291 230 222 226 82 93 85 91 -302 5 2 0 1 318 338 291 226 336 324 236 283 -303 5 2 0 1 82 84 92 91 291 236 283 226 -304 5 2 0 1 226 283 92 91 286 289 83 86 -305 5 2 0 1 96 100 223 272 114 115 281 232 -306 5 2 0 1 286 226 222 241 86 91 85 90 -307 5 2 0 1 280 117 109 284 315 116 107 273 -308 5 2 0 1 318 226 291 338 317 222 230 322 -309 5 2 0 1 249 298 214 210 300 301 217 216 -310 5 2 0 1 4 23 149 64 28 75 254 201 -311 5 2 0 1 312 189 180 295 246 186 171 275 -312 5 2 0 1 187 299 288 185 188 227 229 192 -313 5 2 0 1 216 209 271 300 217 212 250 301 -314 5 2 0 1 312 189 184 270 295 180 175 294 -315 5 2 0 1 48 202 225 98 47 205 268 112 -316 5 2 0 1 68 202 225 152 8 48 98 43 -317 5 2 0 1 189 312 227 188 186 246 229 192 -318 5 2 0 1 114 118 277 232 96 112 268 272 -319 5 2 0 1 38 174 181 37 106 231 243 102 -320 5 2 0 1 215 200 202 205 237 233 225 268 -321 5 2 0 1 252 169 35 111 220 176 34 104 -322 5 2 0 1 169 177 240 252 35 36 97 111 -323 5 2 0 1 183 263 228 191 174 231 243 181 -324 5 2 0 1 37 102 243 181 36 97 240 177 -325 5 2 0 1 254 75 88 239 201 28 27 206 -326 5 2 0 1 183 171 186 191 263 275 246 228 -327 5 2 0 1 18 81 77 17 175 294 270 184 -328 5 2 0 1 309 110 94 276 252 111 104 220 -329 5 2 0 1 225 152 43 98 267 150 42 113 -330 5 2 0 1 287 282 138 136 302 266 137 141 -331 5 2 0 1 110 111 252 309 117 97 240 280 -332 5 2 0 1 97 240 280 117 102 243 284 109 -333 5 2 0 1 162 167 297 242 159 166 296 260 -334 5 2 0 1 143 293 302 141 135 314 266 137 -335 5 2 0 1 200 202 225 233 67 68 152 156 -336 5 2 0 1 100 105 261 223 115 103 313 281 -337 5 2 0 1 267 277 273 221 113 118 107 95 -338 5 2 0 1 6 44 105 33 58 199 261 128 -339 5 2 0 1 118 112 98 113 277 268 225 267 -340 5 2 0 1 315 232 277 273 342 323 325 333 -341 5 2 0 1 96 46 47 112 272 196 205 268 -342 5 2 0 1 166 167 297 296 160 163 311 238 -343 5 2 0 1 303 298 214 207 257 249 210 194 -344 5 2 0 1 252 169 178 304 240 177 190 307 -345 5 2 0 1 67 200 207 66 156 233 303 146 -346 5 2 0 1 45 100 105 44 204 223 261 199 -347 5 2 0 1 315 273 284 280 342 333 320 334 -348 5 2 0 1 249 239 206 210 257 254 201 194 -349 5 2 0 1 154 257 194 65 149 254 201 64 -350 5 2 0 1 240 177 190 307 243 181 191 228 -351 5 2 0 1 127 128 33 32 313 261 105 103 -352 5 2 0 1 204 196 46 45 223 272 96 100 -353 5 2 0 1 154 146 66 65 257 303 207 194 -354 5 2 0 1 77 72 224 270 90 86 286 241 -355 5 2 0 1 7 39 151 63 38 106 231 174 -356 5 2 0 1 299 310 193 187 219 253 173 179 -357 5 2 0 1 157 221 267 150 163 311 245 165 -358 5 2 0 1 157 221 95 41 150 267 113 42 -359 5 2 0 1 237 298 214 215 233 303 207 200 -360 5 2 0 1 326 327 321 319 312 246 229 227 -361 5 2 0 1 176 220 125 53 34 104 29 5 -362 5 2 0 1 75 254 292 89 23 149 158 22 -363 5 2 0 1 298 237 215 214 301 308 211 217 -364 5 2 0 1 233 245 165 156 225 267 150 152 -365 5 2 0 1 141 140 131 136 302 316 248 287 -366 5 2 0 1 106 39 151 231 99 40 144 279 -367 5 2 0 1 165 161 167 163 245 258 297 311 -368 5 2 0 1 14 15 87 79 173 193 310 253 -369 5 2 0 1 288 304 178 185 229 307 190 192 -370 5 2 0 1 197 290 71 26 206 239 88 27 -371 5 2 0 1 268 272 196 205 237 308 211 215 -372 5 2 0 1 146 156 233 303 161 165 245 258 -373 5 2 0 1 221 95 99 279 157 41 40 144 -374 5 2 0 1 94 30 130 276 104 29 125 220 -375 5 2 0 1 208 57 124 247 199 58 128 261 -376 5 2 0 1 17 16 170 184 77 72 224 270 -377 5 2 0 1 160 238 279 144 163 311 221 157 -378 5 2 0 1 280 284 243 240 334 320 228 307 -379 5 2 0 1 243 284 109 102 231 279 99 106 -380 5 2 0 1 89 93 230 292 75 88 239 254 -381 5 2 0 1 293 302 316 265 143 141 140 139 -382 5 2 0 1 273 107 109 284 221 95 99 279 -383 5 2 0 1 193 170 16 15 310 224 72 87 -384 5 2 0 1 31 121 127 32 101 306 313 103 -385 5 2 0 1 192 229 307 190 186 246 228 191 -386 5 2 0 1 310 224 72 87 289 286 86 83 -387 5 2 0 1 287 248 131 136 306 313 127 121 -388 5 2 0 1 277 268 225 267 325 237 233 245 -389 5 2 0 1 130 276 220 125 138 282 305 134 -390 5 2 0 1 271 259 218 209 300 290 197 216 -391 5 2 0 1 251 213 204 223 308 211 196 272 -392 5 2 0 1 292 262 222 230 89 70 85 93 -393 5 2 0 1 213 211 308 251 212 217 301 250 -394 5 2 0 1 213 251 223 204 208 247 261 199 -395 5 2 0 1 268 272 308 237 277 232 323 325 -396 5 2 0 1 125 220 176 53 134 305 182 52 -397 5 2 0 1 94 276 130 30 101 306 121 31 -398 5 2 0 1 87 83 69 79 310 289 274 253 -399 5 2 0 1 137 119 235 266 138 134 305 282 -400 5 2 0 1 59 175 294 148 3 18 81 19 -401 5 2 0 1 281 256 108 115 313 306 101 103 -402 5 2 0 1 309 256 306 276 110 108 101 94 -403 5 2 0 1 303 146 161 258 257 154 164 264 -404 5 2 0 1 287 282 276 306 136 138 130 121 -405 5 2 0 1 193 170 224 310 187 188 227 299 -406 5 2 0 1 292 158 149 254 264 164 154 257 -407 5 2 0 1 290 197 216 300 239 206 210 249 -408 5 2 0 1 258 264 164 161 297 242 162 167 -409 5 2 0 1 277 273 333 325 267 221 311 245 -410 5 2 0 1 297 339 333 311 258 331 325 245 -411 5 2 0 1 258 245 233 303 331 325 237 298 -412 5 2 0 1 26 71 290 197 25 74 259 218 -413 5 2 0 1 25 218 198 24 74 259 244 80 -414 5 2 0 1 250 269 195 212 271 234 203 209 -415 5 2 0 1 262 70 89 292 145 21 22 158 -416 5 2 0 1 182 176 169 178 305 220 252 304 -417 5 2 0 1 270 224 227 312 241 286 319 326 -418 5 2 0 1 195 269 250 212 208 247 251 213 -419 5 2 0 1 231 151 63 174 263 155 62 183 -420 5 2 0 1 133 269 234 120 56 195 203 55 -421 5 2 0 1 300 290 239 249 338 291 230 322 -422 5 2 0 1 247 269 316 248 124 133 140 131 -423 5 2 0 1 221 279 284 273 311 238 320 333 -424 5 2 0 1 242 162 145 262 264 164 158 292 -425 5 2 0 1 244 126 54 198 80 13 2 24 -426 5 2 0 1 224 227 188 170 270 312 189 184 -427 5 2 0 1 145 162 242 262 153 159 260 285 -428 5 2 0 1 234 120 139 265 269 133 140 316 -429 5 2 0 1 247 124 57 208 269 133 56 195 -430 5 2 0 1 127 128 261 313 131 124 247 248 -431 5 2 0 1 12 13 126 132 78 80 244 278 -432 5 2 0 1 317 326 319 318 222 241 286 226 -433 5 2 0 1 62 155 147 61 183 263 275 171 -434 5 2 0 1 239 249 322 230 254 257 264 292 -435 5 2 0 1 19 81 294 148 20 73 285 153 -436 5 2 0 1 286 224 227 319 289 310 299 337 -437 5 2 0 1 264 258 303 257 322 331 298 249 -438 5 2 0 1 334 320 228 307 321 327 246 229 -439 5 2 0 1 291 82 84 236 290 71 74 259 -440 5 2 0 1 223 281 313 261 251 328 248 247 -441 5 2 0 1 61 171 275 147 60 180 295 168 -442 5 2 0 1 144 160 155 151 279 238 263 231 -443 5 2 0 1 126 54 55 120 244 198 203 234 -444 5 2 0 1 52 134 305 182 51 119 235 172 -445 5 2 0 1 85 222 241 90 70 262 285 73 -446 5 2 0 1 323 325 331 340 308 237 298 301 -447 5 2 0 1 60 180 295 168 59 175 294 148 -448 5 2 0 1 153 285 262 145 20 73 70 21 -449 5 2 0 1 271 300 338 324 259 290 291 236 -450 5 2 0 1 139 120 234 265 132 126 244 278 -451 5 2 0 1 74 80 78 84 259 244 278 236 -452 5 2 0 1 342 330 256 315 323 328 281 232 -453 5 2 0 1 271 259 244 234 209 218 198 203 -454 5 2 0 1 308 323 232 272 251 328 281 223 -455 5 2 0 1 14 79 9 1 173 253 123 49 -456 5 2 0 1 123 253 219 142 49 173 179 50 -457 5 2 0 1 238 263 231 279 320 228 243 284 -458 5 2 0 1 322 249 300 338 331 298 301 340 -459 5 2 0 1 182 172 235 305 178 185 288 304 -460 5 2 0 1 337 336 318 319 289 283 226 286 -461 5 2 0 1 143 122 132 139 293 255 278 265 -462 5 2 0 1 73 81 294 285 90 77 270 241 -463 5 2 0 1 10 69 274 129 9 79 253 123 -464 5 2 0 1 304 335 282 305 252 309 276 220 -465 5 2 0 1 263 238 160 155 275 296 166 147 -466 5 2 0 1 287 306 256 330 248 313 281 328 -467 5 2 0 1 342 333 339 343 323 325 331 340 -468 5 2 0 1 278 78 76 255 132 12 11 122 -469 5 2 0 1 317 326 241 222 242 260 285 262 -470 5 2 0 1 339 333 320 327 297 311 238 296 -471 5 2 0 1 317 242 262 222 322 264 292 230 -472 5 2 0 1 321 327 339 343 334 320 333 342 -473 5 2 0 1 322 331 340 338 317 339 343 318 -474 5 2 0 1 142 119 51 50 219 235 172 179 -475 5 2 0 1 280 315 342 334 309 256 330 335 -476 5 2 0 1 329 340 338 324 250 301 300 271 -477 5 2 0 1 253 274 129 123 219 314 135 142 -478 5 2 0 1 335 334 280 309 304 307 240 252 -479 5 2 0 1 275 295 260 296 147 168 159 166 -480 5 2 0 1 339 297 242 317 331 258 264 322 -481 5 2 0 1 317 318 319 326 339 343 321 327 -482 5 2 0 1 287 282 335 330 306 276 309 256 -483 5 2 0 1 246 228 263 275 327 320 238 296 -484 5 2 0 1 278 255 76 78 236 283 92 84 -485 5 2 0 1 295 168 148 294 260 159 153 285 -486 5 2 0 1 329 324 265 316 250 271 234 269 -487 5 2 0 1 330 287 302 341 335 282 266 332 -488 5 2 0 1 244 234 265 278 259 271 324 236 -489 5 2 0 1 274 255 76 69 129 122 11 10 -490 5 2 0 1 274 255 122 129 314 293 143 135 -491 5 2 0 1 341 336 337 332 302 293 314 266 -492 5 2 0 1 293 336 337 314 255 283 289 274 -493 5 2 0 1 327 296 260 326 339 297 242 317 -494 5 2 0 1 187 185 288 299 179 172 235 219 -495 5 2 0 1 337 336 341 332 319 318 343 321 -496 5 2 0 1 336 283 255 293 324 236 278 265 -497 5 2 0 1 293 336 324 265 302 341 329 316 -498 5 2 0 1 274 255 283 289 69 76 92 83 -499 5 2 0 1 260 295 275 296 326 312 246 327 -500 5 2 0 1 316 329 250 269 248 328 251 247 -501 5 2 0 1 270 294 285 241 312 295 260 326 -502 5 2 0 1 329 340 301 250 328 323 308 251 -503 5 2 0 1 332 288 304 335 266 235 305 282 -504 5 2 0 1 235 219 142 119 266 314 135 137 -505 5 2 0 1 337 299 288 332 314 219 235 266 -506 5 2 0 1 310 299 337 289 253 219 314 274 -507 5 2 0 1 343 341 329 340 318 336 324 338 -508 5 2 0 1 329 328 330 341 340 323 342 343 -509 5 2 0 1 337 332 288 299 319 321 229 227 -510 5 2 0 1 304 288 332 335 307 229 321 334 -511 5 2 0 1 287 330 341 302 248 328 329 316 -512 5 2 0 1 334 342 343 321 335 330 341 332 -$EndElements diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 975ed6263..e853c5c62 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -490,7 +490,7 @@ def f(x): print(eoc_rec) assert ( eoc_rec.order_estimate() >= order-order_slack - or eoc_rec.max_error() < 1e-13) + or eoc_rec.max_error() < 1.5e-13) # }}} From 02a3a4824f59aea16a66c66774e3b9c8d70d3e89 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 5 Nov 2020 19:20:10 -0600 Subject: [PATCH 11/51] use mesh order for resampling matrix --- meshmode/discretization/poly_element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index ee9949798..48e09778a 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -393,7 +393,7 @@ def from_mesh_interp_matrix(self): from modepy.modes import legendre_tensor_product_basis meg = self.mesh_el_group return mp.resampling_matrix( - legendre_tensor_product_basis(self.dim, self.order), + legendre_tensor_product_basis(self.dim, meg.order), self.unit_nodes, meg.unit_nodes) From 673dd009b87e6da24016bdcb78b8b2e90bfe27c3 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 5 Nov 2020 20:48:23 -0600 Subject: [PATCH 12/51] fix quad weights for tensor products --- meshmode/discretization/poly_element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 48e09778a..ac54fb833 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -422,7 +422,7 @@ def weights(self): return np.fromiter( (np.prod(w) for w in product(weights, repeat=self.dim)), dtype=np.float, - count=self.order**self.dim) + count=(self.order + 1)**self.dim) class LegendreGaussLobattoTensorProductElementGroup( From a82685f5c043e19125bd9bef9e62beff41286265 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 6 Nov 2020 09:06:24 -0600 Subject: [PATCH 13/51] add back removed mesh (huh) --- test/cubed-cube.msh | 864 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 864 insertions(+) create mode 100644 test/cubed-cube.msh diff --git a/test/cubed-cube.msh b/test/cubed-cube.msh new file mode 100644 index 000000000..dde325a0a --- /dev/null +++ b/test/cubed-cube.msh @@ -0,0 +1,864 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +343 +1 0 0 0 +2 0 0 10 +3 0 10 0 +4 0 10 10 +5 10 0 0 +6 10 0 10 +7 10 10 0 +8 10 10 10 +9 0 0 1.666666666666672 +10 0 0 3.333333333333344 +11 0 0 5.000000000000009 +12 0 0 6.666666666666672 +13 0 0 8.333333333333336 +14 0 1.666666666666672 0 +15 0 3.333333333333344 0 +16 0 5.000000000000009 0 +17 0 6.666666666666672 0 +18 0 8.333333333333336 0 +19 0 10 1.666666666666672 +20 0 10 3.333333333333344 +21 0 10 5.000000000000009 +22 0 10 6.666666666666672 +23 0 10 8.333333333333336 +24 0 1.666666666666672 10 +25 0 3.333333333333344 10 +26 0 5.000000000000009 10 +27 0 6.666666666666672 10 +28 0 8.333333333333336 10 +29 10 0 1.666666666666672 +30 10 0 3.333333333333344 +31 10 0 5.000000000000009 +32 10 0 6.666666666666672 +33 10 0 8.333333333333336 +34 10 1.666666666666672 0 +35 10 3.333333333333344 0 +36 10 5.000000000000009 0 +37 10 6.666666666666672 0 +38 10 8.333333333333336 0 +39 10 10 1.666666666666672 +40 10 10 3.333333333333344 +41 10 10 5.000000000000009 +42 10 10 6.666666666666672 +43 10 10 8.333333333333336 +44 10 1.666666666666672 10 +45 10 3.333333333333344 10 +46 10 5.000000000000009 10 +47 10 6.666666666666672 10 +48 10 8.333333333333336 10 +49 1.666666666666672 0 0 +50 3.333333333333344 0 0 +51 5.000000000000009 0 0 +52 6.666666666666672 0 0 +53 8.333333333333336 0 0 +54 1.666666666666672 0 10 +55 3.333333333333344 0 10 +56 5.000000000000009 0 10 +57 6.666666666666672 0 10 +58 8.333333333333336 0 10 +59 1.666666666666672 10 0 +60 3.333333333333344 10 0 +61 5.000000000000009 10 0 +62 6.666666666666672 10 0 +63 8.333333333333336 10 0 +64 1.666666666666672 10 10 +65 3.333333333333344 10 10 +66 5.000000000000009 10 10 +67 6.666666666666672 10 10 +68 8.333333333333336 10 10 +69 0 1.666666666666668 3.333333333333338 +70 0 8.333333333333336 5.000000000000004 +71 0 5.000000000000005 8.333333333333336 +72 0 5.000000000000006 1.666666666666667 +73 0 8.333333333333334 3.333333333333338 +74 0 3.333333333333338 8.333333333333334 +75 0 8.333333333333336 8.333333333333336 +76 0 1.666666666666668 5.000000000000005 +77 0 6.666666666666671 1.666666666666668 +78 0 1.666666666666668 6.666666666666671 +79 0 1.66666666666667 1.666666666666669 +80 0 1.66666666666667 8.333333333333334 +81 0 8.333333333333334 1.66666666666667 +82 0 5.000000000000003 6.666666666666669 +83 0 3.333333333333337 3.333333333333335 +84 0 3.333333333333335 6.666666666666669 +85 0 6.666666666666669 5.000000000000002 +86 0 5.000000000000003 3.333333333333334 +87 0 3.333333333333339 1.666666666666668 +88 0 6.666666666666671 8.333333333333336 +89 0 8.333333333333336 6.666666666666671 +90 0 6.666666666666668 3.333333333333335 +91 0 5.000000000000002 5.000000000000002 +92 0 3.333333333333335 5.000000000000002 +93 0 6.66666666666667 6.666666666666671 +94 10 1.666666666666668 3.333333333333338 +95 10 8.333333333333336 5.000000000000004 +96 10 5.000000000000005 8.333333333333336 +97 10 5.000000000000004 1.666666666666668 +98 10 8.333333333333337 8.333333333333336 +99 10 8.333333333333334 3.333333333333338 +100 10 3.333333333333338 8.333333333333334 +101 10 1.666666666666667 5.000000000000006 +102 10 6.666666666666671 1.666666666666668 +103 10 1.666666666666668 6.666666666666671 +104 10 1.666666666666669 1.66666666666667 +105 10 1.666666666666669 8.333333333333334 +106 10 8.333333333333334 1.666666666666669 +107 10 6.666666666666669 5.000000000000004 +108 10 3.333333333333333 5.000000000000004 +109 10 6.666666666666668 3.333333333333337 +110 10 3.333333333333334 3.333333333333337 +111 10 3.333333333333337 1.666666666666668 +112 10 6.666666666666671 8.333333333333336 +113 10 8.333333333333337 6.666666666666671 +114 10 5.000000000000003 6.666666666666671 +115 10 3.333333333333335 6.66666666666667 +116 10 5.000000000000001 5.000000000000003 +117 10 5.000000000000002 3.333333333333337 +118 10 6.666666666666671 6.666666666666671 +119 4.963214356275155 0 1.668784957944797 +120 3.317582237776889 0 8.331760196778665 +121 8.315733273001639 0 5.000155748895926 +122 1.623165746041504 0 4.99391383873798 +123 1.636937952053827 0 1.662777951906097 +124 6.654194529315052 0 8.332831488754294 +125 8.323518108830998 0 1.667524092257246 +126 1.656504071342595 0 8.331974911699989 +127 8.323288233378031 0 6.666461769940527 +128 8.327361873953064 0 8.333156029608485 +129 1.625841868965969 0 3.328919181967938 +130 8.316557036413666 0 3.334049509706089 +131 6.638537756742497 0 6.665974027330344 +132 1.642552721132448 0 6.663695878577564 +133 4.982918784821877 0 8.332255412505511 +134 6.64528602688278 0 1.670374112938904 +135 3.278505243836852 0 3.33055000927853 +136 6.62594906800776 0 4.999916920658317 +137 4.937771380378987 0 3.33355951443616 +138 6.62749052424298 0 3.334574749423636 +139 3.286658535168235 0 6.66196704619375 +140 4.958530187990556 0 6.664632353924778 +141 4.943783761023866 0 4.998321579025965 +142 3.296459439269126 0 1.665682739796669 +143 3.270150200748045 0 4.995167900443951 +144 8.319477494274667 10 3.333169972868349 +145 1.642401289952334 10 4.995994122053824 +146 4.990499273962526 10 8.332390464609082 +147 4.974718901986735 10 1.665089100516185 +148 1.650806613158973 10 1.663879441264121 +149 1.666666666666671 10 8.333333333333332 +150 8.326507895108154 10 6.666387147756299 +151 8.324217698423221 10 1.666577446189837 +152 8.330421346341454 10 8.33317246986971 +153 1.637778017620791 10 3.327652077469392 +154 3.324767314036656 10 8.332279894336549 +155 6.64477579136364 10 1.666338457530614 +156 6.660165660811743 10 8.332843092640404 +157 8.321592194848483 10 4.999653680812338 +158 1.651821825244905 10 6.664257878478076 +159 3.28326020484872 10 3.326681050109922 +160 6.634104856166324 10 3.332502948194829 +161 4.97579818627555 10 6.664430588732456 +162 3.284967251501693 10 4.993482825700283 +163 6.641946928853725 10 4.998695732656651 +164 3.308721581528395 10 6.663488840038644 +165 6.651358102186887 10 6.665729135182278 +166 4.956018879142214 10 3.330340837991536 +167 4.961068569266943 10 4.996722349190394 +168 3.306257787745393 10 1.663550909007059 +169 8.332031501668784 3.327655810609839 0 +170 1.661583967397951 4.985925016393396 0 +171 4.99888149235128 8.328781672181401 0 +172 4.993595831504386 1.645761942295224 0 +173 1.660497912810677 1.651563432065564 0 +174 8.332958056923808 8.3315936822691 0 +175 1.665868128287409 8.330677281028967 0 +176 8.333333333333337 1.666666666666669 0 +177 8.332074246359143 4.994289507172165 0 +178 6.662704487145006 3.317607957790754 0 +179 3.323175186310005 1.639797188626136 0 +180 3.332149742463294 8.329006867792472 0 +181 8.332511595730896 6.662870524665738 0 +182 6.664258683158848 1.657777476463054 0 +183 6.665924415024502 8.330016659525276 0 +184 1.664258409514333 6.659472931888667 0 +185 4.991559204424686 3.303632966319623 0 +186 4.996765277527716 6.653675864145596 0 +187 3.32335131306381 3.303903840008758 0 +188 3.325134218047016 4.973125396957162 0 +189 3.329655718808052 6.653762669731537 0 +190 6.663143942894656 4.984579210873332 0 +191 6.66479451352513 6.658229402591791 0 +192 4.992487513742513 4.968476001600282 0 +193 1.660403493758395 3.316764122600103 0 +194 3.332833974604121 8.330588915852831 10 +195 4.990396923629646 1.639192438518225 10 +196 8.331700204977253 4.99348613953155 10 +197 1.66359950713175 4.986809017439535 10 +198 1.659476797213571 1.641166790796788 10 +199 8.331602596644176 1.661100149132812 10 +200 6.666262451736792 8.331036472916939 10 +201 1.666360512036942 8.331775733188479 10 +202 8.333176767431947 8.332461623855972 10 +203 3.324752363659667 1.640654676620429 10 +204 8.331143154549407 3.325601742294733 10 +205 8.332585110528434 6.663185298000827 10 +206 1.665690323527357 6.661884085335791 10 +207 4.999238771603533 8.329411344880642 10 +208 6.660931836157782 1.649309184852187 10 +209 3.322620058061121 3.296169448190562 10 +210 3.331386262994759 6.656799061971727 10 +211 6.662741535306179 4.984497939484297 10 +212 4.988687405622154 3.294823441600734 10 +213 6.660538599871673 3.312736539238734 10 +214 4.997593572036546 6.655043997091015 10 +215 6.66508040097065 6.65903418390889 10 +216 3.327865240662655 4.977063523995539 10 +217 4.993810674308104 4.973188365824303 10 +218 1.660298963080085 3.305559657733506 10 +219 3.323175186310005 1.639797188626136 1.666666665894576 +220 8.333333333333337 1.666666666666669 1.666666666159922 +221 8.333333332883761 8.333333333333336 5.000000000000004 +222 1.666666667399848 6.666666666666669 5.000000000000002 +223 8.333333335154675 3.333333333333338 8.333333333333334 +224 1.666666665276616 5.000000000000006 1.666666666666667 +225 8.333333332534451 8.333333333333337 8.333333333333336 +226 1.66666666735759 5.000000000000002 5.000000000000002 +227 3.325134218047016 4.973125396957162 1.666666665215716 +228 6.66479451352513 6.658229402591791 1.666666668556142 +229 4.992487513742513 4.968476001600282 1.666666662890887 +230 1.666666667636574 6.66666666666667 6.666666666666671 +231 8.333333332812531 8.333333333333334 1.666666666666669 +232 8.333333334638025 5.000000000000003 6.666666666666671 +233 6.666262451736792 8.331036472916939 8.33333333154998 +234 3.317582237776889 1.666666666712437 8.331760196778665 +235 4.963214356275155 1.666666665965598 1.668784957944797 +236 1.666666666034267 3.333333333333335 6.666666666666669 +237 6.66508040097065 6.65903418390889 8.333333333387067 +238 6.634104856166324 8.333333334715782 3.332502948194829 +239 1.666666665693501 6.666666666666671 8.333333333333336 +240 8.332074246359143 4.994289507172165 1.666666666314359 +241 1.666666666304344 6.666666666666668 3.333333333333335 +242 3.284967251501693 8.333333334781612 4.993482825700283 +243 8.332511595730896 6.662870524665738 1.666666666937509 +244 1.659476797213571 1.641166790796788 8.333333332432243 +245 6.651358102186887 8.333333331446589 6.665729135182278 +246 4.996765277527716 6.653675864145596 1.666666669892875 +247 6.660931836157782 1.649309184852187 8.333333333926896 +248 6.638537756742497 1.666666666445091 6.665974027330344 +249 3.331386262994759 6.656799061971727 8.333333334876567 +250 4.988687405622154 3.294823441600734 8.333333333444344 +251 6.660538599871673 3.312736539238734 8.333333333244161 +252 8.332031501668784 3.327655810609839 1.666666663483084 +253 1.660497912810677 1.651563432065564 1.666666666532634 +254 1.666666666666671 8.333333332815677 8.333333333333332 +255 1.623165746041504 1.6666666666034 4.99391383873798 +256 8.333333333273636 3.333333333333333 5.000000000000004 +257 3.332833974604121 8.330588915852831 8.33333333503149 +258 4.97579818627555 8.333333331156744 6.664430588732456 +259 1.660298963080085 3.305559657733506 8.33333333203548 +260 3.28326020484872 8.333333336724735 3.326681050109922 +261 8.331602596644176 1.661100149132812 8.333333335117558 +262 1.666666666603764 8.333333333333336 5.000000000000004 +263 6.665924415024502 8.330016659525276 1.666666667557128 +264 3.308721581528395 8.333333331465084 6.663488840038644 +265 3.286658535168235 1.666666666599032 6.66196704619375 +266 4.937771380378987 1.666666665951131 3.33355951443616 +267 8.326507895108154 8.333333331451792 6.666387147756299 +268 8.332585110528434 6.663185298000827 8.333333332584449 +269 4.982918784821877 1.6666666659068 8.332255412505511 +270 1.664258409514333 6.659472931888667 1.666666667233425 +271 3.322620058061121 3.296169448190562 8.33333333167775 +272 8.33333333683429 5.000000000000005 8.333333333333336 +273 8.333333333060155 6.666666666666669 5.000000000000004 +274 1.666666666327729 1.666666666666668 3.333333333333338 +275 4.99888149235128 8.328781672181401 1.666666669469232 +276 8.333333333095132 1.666666666666668 3.333333333333338 +277 8.333333334111209 6.666666666666671 6.666666666666671 +278 1.666666666231004 1.666666666666668 6.666666666666671 +279 8.319477494274667 8.333333334546891 3.333169972868349 +280 8.333333333713604 5.000000000000002 3.333333333333337 +281 8.333333335124941 3.333333333333335 6.66666666666667 +282 6.62749052424298 1.66666666623669 3.334574749423636 +283 1.66666666508557 3.333333333333335 5.000000000000002 +284 8.333333332051213 6.666666666666668 3.333333333333337 +285 1.637778017620791 8.333333342488848 3.327652077469392 +286 1.666666666089987 5.000000000000003 3.333333333333334 +287 6.62594906800776 1.666666666783335 4.999916920658317 +288 4.991559204424686 3.303632966319623 1.666666666035724 +289 1.666666665979094 3.333333333333337 3.333333333333335 +290 1.666666665657806 5.000000000000005 8.333333333333336 +291 1.666666666853695 5.000000000000003 6.666666666666669 +292 1.651821825244905 8.333333330438013 6.664257878478076 +293 3.270150200748045 1.66666666673903 4.995167900443951 +294 1.666666665802565 8.333333333333334 1.66666666666667 +295 3.332149742463294 8.329006867792472 1.666666666152793 +296 4.956018879142214 8.333333342722435 3.330340837991536 +297 4.961068569266943 8.333333337999449 4.996722349190394 +298 4.997593572036546 6.655043997091015 8.33333333485422 +299 3.32335131306381 3.303903840008758 1.666666665489136 +300 3.327865240662655 4.977063523995539 8.333333333931193 +301 4.993810674308104 4.973188365824303 8.333333333793638 +302 4.943783761023866 1.666666666230662 4.998321579025965 +303 4.990499273962526 8.333333328923789 8.332390464609082 +304 6.662704487145006 3.317607957790754 1.666666665319158 +305 6.64528602688278 1.666666667201381 1.670374112938904 +306 8.315733273001639 1.666666666751797 5.000155748895926 +307 6.663143942894656 4.984579210873332 1.666666668138241 +308 6.662741535306179 4.984497939484297 8.333333333366063 +309 8.333333332710026 3.333333333333334 3.333333333333337 +310 1.660403493758395 3.316764122600103 1.666666666364711 +311 6.641946928853725 8.33333332998944 4.998695732656651 +312 3.329655718808052 6.653762669731537 1.666666670056743 +313 8.333333334210733 1.666666666666668 6.666666666666671 +314 3.278505243836852 1.666666665982843 3.33055000927853 +315 8.333333333196052 5.000000000000001 5.000000000000003 +316 4.958530187990556 1.666666666701269 6.664632353924778 +317 3.333333335751559 6.666666666666669 5.000000000000002 +318 3.33333333945139 5.000000000000002 5.000000000000002 +319 3.325134218047016 4.973125396957162 3.333333331397721 +320 6.66479451352513 6.658229402591791 3.333333336627152 +321 4.992487513742513 4.968476001600282 3.333333328810805 +322 3.333333335727594 6.66666666666667 6.666666666666671 +323 6.66666666911012 5.000000000000003 6.666666666666671 +324 3.333333333488395 3.333333333333335 6.666666666666669 +325 6.66508040097065 6.65903418390889 6.666666667932162 +326 3.333333334384611 6.666666666666668 3.333333333333335 +327 4.996765277527716 6.653675864145596 3.333333337116327 +328 6.638537756742497 3.333333332536744 6.665974027330344 +329 4.988687405622154 3.294823441600734 6.666666666679443 +330 6.666666667093498 3.333333333333333 5.000000000000004 +331 4.97579818627555 6.666666665298601 6.664430588732456 +332 4.937771380378987 3.333333331698917 3.33355951443616 +333 6.666666667487235 6.666666666666669 5.000000000000004 +334 6.666666667036193 5.000000000000002 3.333333333333337 +335 6.62749052424298 3.333333332299917 3.334574749423636 +336 3.33333333143754 3.333333333333335 5.000000000000002 +337 3.333333331313443 3.333333333333337 3.333333333333335 +338 3.333333334724964 5.000000000000003 6.666666666666669 +339 4.961068569266943 6.666666676876886 4.996722349190394 +340 4.993810674308104 4.973188365824303 6.666666667605718 +341 4.943783761023866 3.33333333244255 4.998321579025965 +342 6.666666666438698 5.000000000000001 5.000000000000003 +343 5.000000007071594 5.000000000000002 5.000000000000002 +$EndNodes +$Elements +512 +1 15 2 0 1 1 +2 15 2 0 2 2 +3 15 2 0 3 3 +4 15 2 0 4 4 +5 15 2 0 5 5 +6 15 2 0 6 6 +7 15 2 0 7 7 +8 15 2 0 8 8 +9 1 2 0 1 1 9 +10 1 2 0 1 9 10 +11 1 2 0 1 10 11 +12 1 2 0 1 11 12 +13 1 2 0 1 12 13 +14 1 2 0 1 13 2 +15 1 2 0 2 1 14 +16 1 2 0 2 14 15 +17 1 2 0 2 15 16 +18 1 2 0 2 16 17 +19 1 2 0 2 17 18 +20 1 2 0 2 18 3 +21 1 2 0 3 3 19 +22 1 2 0 3 19 20 +23 1 2 0 3 20 21 +24 1 2 0 3 21 22 +25 1 2 0 3 22 23 +26 1 2 0 3 23 4 +27 1 2 0 4 2 24 +28 1 2 0 4 24 25 +29 1 2 0 4 25 26 +30 1 2 0 4 26 27 +31 1 2 0 4 27 28 +32 1 2 0 4 28 4 +33 1 2 0 5 5 29 +34 1 2 0 5 29 30 +35 1 2 0 5 30 31 +36 1 2 0 5 31 32 +37 1 2 0 5 32 33 +38 1 2 0 5 33 6 +39 1 2 0 6 5 34 +40 1 2 0 6 34 35 +41 1 2 0 6 35 36 +42 1 2 0 6 36 37 +43 1 2 0 6 37 38 +44 1 2 0 6 38 7 +45 1 2 0 7 7 39 +46 1 2 0 7 39 40 +47 1 2 0 7 40 41 +48 1 2 0 7 41 42 +49 1 2 0 7 42 43 +50 1 2 0 7 43 8 +51 1 2 0 8 6 44 +52 1 2 0 8 44 45 +53 1 2 0 8 45 46 +54 1 2 0 8 46 47 +55 1 2 0 8 47 48 +56 1 2 0 8 48 8 +57 1 2 0 9 1 49 +58 1 2 0 9 49 50 +59 1 2 0 9 50 51 +60 1 2 0 9 51 52 +61 1 2 0 9 52 53 +62 1 2 0 9 53 5 +63 1 2 0 10 2 54 +64 1 2 0 10 54 55 +65 1 2 0 10 55 56 +66 1 2 0 10 56 57 +67 1 2 0 10 57 58 +68 1 2 0 10 58 6 +69 1 2 0 11 3 59 +70 1 2 0 11 59 60 +71 1 2 0 11 60 61 +72 1 2 0 11 61 62 +73 1 2 0 11 62 63 +74 1 2 0 11 63 7 +75 1 2 0 12 4 64 +76 1 2 0 12 64 65 +77 1 2 0 12 65 66 +78 1 2 0 12 66 67 +79 1 2 0 12 67 68 +80 1 2 0 12 68 8 +81 3 2 0 1 82 71 88 93 +82 3 2 0 1 82 93 85 91 +83 3 2 0 1 82 84 92 91 +84 3 2 0 1 91 86 83 92 +85 3 2 0 1 86 91 85 90 +86 3 2 0 1 4 28 75 23 +87 3 2 0 1 75 28 27 88 +88 3 2 0 1 18 81 77 17 +89 3 2 0 1 77 90 86 72 +90 3 2 0 1 75 23 22 89 +91 3 2 0 1 14 15 87 79 +92 3 2 0 1 26 27 88 71 +93 3 2 0 1 17 77 72 16 +94 3 2 0 1 89 75 88 93 +95 3 2 0 1 15 87 72 16 +96 3 2 0 1 87 83 86 72 +97 3 2 0 1 89 70 85 93 +98 3 2 0 1 87 83 69 79 +99 3 2 0 1 3 18 81 19 +100 3 2 0 1 26 25 74 71 +101 3 2 0 1 25 74 80 24 +102 3 2 0 1 70 21 22 89 +103 3 2 0 1 80 13 2 24 +104 3 2 0 1 12 78 80 13 +105 3 2 0 1 19 20 73 81 +106 3 2 0 1 82 71 74 84 +107 3 2 0 1 85 70 73 90 +108 3 2 0 1 20 73 70 21 +109 3 2 0 1 74 80 78 84 +110 3 2 0 1 14 79 9 1 +111 3 2 0 1 73 90 77 81 +112 3 2 0 1 10 9 79 69 +113 3 2 0 1 78 12 11 76 +114 3 2 0 1 78 84 92 76 +115 3 2 0 1 69 10 11 76 +116 3 2 0 1 69 76 92 83 +117 3 2 0 2 114 116 107 118 +118 3 2 0 2 108 110 117 116 +119 3 2 0 2 115 114 116 108 +120 3 2 0 2 96 100 115 114 +121 3 2 0 2 117 109 107 116 +122 3 2 0 2 48 98 112 47 +123 3 2 0 2 8 43 98 48 +124 3 2 0 2 114 118 112 96 +125 3 2 0 2 38 37 102 106 +126 3 2 0 2 111 35 34 104 +127 3 2 0 2 35 111 97 36 +128 3 2 0 2 37 102 97 36 +129 3 2 0 2 110 94 104 111 +130 3 2 0 2 98 43 42 113 +131 3 2 0 2 110 111 97 117 +132 3 2 0 2 97 117 109 102 +133 3 2 0 2 100 105 103 115 +134 3 2 0 2 113 95 107 118 +135 3 2 0 2 6 33 105 44 +136 3 2 0 2 118 113 98 112 +137 3 2 0 2 96 112 47 46 +138 3 2 0 2 45 44 105 100 +139 3 2 0 2 32 33 105 103 +140 3 2 0 2 45 46 96 100 +141 3 2 0 2 7 39 106 38 +142 3 2 0 2 41 95 113 42 +143 3 2 0 2 34 5 29 104 +144 3 2 0 2 106 39 40 99 +145 3 2 0 2 95 99 40 41 +146 3 2 0 2 94 30 29 104 +147 3 2 0 2 102 109 99 106 +148 3 2 0 2 107 109 99 95 +149 3 2 0 2 31 32 103 101 +150 3 2 0 2 94 30 31 101 +151 3 2 0 2 115 108 101 103 +152 3 2 0 2 110 94 101 108 +153 3 2 0 3 136 141 137 138 +154 3 2 0 3 143 135 137 141 +155 3 2 0 3 6 58 128 33 +156 3 2 0 3 127 128 33 32 +157 3 2 0 3 53 5 29 125 +158 3 2 0 3 141 140 131 136 +159 3 2 0 3 30 29 125 130 +160 3 2 0 3 57 58 128 124 +161 3 2 0 3 143 141 140 139 +162 3 2 0 3 31 121 127 32 +163 3 2 0 3 136 121 127 131 +164 3 2 0 3 130 138 134 125 +165 3 2 0 3 125 134 52 53 +166 3 2 0 3 30 31 121 130 +167 3 2 0 3 137 138 134 119 +168 3 2 0 3 136 138 130 121 +169 3 2 0 3 133 56 55 120 +170 3 2 0 3 124 133 140 131 +171 3 2 0 3 126 13 2 54 +172 3 2 0 3 120 133 140 139 +173 3 2 0 3 124 133 56 57 +174 3 2 0 3 127 131 124 128 +175 3 2 0 3 12 13 126 132 +176 3 2 0 3 126 54 55 120 +177 3 2 0 3 52 51 119 134 +178 3 2 0 3 139 132 126 120 +179 3 2 0 3 1 49 123 9 +180 3 2 0 3 123 49 50 142 +181 3 2 0 3 143 122 132 139 +182 3 2 0 3 10 9 123 129 +183 3 2 0 3 132 12 11 122 +184 3 2 0 3 142 119 51 50 +185 3 2 0 3 123 142 135 129 +186 3 2 0 3 129 122 11 10 +187 3 2 0 3 129 135 143 122 +188 3 2 0 3 119 137 135 142 +189 3 2 0 4 4 64 149 23 +190 3 2 0 4 68 152 43 8 +191 3 2 0 4 152 43 42 150 +192 3 2 0 4 162 167 166 159 +193 3 2 0 4 67 156 152 68 +194 3 2 0 4 166 167 163 160 +195 3 2 0 4 67 66 146 156 +196 3 2 0 4 154 65 64 149 +197 3 2 0 4 154 65 66 146 +198 3 2 0 4 7 63 151 39 +199 3 2 0 4 157 150 165 163 +200 3 2 0 4 157 41 42 150 +201 3 2 0 4 23 22 158 149 +202 3 2 0 4 156 165 150 152 +203 3 2 0 4 39 151 144 40 +204 3 2 0 4 165 163 167 161 +205 3 2 0 4 146 156 165 161 +206 3 2 0 4 157 144 40 41 +207 3 2 0 4 160 144 157 163 +208 3 2 0 4 59 148 19 3 +209 3 2 0 4 146 161 164 154 +210 3 2 0 4 158 149 154 164 +211 3 2 0 4 161 164 162 167 +212 3 2 0 4 145 158 22 21 +213 3 2 0 4 151 63 62 155 +214 3 2 0 4 162 145 158 164 +215 3 2 0 4 145 162 159 153 +216 3 2 0 4 62 61 147 155 +217 3 2 0 4 19 148 153 20 +218 3 2 0 4 61 147 168 60 +219 3 2 0 4 144 151 155 160 +220 3 2 0 4 60 168 148 59 +221 3 2 0 4 153 145 21 20 +222 3 2 0 4 155 160 166 147 +223 3 2 0 4 147 166 159 168 +224 3 2 0 4 168 148 153 159 +225 3 2 0 5 189 180 171 186 +226 3 2 0 5 187 185 192 188 +227 3 2 0 5 189 184 175 180 +228 3 2 0 5 189 188 192 186 +229 3 2 0 5 38 37 181 174 +230 3 2 0 5 169 35 34 176 +231 3 2 0 5 169 177 36 35 +232 3 2 0 5 183 191 181 174 +233 3 2 0 5 37 181 177 36 +234 3 2 0 5 183 191 186 171 +235 3 2 0 5 18 17 184 175 +236 3 2 0 5 169 178 190 177 +237 3 2 0 5 177 190 191 181 +238 3 2 0 5 7 63 174 38 +239 3 2 0 5 187 193 173 179 +240 3 2 0 5 176 53 5 34 +241 3 2 0 5 14 15 193 173 +242 3 2 0 5 185 178 190 192 +243 3 2 0 5 17 184 170 16 +244 3 2 0 5 193 15 16 170 +245 3 2 0 5 192 190 191 186 +246 3 2 0 5 53 176 182 52 +247 3 2 0 5 59 175 18 3 +248 3 2 0 5 193 170 188 187 +249 3 2 0 5 182 178 169 176 +250 3 2 0 5 174 63 62 183 +251 3 2 0 5 170 188 189 184 +252 3 2 0 5 62 61 171 183 +253 3 2 0 5 61 171 180 60 +254 3 2 0 5 52 182 172 51 +255 3 2 0 5 60 180 175 59 +256 3 2 0 5 14 1 49 173 +257 3 2 0 5 49 50 179 173 +258 3 2 0 5 182 172 185 178 +259 3 2 0 5 50 51 172 179 +260 3 2 0 5 187 185 172 179 +261 3 2 0 6 210 216 217 214 +262 3 2 0 6 4 28 201 64 +263 3 2 0 6 216 217 212 209 +264 3 2 0 6 48 47 205 202 +265 3 2 0 6 68 8 48 202 +266 3 2 0 6 215 200 202 205 +267 3 2 0 6 201 28 27 206 +268 3 2 0 6 200 67 68 202 +269 3 2 0 6 6 58 199 44 +270 3 2 0 6 46 196 205 47 +271 3 2 0 6 207 194 210 214 +272 3 2 0 6 67 200 207 66 +273 3 2 0 6 45 204 199 44 +274 3 2 0 6 210 194 201 206 +275 3 2 0 6 65 64 201 194 +276 3 2 0 6 204 196 46 45 +277 3 2 0 6 65 194 207 66 +278 3 2 0 6 215 200 207 214 +279 3 2 0 6 214 217 211 215 +280 3 2 0 6 197 206 27 26 +281 3 2 0 6 205 215 211 196 +282 3 2 0 6 208 199 58 57 +283 3 2 0 6 209 216 197 218 +284 3 2 0 6 213 211 196 204 +285 3 2 0 6 213 212 217 211 +286 3 2 0 6 213 208 199 204 +287 3 2 0 6 197 206 210 216 +288 3 2 0 6 26 25 218 197 +289 3 2 0 6 25 218 198 24 +290 3 2 0 6 212 209 203 195 +291 3 2 0 6 195 208 213 212 +292 3 2 0 6 56 195 203 55 +293 3 2 0 6 198 24 2 54 +294 3 2 0 6 208 195 56 57 +295 3 2 0 6 54 198 203 55 +296 3 2 0 6 209 218 198 203 +297 5 2 0 1 114 118 107 116 232 277 273 315 +298 5 2 0 1 108 116 117 110 256 315 280 309 +299 5 2 0 1 281 232 114 115 256 315 116 108 +300 5 2 0 1 82 291 230 93 71 290 239 88 +301 5 2 0 1 291 230 222 226 82 93 85 91 +302 5 2 0 1 318 338 291 226 336 324 236 283 +303 5 2 0 1 82 84 92 91 291 236 283 226 +304 5 2 0 1 226 283 92 91 286 289 83 86 +305 5 2 0 1 96 100 223 272 114 115 281 232 +306 5 2 0 1 286 226 222 241 86 91 85 90 +307 5 2 0 1 280 117 109 284 315 116 107 273 +308 5 2 0 1 318 226 291 338 317 222 230 322 +309 5 2 0 1 249 298 214 210 300 301 217 216 +310 5 2 0 1 4 23 149 64 28 75 254 201 +311 5 2 0 1 312 189 180 295 246 186 171 275 +312 5 2 0 1 187 299 288 185 188 227 229 192 +313 5 2 0 1 216 209 271 300 217 212 250 301 +314 5 2 0 1 312 189 184 270 295 180 175 294 +315 5 2 0 1 48 202 225 98 47 205 268 112 +316 5 2 0 1 68 202 225 152 8 48 98 43 +317 5 2 0 1 189 312 227 188 186 246 229 192 +318 5 2 0 1 114 118 277 232 96 112 268 272 +319 5 2 0 1 38 174 181 37 106 231 243 102 +320 5 2 0 1 215 200 202 205 237 233 225 268 +321 5 2 0 1 252 169 35 111 220 176 34 104 +322 5 2 0 1 169 177 240 252 35 36 97 111 +323 5 2 0 1 183 263 228 191 174 231 243 181 +324 5 2 0 1 37 102 243 181 36 97 240 177 +325 5 2 0 1 254 75 88 239 201 28 27 206 +326 5 2 0 1 183 171 186 191 263 275 246 228 +327 5 2 0 1 18 81 77 17 175 294 270 184 +328 5 2 0 1 309 110 94 276 252 111 104 220 +329 5 2 0 1 225 152 43 98 267 150 42 113 +330 5 2 0 1 287 282 138 136 302 266 137 141 +331 5 2 0 1 110 111 252 309 117 97 240 280 +332 5 2 0 1 97 240 280 117 102 243 284 109 +333 5 2 0 1 162 167 297 242 159 166 296 260 +334 5 2 0 1 143 293 302 141 135 314 266 137 +335 5 2 0 1 200 202 225 233 67 68 152 156 +336 5 2 0 1 100 105 261 223 115 103 313 281 +337 5 2 0 1 267 277 273 221 113 118 107 95 +338 5 2 0 1 6 44 105 33 58 199 261 128 +339 5 2 0 1 118 112 98 113 277 268 225 267 +340 5 2 0 1 315 232 277 273 342 323 325 333 +341 5 2 0 1 96 46 47 112 272 196 205 268 +342 5 2 0 1 166 167 297 296 160 163 311 238 +343 5 2 0 1 303 298 214 207 257 249 210 194 +344 5 2 0 1 252 169 178 304 240 177 190 307 +345 5 2 0 1 67 200 207 66 156 233 303 146 +346 5 2 0 1 45 100 105 44 204 223 261 199 +347 5 2 0 1 315 273 284 280 342 333 320 334 +348 5 2 0 1 249 239 206 210 257 254 201 194 +349 5 2 0 1 154 257 194 65 149 254 201 64 +350 5 2 0 1 240 177 190 307 243 181 191 228 +351 5 2 0 1 127 128 33 32 313 261 105 103 +352 5 2 0 1 204 196 46 45 223 272 96 100 +353 5 2 0 1 154 146 66 65 257 303 207 194 +354 5 2 0 1 77 72 224 270 90 86 286 241 +355 5 2 0 1 7 39 151 63 38 106 231 174 +356 5 2 0 1 299 310 193 187 219 253 173 179 +357 5 2 0 1 157 221 267 150 163 311 245 165 +358 5 2 0 1 157 221 95 41 150 267 113 42 +359 5 2 0 1 237 298 214 215 233 303 207 200 +360 5 2 0 1 326 327 321 319 312 246 229 227 +361 5 2 0 1 176 220 125 53 34 104 29 5 +362 5 2 0 1 75 254 292 89 23 149 158 22 +363 5 2 0 1 298 237 215 214 301 308 211 217 +364 5 2 0 1 233 245 165 156 225 267 150 152 +365 5 2 0 1 141 140 131 136 302 316 248 287 +366 5 2 0 1 106 39 151 231 99 40 144 279 +367 5 2 0 1 165 161 167 163 245 258 297 311 +368 5 2 0 1 14 15 87 79 173 193 310 253 +369 5 2 0 1 288 304 178 185 229 307 190 192 +370 5 2 0 1 197 290 71 26 206 239 88 27 +371 5 2 0 1 268 272 196 205 237 308 211 215 +372 5 2 0 1 146 156 233 303 161 165 245 258 +373 5 2 0 1 221 95 99 279 157 41 40 144 +374 5 2 0 1 94 30 130 276 104 29 125 220 +375 5 2 0 1 208 57 124 247 199 58 128 261 +376 5 2 0 1 17 16 170 184 77 72 224 270 +377 5 2 0 1 160 238 279 144 163 311 221 157 +378 5 2 0 1 280 284 243 240 334 320 228 307 +379 5 2 0 1 243 284 109 102 231 279 99 106 +380 5 2 0 1 89 93 230 292 75 88 239 254 +381 5 2 0 1 293 302 316 265 143 141 140 139 +382 5 2 0 1 273 107 109 284 221 95 99 279 +383 5 2 0 1 193 170 16 15 310 224 72 87 +384 5 2 0 1 31 121 127 32 101 306 313 103 +385 5 2 0 1 192 229 307 190 186 246 228 191 +386 5 2 0 1 310 224 72 87 289 286 86 83 +387 5 2 0 1 287 248 131 136 306 313 127 121 +388 5 2 0 1 277 268 225 267 325 237 233 245 +389 5 2 0 1 130 276 220 125 138 282 305 134 +390 5 2 0 1 271 259 218 209 300 290 197 216 +391 5 2 0 1 251 213 204 223 308 211 196 272 +392 5 2 0 1 292 262 222 230 89 70 85 93 +393 5 2 0 1 213 211 308 251 212 217 301 250 +394 5 2 0 1 213 251 223 204 208 247 261 199 +395 5 2 0 1 268 272 308 237 277 232 323 325 +396 5 2 0 1 125 220 176 53 134 305 182 52 +397 5 2 0 1 94 276 130 30 101 306 121 31 +398 5 2 0 1 87 83 69 79 310 289 274 253 +399 5 2 0 1 137 119 235 266 138 134 305 282 +400 5 2 0 1 59 175 294 148 3 18 81 19 +401 5 2 0 1 281 256 108 115 313 306 101 103 +402 5 2 0 1 309 256 306 276 110 108 101 94 +403 5 2 0 1 303 146 161 258 257 154 164 264 +404 5 2 0 1 287 282 276 306 136 138 130 121 +405 5 2 0 1 193 170 224 310 187 188 227 299 +406 5 2 0 1 292 158 149 254 264 164 154 257 +407 5 2 0 1 290 197 216 300 239 206 210 249 +408 5 2 0 1 258 264 164 161 297 242 162 167 +409 5 2 0 1 277 273 333 325 267 221 311 245 +410 5 2 0 1 297 339 333 311 258 331 325 245 +411 5 2 0 1 258 245 233 303 331 325 237 298 +412 5 2 0 1 26 71 290 197 25 74 259 218 +413 5 2 0 1 25 218 198 24 74 259 244 80 +414 5 2 0 1 250 269 195 212 271 234 203 209 +415 5 2 0 1 262 70 89 292 145 21 22 158 +416 5 2 0 1 182 176 169 178 305 220 252 304 +417 5 2 0 1 270 224 227 312 241 286 319 326 +418 5 2 0 1 195 269 250 212 208 247 251 213 +419 5 2 0 1 231 151 63 174 263 155 62 183 +420 5 2 0 1 133 269 234 120 56 195 203 55 +421 5 2 0 1 300 290 239 249 338 291 230 322 +422 5 2 0 1 247 269 316 248 124 133 140 131 +423 5 2 0 1 221 279 284 273 311 238 320 333 +424 5 2 0 1 242 162 145 262 264 164 158 292 +425 5 2 0 1 244 126 54 198 80 13 2 24 +426 5 2 0 1 224 227 188 170 270 312 189 184 +427 5 2 0 1 145 162 242 262 153 159 260 285 +428 5 2 0 1 234 120 139 265 269 133 140 316 +429 5 2 0 1 247 124 57 208 269 133 56 195 +430 5 2 0 1 127 128 261 313 131 124 247 248 +431 5 2 0 1 12 13 126 132 78 80 244 278 +432 5 2 0 1 317 326 319 318 222 241 286 226 +433 5 2 0 1 62 155 147 61 183 263 275 171 +434 5 2 0 1 239 249 322 230 254 257 264 292 +435 5 2 0 1 19 81 294 148 20 73 285 153 +436 5 2 0 1 286 224 227 319 289 310 299 337 +437 5 2 0 1 264 258 303 257 322 331 298 249 +438 5 2 0 1 334 320 228 307 321 327 246 229 +439 5 2 0 1 291 82 84 236 290 71 74 259 +440 5 2 0 1 223 281 313 261 251 328 248 247 +441 5 2 0 1 61 171 275 147 60 180 295 168 +442 5 2 0 1 144 160 155 151 279 238 263 231 +443 5 2 0 1 126 54 55 120 244 198 203 234 +444 5 2 0 1 52 134 305 182 51 119 235 172 +445 5 2 0 1 85 222 241 90 70 262 285 73 +446 5 2 0 1 323 325 331 340 308 237 298 301 +447 5 2 0 1 60 180 295 168 59 175 294 148 +448 5 2 0 1 153 285 262 145 20 73 70 21 +449 5 2 0 1 271 300 338 324 259 290 291 236 +450 5 2 0 1 139 120 234 265 132 126 244 278 +451 5 2 0 1 74 80 78 84 259 244 278 236 +452 5 2 0 1 342 330 256 315 323 328 281 232 +453 5 2 0 1 271 259 244 234 209 218 198 203 +454 5 2 0 1 308 323 232 272 251 328 281 223 +455 5 2 0 1 14 79 9 1 173 253 123 49 +456 5 2 0 1 123 253 219 142 49 173 179 50 +457 5 2 0 1 238 263 231 279 320 228 243 284 +458 5 2 0 1 322 249 300 338 331 298 301 340 +459 5 2 0 1 182 172 235 305 178 185 288 304 +460 5 2 0 1 337 336 318 319 289 283 226 286 +461 5 2 0 1 143 122 132 139 293 255 278 265 +462 5 2 0 1 73 81 294 285 90 77 270 241 +463 5 2 0 1 10 69 274 129 9 79 253 123 +464 5 2 0 1 304 335 282 305 252 309 276 220 +465 5 2 0 1 263 238 160 155 275 296 166 147 +466 5 2 0 1 287 306 256 330 248 313 281 328 +467 5 2 0 1 342 333 339 343 323 325 331 340 +468 5 2 0 1 278 78 76 255 132 12 11 122 +469 5 2 0 1 317 326 241 222 242 260 285 262 +470 5 2 0 1 339 333 320 327 297 311 238 296 +471 5 2 0 1 317 242 262 222 322 264 292 230 +472 5 2 0 1 321 327 339 343 334 320 333 342 +473 5 2 0 1 322 331 340 338 317 339 343 318 +474 5 2 0 1 142 119 51 50 219 235 172 179 +475 5 2 0 1 280 315 342 334 309 256 330 335 +476 5 2 0 1 329 340 338 324 250 301 300 271 +477 5 2 0 1 253 274 129 123 219 314 135 142 +478 5 2 0 1 335 334 280 309 304 307 240 252 +479 5 2 0 1 275 295 260 296 147 168 159 166 +480 5 2 0 1 339 297 242 317 331 258 264 322 +481 5 2 0 1 317 318 319 326 339 343 321 327 +482 5 2 0 1 287 282 335 330 306 276 309 256 +483 5 2 0 1 246 228 263 275 327 320 238 296 +484 5 2 0 1 278 255 76 78 236 283 92 84 +485 5 2 0 1 295 168 148 294 260 159 153 285 +486 5 2 0 1 329 324 265 316 250 271 234 269 +487 5 2 0 1 330 287 302 341 335 282 266 332 +488 5 2 0 1 244 234 265 278 259 271 324 236 +489 5 2 0 1 274 255 76 69 129 122 11 10 +490 5 2 0 1 274 255 122 129 314 293 143 135 +491 5 2 0 1 341 336 337 332 302 293 314 266 +492 5 2 0 1 293 336 337 314 255 283 289 274 +493 5 2 0 1 327 296 260 326 339 297 242 317 +494 5 2 0 1 187 185 288 299 179 172 235 219 +495 5 2 0 1 337 336 341 332 319 318 343 321 +496 5 2 0 1 336 283 255 293 324 236 278 265 +497 5 2 0 1 293 336 324 265 302 341 329 316 +498 5 2 0 1 274 255 283 289 69 76 92 83 +499 5 2 0 1 260 295 275 296 326 312 246 327 +500 5 2 0 1 316 329 250 269 248 328 251 247 +501 5 2 0 1 270 294 285 241 312 295 260 326 +502 5 2 0 1 329 340 301 250 328 323 308 251 +503 5 2 0 1 332 288 304 335 266 235 305 282 +504 5 2 0 1 235 219 142 119 266 314 135 137 +505 5 2 0 1 337 299 288 332 314 219 235 266 +506 5 2 0 1 310 299 337 289 253 219 314 274 +507 5 2 0 1 343 341 329 340 318 336 324 338 +508 5 2 0 1 329 328 330 341 340 323 342 343 +509 5 2 0 1 337 332 288 299 319 321 229 227 +510 5 2 0 1 304 288 332 335 307 229 321 334 +511 5 2 0 1 287 330 341 302 248 328 329 316 +512 5 2 0 1 334 342 343 321 335 330 341 332 +$EndElements From 06711ac90cbfba79f58047990729e2be36ac8a19 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 6 Nov 2020 09:28:10 -0600 Subject: [PATCH 14/51] rework tensor product inheritance a bit --- meshmode/discretization/poly_element.py | 99 +++++++++++++------------ 1 file changed, 53 insertions(+), 46 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index ac54fb833..1d37bcbc6 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -42,7 +42,7 @@ .. autoclass:: PolynomialEquidistantSimplexElementGroup .. autoclass:: PolynomialGivenNodesElementGroup -.. autoclass:: InterpolatoryQuadratureTensorProductElementGroup +.. autoclass:: GaussLegendreTensorProductElementGroup .. autoclass:: LegendreGaussLobattoTensorProductElementGroup .. autoclass:: EquidistantTensorProductElementGroup @@ -59,7 +59,7 @@ .. autoclass:: PolynomialEquidistantSimplexGroupFactory .. autoclass:: PolynomialGivenNodesGroupFactory -.. autoclass:: InterpolatoryQuadratureTensorProductGroupFactory +.. autoclass:: GaussLegendreTensorProductGroupFactory .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ @@ -353,40 +353,33 @@ def unit_nodes(self): # {{{ concrete element groups for tensor product elements class _TensorProductElementGroupBase(PolynomialElementGroupBase): - # {{{ 1D basis - - def _basis_1d(self): - from modepy.modes import jacobi - from functools import partial - return tuple(partial(jacobi, 0, 0, i) for i in range(self.order + 1)) - - def _grad_basis_1d(self): - from modepy.modes import grad_jacobi - from functools import partial - return tuple(partial(grad_jacobi, 0, 0, i) for i in range(self.order + 1)) - - def _unit_nodes_1d(self): - raise NotImplementedError + def __init__(self, mesh_el_group, order, index, *, + basis_1d=None, grad_basis_1d=None, + unit_nodes_1d=None, quad_weights_1d=None): + super().__init__(mesh_el_group, order, index) - # }}} + self._basis_1d = basis_1d + self._grad_basis_1d = grad_basis_1d + self._unit_nodes_1d = unit_nodes_1d + self._quad_weights_1d = quad_weights_1d def is_orthogonal_basis(self): return True def basis(self): from modepy.modes import tensor_product_basis - return tensor_product_basis(self.dim, self._basis_1d()) + return tensor_product_basis(self.dim, self._basis_1d) def grad_basis(self): from modepy.modes import grad_tensor_product_basis return grad_tensor_product_basis(self.dim, - self._basis_1d(), self._grad_basis_1d()) + self._basis_1d, self._grad_basis_1d) @property @memoize_method def unit_nodes(self): from modepy.nodes import tensor_product_nodes - return tensor_product_nodes(self.dim, self._unit_nodes_1d()) + return tensor_product_nodes(self.dim, self._unit_nodes_1d) @memoize_method def from_mesh_interp_matrix(self): @@ -397,36 +390,49 @@ def from_mesh_interp_matrix(self): self.unit_nodes, meg.unit_nodes) + @property + def weights(self): + from itertools import product + weights = self._quad_weights_1d + return np.fromiter( + (np.prod(w) for w in product(weights, repeat=self.dim)), + dtype=np.float, + count=(self.order + 1)**self.dim) + -class InterpolatoryQuadratureTensorProductElementGroup( - _TensorProductElementGroupBase): +class _LegendreTensorProductElementGroup(_TensorProductElementGroupBase): + def __init__(self, mesh_el_group, order, index, *, + unit_nodes_1d=None, quad_weights_1d=None): + from modepy.modes import jacobi, grad_jacobi + from functools import partial + + super().__init__(mesh_el_group, order, index, + basis_1d=tuple( + partial(jacobi, 0, 0, i) for i in range(order + 1)), + grad_basis_1d=tuple( + partial(grad_jacobi, 0, 0, i) for i in range(order + 1)), + unit_nodes_1d=unit_nodes_1d, + quad_weights_1d=quad_weights_1d) + + +class GaussLegendreTensorProductElementGroup(_LegendreTensorProductElementGroup): """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. No interpolation nodes are present on the boundary of the hypercube. """ - - @memoize_method - def _quadrature_rule_1d(self): + def __init__(self, mesh_el_group, order, index): import modepy as mp - return mp.LegendreGaussQuadrature(self.order) - - def _unit_nodes_1d(self): - return self._quadrature_rule_1d().nodes + quad = mp.LegendreGaussQuadrature(self.order) - @property - def weights(self): - from itertools import product - weights = self._quadrature_rule_1d().weights - return np.fromiter( - (np.prod(w) for w in product(weights, repeat=self.dim)), - dtype=np.float, - count=(self.order + 1)**self.dim) + super().__init__(mesh_el_group, order, index, + unit_nodes_1d=quad.nodes, + quad_nodes_1d=quad.weights) class LegendreGaussLobattoTensorProductElementGroup( - _TensorProductElementGroupBase): + _LegendreTensorProductElementGroup): """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. @@ -435,12 +441,13 @@ class LegendreGaussLobattoTensorProductElementGroup( Uses :func:`~modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes`. """ - def _unit_nodes_1d(self): + def __init__(self, mesh_el_group, order, index): from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes - return legendre_gauss_lobatto_nodes(self.order) + super().__init__(mesh_el_group, order, index, + unit_nodes_1d=legendre_gauss_lobatto_nodes(order)) -class EquidistantTensorProductElementGroup(_TensorProductElementGroupBase): +class EquidistantTensorProductElementGroup(_LegendreTensorProductElementGroup): """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. @@ -449,9 +456,10 @@ class EquidistantTensorProductElementGroup(_TensorProductElementGroupBase): Uses :func:`~modepy.equidistant_nodes`. """ - def _unit_nodes_1d(self): + def __init__(self, mesh_el_group, order, index): from modepy.nodes import equidistant_nodes - return equidistant_nodes(1, self.order)[0] + super().__init__(mesh_el_group, order, index, + unit_nodes_1d=equidistant_nodes(1, order)[0]) # }}} @@ -555,10 +563,9 @@ def __call__(self, mesh_el_group, index): # {{{ group factories for tensor products -class InterpolatoryQuadratureTensorProductGroupFactory( - HomogeneousOrderBasedGroupFactory): +class GaussLegendreTensorProductGroupFactory(HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshTensorProductElementGroup - group_class = InterpolatoryQuadratureTensorProductElementGroup + group_class = GaussLegendreTensorProductElementGroup class LegendreGaussLobattoTensorProductGroupFactory( From 3356de65cbec592059d1897102663fe3b192eb0f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 7 Nov 2020 15:06:54 -0600 Subject: [PATCH 15/51] add gmsh quad test --- test/test_meshmode.py | 77 ++++++++++++++++++++++++++++++------------- 1 file changed, 54 insertions(+), 23 deletions(-) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index e853c5c62..73bb0345e 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1261,32 +1261,63 @@ def test_quad_mesh_2d(ambient_dim, filename, visualize=False): # {{{ test_quad_mesh_3d -# This currently (gmsh 2.13.2) crashes gmsh. A massaged version of this using -# 'cube.step' succeeded in generating 'hybrid-cube.msh' and 'cubed-cube.msh'. -def no_test_quad_mesh_3d(): - from meshmode.mesh.io import generate_gmsh, ScriptWithFilesSource - print("BEGIN GEN") - mesh = generate_gmsh( - ScriptWithFilesSource( - """ - Merge "ball-radius-1.step"; - // Mesh.CharacteristicLengthMax = 0.1; +@pytest.mark.parametrize("mesh_name", [ + # this currently (2020-11-05 with gmsh 4.6.0) does not recombine anything + # or flat out crashes gmsh + # "ball", - Mesh.RecombineAll=1; - Mesh.Recombine3DAll=1; - Mesh.Algorithm = 8; - Mesh.Algorithm3D = 9; - // Mesh.Smoothing = 0; + "cube", + ]) +def test_quad_mesh_3d(mesh_name, visualize=False): + if mesh_name == "ball": + from meshmode.mesh.io import ScriptWithFilesSource + script = ScriptWithFilesSource( + """ + Merge "ball-radius-1.step"; + // Mesh.CharacteristicLengthMax = 0.1; + + Mesh.RecombineAll = 1; + Mesh.Recombine3DAll = 1; + Mesh.Recombine3DLevel = 2; + + Mesh.Algorithm = 8; + Mesh.Algorithm3D = 8; + + Mesh 3; + Save "output.msh"; + """, + ["ball-radius-1.step"]) + + with open("ball-quad.geo", "w") as f: + f.write(script.source) + + elif mesh_name == "cube": + from meshmode.mesh.io import ScriptSource + script = ScriptSource( + """ + SetFactory("OpenCASCADE"); + Box(1) = {0, 0, 0, 1, 1, 1}; + + Transfinite Line "*" = 8; + Transfinite Surface "*"; + Transfinite Volume "*"; + + Mesh.RecombineAll = 1; + Mesh.Recombine3DAll = 1; + Mesh.Recombine3DLevel = 2; + """, "geo") + else: + raise ValueError(f"unknown mesh name: '{mesh_name}'") - // Mesh.ElementOrder = 3; + from meshmode.mesh.io import generate_gmsh + logger.info("BEGIN GEN") + # FIXME: this fails when order = 2 + mesh = generate_gmsh(script, 3, order=1, target_unit="MM") + logger.info("END GEN") - Mesh 3; - Save "output.msh"; - """, - ["ball-radius-1.step"]), - ) - print("END GEN") - print(mesh.nelements) + if visualize: + from meshmode.mesh.visualization import write_vertex_vtk_file + write_vertex_vtk_file(mesh, f"quad_mesh_3d_{mesh_name}.vtu", overwrite=True) # }}} From 165b8440ac75daf6e5a9b4c059619f287db179d6 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 8 Nov 2020 16:04:44 -0600 Subject: [PATCH 16/51] move vertex shuffle to separate function --- meshmode/mesh/io.py | 26 ++++++++++++++------------ test/test_meshmode.py | 7 ++++--- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index cf34f407c..cabbe043d 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -46,6 +46,16 @@ # {{{ gmsh receiver +def _gmsh_node_shuffle(cls, order): + group = cls(order=order) + gmsh_node_dict = {gvt: i for i, gvt in enumerate(group.gmsh_node_tuples())} + + from pytools import generate_nonnegative_integer_tuples_below as gnitb + return np.array([ + gmsh_node_dict[nt] for nt in gnitb(order + 1, group.dimensions) + ]) + + class GmshMeshReceiver(GmshMeshReceiverBase): def __init__(self, mesh_construction_kwargs=None): # Use data fields similar to meshpy.triangle.MeshInfo and @@ -191,8 +201,9 @@ def get_mesh(self): continue nodes[:, i] = self.points[el_nodes].T - vertex_indices[i] = [vertex_gmsh_index_to_mine[v_nr] - for v_nr in el_vertices] + vertex_indices[i] = [ + vertex_gmsh_index_to_mine[v_nr] for v_nr in el_vertices + ] i += 1 @@ -213,16 +224,7 @@ def get_mesh(self): np.ones(ngroup_elements, np.bool)) elif isinstance(group_el_type, GmshTensorProductElementBase): - gmsh_vertex_tuples = type(group_el_type)(order=1).gmsh_node_tuples() - gmsh_vertex_tuples_loc_dict = { - gvt: i - for i, gvt in enumerate(gmsh_vertex_tuples)} - - from pytools import ( - generate_nonnegative_integer_tuples_below as gnitb) - vertex_shuffle = np.array([ - gmsh_vertex_tuples_loc_dict[vt] - for vt in gnitb(2, group_el_type.dimensions)]) + vertex_shuffle = _gmsh_node_shuffle(type(group_el_type), 1) group = TensorProductElementGroup( group_el_type.order, diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 73bb0345e..5c2feaf21 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1268,7 +1268,7 @@ def test_quad_mesh_2d(ambient_dim, filename, visualize=False): "cube", ]) -def test_quad_mesh_3d(mesh_name, visualize=False): +def test_quad_mesh_3d(mesh_name, order=1, visualize=False): if mesh_name == "ball": from meshmode.mesh.io import ScriptWithFilesSource script = ScriptWithFilesSource( @@ -1298,7 +1298,7 @@ def test_quad_mesh_3d(mesh_name, visualize=False): SetFactory("OpenCASCADE"); Box(1) = {0, 0, 0, 1, 1, 1}; - Transfinite Line "*" = 8; + Transfinite Line "*" = 1; Transfinite Surface "*"; Transfinite Volume "*"; @@ -1309,10 +1309,11 @@ def test_quad_mesh_3d(mesh_name, visualize=False): else: raise ValueError(f"unknown mesh name: '{mesh_name}'") + np.set_printoptions(linewidth=200) from meshmode.mesh.io import generate_gmsh logger.info("BEGIN GEN") # FIXME: this fails when order = 2 - mesh = generate_gmsh(script, 3, order=1, target_unit="MM") + mesh = generate_gmsh(script, 3, order=order, target_unit="MM") logger.info("END GEN") if visualize: From 468a62548690a3645e6dd3d75941247351d68d37 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 9 Nov 2020 13:20:57 -0600 Subject: [PATCH 17/51] generate high-order gmsh quad mesh in test --- meshmode/mesh/io.py | 13 ++----------- requirements.txt | 2 +- test/test_meshmode.py | 3 +-- 3 files changed, 4 insertions(+), 14 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index cabbe043d..7c345b04c 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -46,16 +46,6 @@ # {{{ gmsh receiver -def _gmsh_node_shuffle(cls, order): - group = cls(order=order) - gmsh_node_dict = {gvt: i for i, gvt in enumerate(group.gmsh_node_tuples())} - - from pytools import generate_nonnegative_integer_tuples_below as gnitb - return np.array([ - gmsh_node_dict[nt] for nt in gnitb(order + 1, group.dimensions) - ]) - - class GmshMeshReceiver(GmshMeshReceiverBase): def __init__(self, mesh_construction_kwargs=None): # Use data fields similar to meshpy.triangle.MeshInfo and @@ -224,7 +214,8 @@ def get_mesh(self): np.ones(ngroup_elements, np.bool)) elif isinstance(group_el_type, GmshTensorProductElementBase): - vertex_shuffle = _gmsh_node_shuffle(type(group_el_type), 1) + vertex_shuffle = type(group_el_type)( + order=1).get_lexicographic_gmsh_node_indices() group = TensorProductElementGroup( group_el_type.order, diff --git a/requirements.txt b/requirements.txt index 7027fed52..2f1803a26 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ numpy recursivenodes git+https://github.com/inducer/pytools.git#egg=pytools -git+https://github.com/inducer/gmsh_interop.git#egg=gmsh_interop +git+https://github.com/alexfikl/gmsh_interop.git@contrib-node-order#egg=gmsh_interop git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile git+https://github.com/inducer/modepy.git#egg=modepy git+https://github.com/inducer/pyopencl.git#egg=pyopencl diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 5c2feaf21..c212b011b 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1268,7 +1268,7 @@ def test_quad_mesh_2d(ambient_dim, filename, visualize=False): "cube", ]) -def test_quad_mesh_3d(mesh_name, order=1, visualize=False): +def test_quad_mesh_3d(mesh_name, order=3, visualize=False): if mesh_name == "ball": from meshmode.mesh.io import ScriptWithFilesSource script = ScriptWithFilesSource( @@ -1312,7 +1312,6 @@ def test_quad_mesh_3d(mesh_name, order=1, visualize=False): np.set_printoptions(linewidth=200) from meshmode.mesh.io import generate_gmsh logger.info("BEGIN GEN") - # FIXME: this fails when order = 2 mesh = generate_gmsh(script, 3, order=order, target_unit="MM") logger.info("END GEN") From 53e819ef88795087500c24c24c5763ded3a7b484 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 9 Nov 2020 13:24:55 -0600 Subject: [PATCH 18/51] use more elements in gmsh quad test --- test/test_meshmode.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index c212b011b..c573cdf67 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1298,7 +1298,7 @@ def test_quad_mesh_3d(mesh_name, order=3, visualize=False): SetFactory("OpenCASCADE"); Box(1) = {0, 0, 0, 1, 1, 1}; - Transfinite Line "*" = 1; + Transfinite Line "*" = 8; Transfinite Surface "*"; Transfinite Volume "*"; From 1563970adf5ccecdac9d8b8a131554e8d8e2ff9f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 9 Nov 2020 15:07:23 -0600 Subject: [PATCH 19/51] fix typo in element group --- meshmode/discretization/poly_element.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 1d37bcbc6..6f17ef7be 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -424,11 +424,11 @@ class GaussLegendreTensorProductElementGroup(_LegendreTensorProductElementGroup) """ def __init__(self, mesh_el_group, order, index): import modepy as mp - quad = mp.LegendreGaussQuadrature(self.order) + quad = mp.LegendreGaussQuadrature(order) super().__init__(mesh_el_group, order, index, unit_nodes_1d=quad.nodes, - quad_nodes_1d=quad.weights) + quad_weights_1d=quad.weights) class LegendreGaussLobattoTensorProductElementGroup( From d6467f414854641be913f3908ddfdc9825dc97ac Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 12 Nov 2020 21:10:17 -0600 Subject: [PATCH 20/51] use mass matrix is quadrature weights are not provided --- meshmode/discretization/poly_element.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 6f17ef7be..36fb6378d 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -346,7 +346,6 @@ def unit_nodes(self): return self._unit_nodes - # }}} @@ -392,8 +391,14 @@ def from_mesh_interp_matrix(self): @property def weights(self): + if self._quad_weights_1d is None: + import modepy as mp + mm = mp.mass_matrix(self._basis_1d(), self._unit_nodes_1d) + weights = mm @ np.ones(len(self._unit_nodes_1d)) + else: + weights = self._quad_weights_1d + from itertools import product - weights = self._quad_weights_1d return np.fromiter( (np.prod(w) for w in product(weights, repeat=self.dim)), dtype=np.float, From 72b2e5199e594d03588f139a7e15cce3f00b52f5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 12 Nov 2020 21:57:51 -0600 Subject: [PATCH 21/51] use same order for the mesh and discr --- meshmode/discretization/connection/face.py | 6 ++---- test/test_meshmode.py | 4 ++-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 7cb7a33ea..5be724b04 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -337,9 +337,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # Code does not.) face_offset = face_vertex_unit_coordinates[0] - face_basis = ( - face_vertex_unit_coordinates[1:mgrp.dim] - - face_offset) + face_basis = face_vertex_unit_coordinates[1:] - face_offset if isinstance(mgrp, SimplexElementGroup): # no further action required @@ -353,7 +351,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, elif mgrp.dim == 3: # FIXME: This makes the tests pass, but I do not understand # why this reversal is necessary. - face_basis = face_basis[::-1] + face_basis = face_basis[-2::-1] # assert that all basis vectors are axis-aligned for bvec in face_basis: diff --git a/test/test_meshmode.py b/test/test_meshmode.py index c573cdf67..e6e59c0ef 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -440,7 +440,7 @@ def f(x): force_ambient_dim=2) elif mesh_name == "warp": from meshmode.mesh.generation import generate_warped_rect_mesh - mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par, + mesh = generate_warped_rect_mesh(dim, order=order, n=mesh_par, group_cls=group_cls) h = 1/mesh_par @@ -448,7 +448,7 @@ def f(x): elif mesh_name == "rect": from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh(a=(0,)*dim, b=(1,)*dim, - order=4, n=(mesh_par,)*dim, group_cls=group_cls) + order=order, n=(mesh_par,)*dim, group_cls=group_cls) h = 1/mesh_par else: From 90ffb1e3181a1578f3d5df624ac0fb4f8f08de78 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 12 Nov 2020 21:59:35 -0600 Subject: [PATCH 22/51] point back at gmsh_interop master --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 2f1803a26..7027fed52 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ numpy recursivenodes git+https://github.com/inducer/pytools.git#egg=pytools -git+https://github.com/alexfikl/gmsh_interop.git@contrib-node-order#egg=gmsh_interop +git+https://github.com/inducer/gmsh_interop.git#egg=gmsh_interop git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile git+https://github.com/inducer/modepy.git#egg=modepy git+https://github.com/inducer/pyopencl.git#egg=pyopencl From 73bea32d0abc0897dd287b78585741c4b42e4eb0 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 12 Nov 2020 22:04:19 -0600 Subject: [PATCH 23/51] memoize tensor product weights --- meshmode/discretization/poly_element.py | 1 + 1 file changed, 1 insertion(+) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 36fb6378d..80f0ca64b 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -390,6 +390,7 @@ def from_mesh_interp_matrix(self): meg.unit_nodes) @property + @memoize_method def weights(self): if self._quad_weights_1d is None: import modepy as mp From 5f27b785d0fe732a0ceef3707119d60ba7403d3e Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 18 Nov 2020 09:39:58 -0600 Subject: [PATCH 24/51] fix typo in tensor product mass matrix --- meshmode/discretization/poly_element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 80f0ca64b..538ff1985 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -394,7 +394,7 @@ def from_mesh_interp_matrix(self): def weights(self): if self._quad_weights_1d is None: import modepy as mp - mm = mp.mass_matrix(self._basis_1d(), self._unit_nodes_1d) + mm = mp.mass_matrix(self._basis_1d, self._unit_nodes_1d) weights = mm @ np.ones(len(self._unit_nodes_1d)) else: weights = self._quad_weights_1d From 6ad559bd94ce2f7394858b9d40895d9c377af49a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 19 Nov 2020 00:25:49 -0600 Subject: [PATCH 25/51] Explain the tensor product face restriction FIXME --- meshmode/discretization/connection/face.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 5be724b04..3312785d9 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -349,9 +349,24 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, pass elif mgrp.dim == 3: - # FIXME: This makes the tests pass, but I do not understand - # why this reversal is necessary. - face_basis = face_basis[-2::-1] + # The 3D vertices that generate the face basis for the + # volume cube are [a, b, c]. + + # ^ [1, 0] + # | + # b----c + # | | + # | | + # 0----a--> [0,1] + + # First, drop 'c': He's linearly dependent on a, b. + face_basis = face_basis[:-1] + + # Somehow, we've ended up with a 'last dimension varies + # fastest' convention for vertex numbering, which means + # that the multiplication by the face basis will swap the + # axes. Not good, undo: + face_basis = face_basis[::-1] # assert that all basis vectors are axis-aligned for bvec in face_basis: From 82be9e81f29a20fcbfe2455e20f6ccaa172737fa Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 1 Dec 2020 01:01:58 -0600 Subject: [PATCH 26/51] Make use of some shapes, spaces facilities from modepy in meshmode.mesh --- meshmode/discretization/connection/face.py | 6 - meshmode/mesh/__init__.py | 165 ++++----------------- meshmode/mesh/generation.py | 8 +- requirements.txt | 2 +- 4 files changed, 36 insertions(+), 145 deletions(-) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 3312785d9..86a2b2906 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -362,12 +362,6 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # First, drop 'c': He's linearly dependent on a, b. face_basis = face_basis[:-1] - # Somehow, we've ended up with a 'last dimension varies - # fastest' convention for vertex numbering, which means - # that the multiplication by the face basis will swap the - # axes. Not good, undo: - face_basis = face_basis[::-1] - # assert that all basis vectors are axis-aligned for bvec in face_basis: assert sum(abs(entry) > 1e-13 for entry in bvec) == 1 diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 87d285051..9f589337e 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -290,13 +290,9 @@ def __ne__(self, other): # }}} -# {{{ simplex - -class SimplexElementGroup(MeshElementGroup): - """ - .. automethod:: __init__ - """ +# {{{ modepy-based element group +class _ModepyElementGroup(MeshElementGroup): def __init__(self, order, vertex_indices, nodes, element_nr_base=None, node_nr_base=None, unit_nodes=None, dim=None): @@ -317,156 +313,57 @@ def __init__(self, order, vertex_indices, nodes, automatically assigned. """ - if unit_nodes is None: + if unit_nodes is not None: + _dim = unit_nodes.shape[0] + if dim is not None and _dim != dim: + raise ValueError("dim does not match unit_nodes") + else: + dim = _dim + else: if dim is None: raise TypeError("'dim' must be passed " "if 'unit_nodes' is not passed") - if dim <= 3: - unit_nodes = mp.warp_and_blend_nodes(dim, order) - else: - unit_nodes = mp.equidistant_nodes(dim, order) + # dim is now usable + shape = self._modepy_shape = self._modepy_shape_cls(dim) + space = self._modepy_space = mp.space_for_shape(self._modepy_shape, order) - dims = unit_nodes.shape[0] + if unit_nodes is None: + unit_nodes = mp.edge_clustered_nodes_for_space(space, shape) if vertex_indices is not None: if not issubclass(vertex_indices.dtype.type, np.integer): raise TypeError("vertex_indices must be integral") - if vertex_indices.shape[-1] != dims+1: + if vertex_indices.shape[-1] != shape.nvertices: raise ValueError("vertex_indices has wrong number of vertices per " - "element. expected: %d, got: %d" % (dims+1, + "element. expected: %d, got: %d" % (shape.nvertices, vertex_indices.shape[-1])) super().__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 ( - (0,), - (1,), - ) - elif self.dim == 2: - return ( - (0, 1), - (2, 0), - (1, 2), - ) - elif self.dim == 3: - return ( - (0, 2, 1), - (0, 1, 3), - (0, 3, 2), - (1, 2, 3) - ) - else: - raise NotImplementedError(f"dimension {self.dim}") + return tuple(face.volume_vertex_indices + for face in mp.faces_for_shape(self._modepy_shape)) def vertex_unit_coordinates(self): - from modepy.tools import unit_vertices - return unit_vertices(self.dim) + return mp.biunit_vertices_for_shape(self._modepy_shape).T # }}} -# {{{ tensor-product - -class TensorProductElementGroup(MeshElementGroup): - """ - .. automethod:: __init__ - """ - - def __init__(self, order, vertex_indices, nodes, - element_nr_base=None, node_nr_base=None, - unit_nodes=None): - r""" - :arg order: the maximum total degree used for interpolation. - :arg nodes: ``[ambient_dim, nelements, nunit_nodes]`` - The nodes are assumed to be mapped versions of *unit_nodes*. - :arg unit_nodes: ``[dim, nunit_nodes]`` - The unit nodes of which *nodes* is a mapped - version. - - The vertex order in 3D is as follows :: - +class SimplexElementGroup(_ModepyElementGroup): + _modepy_shape_cls = mp.Simplex - ^ s - | 6-----------7 - | /| /| - | / | ^ / | - |/ | / t / | - 2-----------3 | - | |/ | | - | 4-------|---5 - | / | / - | / | / - |/ |/ - 0-----------1---->r - - For general dimensions, follow binary upward counting: - 000, 001, 010, 011, ... - - Do not supply *element_nr_base* and *node_nr_base*, they will be - automatically assigned. - """ - - dims = unit_nodes.shape[0] - - if vertex_indices is not None: - if not issubclass(vertex_indices.dtype.type, np.integer): - raise TypeError("vertex_indices must be integral") - - if vertex_indices.shape[-1] != 2**dims: - raise ValueError("vertex_indices has wrong number of vertices per " - "element. expected: %d, got: %d" % (2**dims, - vertex_indices.shape[-1])) - - super().__init__(order, vertex_indices, nodes, - element_nr_base, node_nr_base, unit_nodes) - - def face_vertex_indices(self): - if self.dim == 1: - return ( - (0,), - (1,), - ) - elif self.dim == 2: - return ( - (0, 1), - (3, 2), - (2, 0), - (1, 3), - ) - elif self.dim == 3: - return ( - # binary is tsr - - # rs-aligned - (0b000, 0b001, 0b010, 0b011,), - (0b100, 0b101, 0b110, 0b111,), - - # st-aligned - (0b000, 0b010, 0b100, 0b110,), - (0b001, 0b011, 0b101, 0b111,), - - # rt-aligned - (0b000, 0b001, 0b100, 0b101,), - (0b010, 0b011, 0b110, 0b111,), - ) - else: - raise NotImplementedError(f"dimension {self.dim}") + @property + @memoize_method + def is_affine(self): + return is_affine_simplex_group(self) - def vertex_unit_coordinates(self): - from pytools import generate_nonnegative_integer_tuples_below as gnitb - return np.array(list(gnitb(2, self.dim)), dtype=np.float64)*2.0 - 1.0 -# }}} +class TensorProductElementGroup(_ModepyElementGroup): + _modepy_shape_cls = mp.Hypercube # }}} @@ -1029,10 +926,8 @@ def _test_node_vertex_consistency_resampling(mesh, mgrp, tol): if mgrp.nelements == 0: return True - if isinstance(mgrp, SimplexElementGroup): - basis = mp.simplex_best_available_basis(mgrp.dim, mgrp.order) - elif isinstance(mgrp, TensorProductElementGroup): - basis = mp.legendre_tensor_product_basis(mgrp.dim, mgrp.order) + if isinstance(mgrp, _ModepyElementGroup): + basis = mp.basis_for_space(mgrp._modepy_space).functions else: raise TypeError(f"unsupported group type: {type(mgrp).__name__}") @@ -1071,7 +966,7 @@ def _test_node_vertex_consistency(mesh, tol): """ for mgrp in mesh.groups: - if isinstance(mgrp, (SimplexElementGroup, TensorProductElementGroup)): + if isinstance(mgrp, _ModepyElementGroup): assert _test_node_vertex_consistency_resampling(mesh, mgrp, tol) else: from warnings import warn diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index dbf16fd6c..15bad928d 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -362,7 +362,9 @@ def make_group_from_vertices(vertices, vertex_indices, order, _, nnodes = unit_nodes.shape from pytools import generate_nonnegative_integer_tuples_below as gnitb - id_tuples = list(gnitb(2, dim)) + # Modepy's convention is 'first coordinate varies fastest'. Adapt to it + # for the vertices. It doesn't matter for the basis. + id_tuples = [id_tuple[::-1] for id_tuple in gnitb(2, dim)] assert len(id_tuples) == nvertices vdm = np.empty((nvertices, nvertices)) @@ -849,8 +851,8 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64, if is_tp: el_vertices.append( - (a000, a001, a010, a011, - a100, a101, a110, a111)) + (a000, a100, a010, a110, + a001, a101, a011, a111)) else: el_vertices.append((a000, a100, a010, a001)) diff --git a/requirements.txt b/requirements.txt index 7027fed52..e5d58fdde 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,7 +3,7 @@ recursivenodes git+https://github.com/inducer/pytools.git#egg=pytools git+https://github.com/inducer/gmsh_interop.git#egg=gmsh_interop git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile -git+https://github.com/inducer/modepy.git#egg=modepy +git+https://github.com/inducer/modepy.git@refactor-shapes#egg=modepy git+https://github.com/inducer/pyopencl.git#egg=pyopencl git+https://github.com/inducer/islpy.git#egg=islpy From d79952d54460f94ae079bc50a123ad3862fed0c9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 19:14:38 -0600 Subject: [PATCH 27/51] add more shape checks to MeshElementGroup --- meshmode/mesh/__init__.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 9f589337e..405431dac 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -316,13 +316,12 @@ def __init__(self, order, vertex_indices, nodes, if unit_nodes is not None: _dim = unit_nodes.shape[0] if dim is not None and _dim != dim: - raise ValueError("dim does not match unit_nodes") + raise ValueError("'dim' does not match 'unit_nodes' dimension") else: dim = _dim else: if dim is None: - raise TypeError("'dim' must be passed " - "if 'unit_nodes' is not passed") + raise TypeError("'dim' must be passed if 'unit_nodes' is not passed") # dim is now usable shape = self._modepy_shape = self._modepy_shape_cls(dim) @@ -331,14 +330,21 @@ def __init__(self, order, vertex_indices, nodes, if unit_nodes is None: unit_nodes = mp.edge_clustered_nodes_for_space(space, shape) + if unit_nodes.shape[-1] != nodes.shape[-1]: + raise ValueError( + "'nodes' has wrong number of unit nodes per element." + f" expected {unit_nodes.shape[-1]}, " + f" but got {nodes.shape[-1]}.") + if vertex_indices is not None: if not issubclass(vertex_indices.dtype.type, np.integer): - raise TypeError("vertex_indices must be integral") + raise TypeError("'vertex_indices' must be integral") if vertex_indices.shape[-1] != shape.nvertices: - raise ValueError("vertex_indices has wrong number of vertices per " - "element. expected: %d, got: %d" % (shape.nvertices, - vertex_indices.shape[-1])) + raise ValueError( + "'vertex_indices' has wrong number of vertices per element." + f" expected {shape.nvertices}," + f" got {vertex_indices.shape[-1]}") super().__init__(order, vertex_indices, nodes, element_nr_base, node_nr_base, unit_nodes, dim) From dae30bb0cc7340ff9d4012739cdd9f1131e74359 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 19:56:48 -0600 Subject: [PATCH 28/51] expand test_sanity_single_element for tensor products --- test/test_meshmode.py | 62 ++++++++++++++++++++++++++----------------- 1 file changed, 38 insertions(+), 24 deletions(-) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index a91ff9505..cc3080dfe 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -873,43 +873,57 @@ def test_merge_and_map(actx_factory, group_cls, visualize=False): # {{{ sanity checks: single element @pytest.mark.parametrize("dim", [2, 3]) -@pytest.mark.parametrize("order", [1, 3]) -def test_sanity_single_element(actx_factory, dim, order, visualize=False): +@pytest.mark.parametrize("mesh_order", [1, 3]) +@pytest.mark.parametrize("group_cls", [ + SimplexElementGroup, + TensorProductElementGroup, + ]) +def test_sanity_single_element(actx_factory, dim, mesh_order, group_cls, + visualize=False): pytest.importorskip("pytential") actx = actx_factory() - from modepy.tools import unit_vertices - vertices = unit_vertices(dim).T.copy() + if group_cls is SimplexElementGroup: + from meshmode.discretization.poly_element import \ + PolynomialWarpAndBlendGroupFactory as GroupFactory + elif group_cls is TensorProductElementGroup: + from meshmode.discretization.poly_element import \ + LegendreGaussLobattoTensorProductGroupFactory as GroupFactory + else: + raise TypeError + + import modepy as mp + shape = group_cls._modepy_shape_cls(dim) + space = mp.space_for_shape(shape, mesh_order) + + vertices = mp.biunit_vertices_for_shape(shape) + nodes = mp.edge_clustered_nodes_for_space(space, shape).reshape(dim, 1, -1) + vertex_indices = np.arange(shape.nvertices, dtype=np.int32).reshape(1, -1) center = np.empty(dim, np.float64) center.fill(-0.5) - import modepy as mp - from meshmode.mesh import SimplexElementGroup, Mesh, BTAG_ALL - mg = SimplexElementGroup( - order=order, - vertex_indices=np.arange(dim+1, dtype=np.int32).reshape(1, -1), - nodes=mp.warp_and_blend_nodes(dim, order).reshape(dim, 1, -1), - dim=dim) - + from meshmode.mesh import Mesh, BTAG_ALL + mg = group_cls(mesh_order, vertex_indices, nodes, dim=dim) mesh = Mesh(vertices, [mg], is_conforming=True) from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory - vol_discr = Discretization(actx, mesh, - PolynomialWarpAndBlendGroupFactory(order+3)) + vol_discr = Discretization(actx, mesh, GroupFactory(mesh_order + 3)) # {{{ volume calculation check - vol_x = thaw(actx, vol_discr.nodes()) - - vol_one = vol_x[0] * 0 + 1 - from pytential import norm, integral # noqa + if isinstance(mg, SimplexElementGroup): + from pytools import factorial + true_vol = 1/factorial(dim) * 2**dim + elif isinstance(mg, TensorProductElementGroup): + true_vol = 2**dim + else: + raise TypeError - from pytools import factorial - true_vol = 1/factorial(dim) * 2**dim + nodes = thaw(actx, vol_discr.nodes()) + vol_one = 1 + 0 * nodes[0] + from pytential import norm, integral # noqa comp_vol = integral(vol_discr, vol_one) rel_vol_err = abs(true_vol - comp_vol) / true_vol @@ -921,14 +935,14 @@ def test_sanity_single_element(actx_factory, dim, order, visualize=False): from meshmode.discretization.connection import make_face_restriction bdry_connection = make_face_restriction( - actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3), + actx, vol_discr, GroupFactory(mesh_order + 3), BTAG_ALL) bdry_discr = bdry_connection.to_discr # }}} from pytential import bind, sym - bdry_normals = bind(bdry_discr, sym.normal(dim))(actx).as_vector(dtype=object) + bdry_normals = bind(bdry_discr, sym.normal(dim).as_vector())(actx) if visualize: from meshmode.discretization.visualization import make_visualizer From 8cd7a3857cf6742b1daf337d67a7e1dc4a1d4d8e Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 20:48:37 -0600 Subject: [PATCH 29/51] use modepy unit nodes in gmsh reader --- meshmode/mesh/io.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 7c345b04c..e8a24d941 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -197,8 +197,17 @@ def get_mesh(self): i += 1 - unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(), - dtype=np.float64).T/group_el_type.order)*2 - 1 + import modepy as mp + if isinstance(group_el_type, GmshSimplexElementBase): + shape = mp.Simplex(group_el_type.dimensions) + elif isinstance(group_el_type, GmshTensorProductElementBase): + shape = mp.Hypercube(group_el_type.dimensions) + else: + raise NotImplementedError( + f"gmsh element type: {type(group_el_type).__name__}") + + space = mp.space_for_shape(shape, group_el_type.order) + unit_nodes = mp.equispaced_nodes_for_space(space, shape) if isinstance(group_el_type, GmshSimplexElementBase): group = SimplexElementGroup( @@ -224,8 +233,7 @@ def get_mesh(self): unit_nodes=unit_nodes ) else: - raise NotImplementedError("gmsh element type: %s" - % type(group_el_type).__name__) + pass groups.append(group) From 5a06cbf8178ce7d1151d7276372771e16bce4061 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 20:49:14 -0600 Subject: [PATCH 30/51] point requirements.txt to modepy branch --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index e5d58fdde..c65ea9f71 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,7 +3,7 @@ recursivenodes git+https://github.com/inducer/pytools.git#egg=pytools git+https://github.com/inducer/gmsh_interop.git#egg=gmsh_interop git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile -git+https://github.com/inducer/modepy.git@refactor-shapes#egg=modepy +git+https://github.com/alexfikl/modepy.git@tensor-face-mass-matrix#egg=modepy git+https://github.com/inducer/pyopencl.git#egg=pyopencl git+https://github.com/inducer/islpy.git#egg=islpy From bb5080efd5bfe1afe45dda81b145ed915c9250ad Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 21:11:41 -0600 Subject: [PATCH 31/51] update simple-dg example --- examples/simple-dg.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/simple-dg.py b/examples/simple-dg.py index 324f2a925..60a747ebe 100644 --- a/examples/simple-dg.py +++ b/examples/simple-dg.py @@ -273,12 +273,13 @@ def get_local_face_mass_matrix(self, afgrp, volgrp, dtype): afgrp.nunit_dofs), dtype=dtype) - from modepy.tools import UNIT_VERTICES import modepy as mp - for iface, fvi in enumerate( - volgrp.mesh_el_group.face_vertex_indices()): - face_vertices = UNIT_VERTICES[volgrp.dim][np.array(fvi)].T - matrix[:, iface, :] = mp.nodal_face_mass_matrix( + shape = mp.Simplex(volgrp.dim) + unit_vertices = mp.biunit_vertices_for_shape(shape).T + + for face in mp.faces_for_shape(shape): + face_vertices = unit_vertices[np.array(face.volume_vertex_indices)].T + matrix[:, face.face_index, :] = mp.nodal_face_mass_matrix( volgrp.basis(), volgrp.unit_nodes, afgrp.unit_nodes, volgrp.order, face_vertices) From 5912ccbbdba32d69989580dc8b4359a8205ec5ca Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 21:23:32 -0600 Subject: [PATCH 32/51] allow missing nodes for firedrake --- meshmode/mesh/__init__.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 405431dac..f0dbf5826 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -330,11 +330,12 @@ def __init__(self, order, vertex_indices, nodes, if unit_nodes is None: unit_nodes = mp.edge_clustered_nodes_for_space(space, shape) - if unit_nodes.shape[-1] != nodes.shape[-1]: - raise ValueError( - "'nodes' has wrong number of unit nodes per element." - f" expected {unit_nodes.shape[-1]}, " - f" but got {nodes.shape[-1]}.") + if nodes is not None: + if unit_nodes.shape[-1] != nodes.shape[-1]: + raise ValueError( + "'nodes' has wrong number of unit nodes per element." + f" expected {unit_nodes.shape[-1]}, " + f" but got {nodes.shape[-1]}.") if vertex_indices is not None: if not issubclass(vertex_indices.dtype.type, np.integer): From 7d0d861bb9cff9f854c72ba47e7de50679e71e05 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 21:47:06 -0600 Subject: [PATCH 33/51] implement get_copy_kwargs --- meshmode/mesh/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index f0dbf5826..b0cb0168d 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -219,12 +219,13 @@ def __init__(self, order, vertex_indices, nodes, unit_nodes=unit_nodes, element_nr_base=element_nr_base, node_nr_base=node_nr_base) - def copy(self, **kwargs): + def get_copy_kwargs(self, **kwargs): if "element_nr_base" not in kwargs: kwargs["element_nr_base"] = None if "node_nr_base" not in kwargs: kwargs["node_nr_base"] = None - return Record.copy(self, **kwargs) + + return super().get_copy_kwargs(**kwargs) @property def dim(self): From d53887dec19630bcdfadc900e379d0615008dfd5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 1 Dec 2020 22:27:39 -0600 Subject: [PATCH 34/51] initialize MeshElementGroup fully on unpickling --- meshmode/mesh/__init__.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index b0cb0168d..3607d2bec 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -351,6 +351,13 @@ def __init__(self, order, vertex_indices, nodes, super().__init__(order, vertex_indices, nodes, element_nr_base, node_nr_base, unit_nodes, dim) + def __setstate__(self, valuedict): + super().__setstate__(valuedict) + + import modepy as mp + self._modepy_shape = self._modepy_shape_cls(self.dim) + self._modepy_space = mp.space_for_shape(self._modepy_shape, self.order) + def face_vertex_indices(self): return tuple(face.volume_vertex_indices for face in mp.faces_for_shape(self._modepy_shape)) From 6fa04a8e8b6924da76382e1f6947ee37a2233543 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Dec 2020 11:03:51 -0600 Subject: [PATCH 35/51] update basis_for_space calls --- meshmode/mesh/__init__.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 3607d2bec..ea349cdf4 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -942,7 +942,7 @@ def _test_node_vertex_consistency_resampling(mesh, mgrp, tol): return True if isinstance(mgrp, _ModepyElementGroup): - basis = mp.basis_for_space(mgrp._modepy_space).functions + basis = mp.basis_for_space(mgrp._modepy_space, mgrp._modepy_shape).functions else: raise TypeError(f"unsupported group type: {type(mgrp).__name__}") @@ -1403,11 +1403,10 @@ def is_affine_simplex_group(group, abs_tol=None): 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) + basis = mp.basis_for_space(group._modepy_space, group._modepy_shape) + vinv = la.inv(mp.vandermonde(basis.functions, group.unit_nodes)) + diff = mp.differentiation_matrices( + basis.functions, basis.gradients, group.unit_nodes) if not isinstance(diff, tuple): diff = (diff,) From 0520f71028f5b95d439398d61ab778ad234af86a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Dec 2020 11:11:16 -0600 Subject: [PATCH 36/51] rename biunit_vertices_for_space after modepy --- examples/simple-dg.py | 2 +- meshmode/discretization/__init__.py | 8 ++++---- meshmode/interop/nodal_dg.py | 4 ++-- meshmode/mesh/__init__.py | 2 +- test/test_meshmode.py | 2 +- 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/simple-dg.py b/examples/simple-dg.py index 60a747ebe..ac109b6be 100644 --- a/examples/simple-dg.py +++ b/examples/simple-dg.py @@ -275,7 +275,7 @@ def get_local_face_mass_matrix(self, afgrp, volgrp, dtype): import modepy as mp shape = mp.Simplex(volgrp.dim) - unit_vertices = mp.biunit_vertices_for_shape(shape).T + unit_vertices = mp.unit_vertices_for_shape(shape).T for face in mp.faces_for_shape(shape): face_vertices = unit_vertices[np.array(face.volume_vertex_indices)].T diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py index f164d5d72..7ec816b23 100644 --- a/meshmode/discretization/__init__.py +++ b/meshmode/discretization/__init__.py @@ -91,15 +91,15 @@ def nelements(self): @property def nunit_dofs(self): - """The number of (for now: nodal) degrees of freedom ("DOFs") associated with a - single element. + """The number of (for now: nodal) degrees of freedom ("DOFs") + associated with a single element. """ return self.unit_nodes.shape[-1] @property def ndofs(self): - """The total number of (for now: nodal) degrees of freedom ("DOFs") associated with - the element group. + """The total number of (for now: nodal) degrees of freedom ("DOFs") + associated with the element group. """ return self.nunit_dofs * self.nelements diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py index 98fc278e4..8fe42fd02 100644 --- a/meshmode/interop/nodal_dg.py +++ b/meshmode/interop/nodal_dg.py @@ -113,11 +113,11 @@ def get_discr(self, actx) -> meshmode.discretization.Discretization: unit_nodes_arrays = self.octave.eval( f"Nodes{dim}D(N)", nout=dim, verbose=False) - equilat_to_biunit_func_name = ( + equilat_to_unit_func_name = ( "".join(self.AXES[:dim] + ["to"] + self.REF_AXES[:dim])) unit_nodes_arrays = self.octave.feval( - equilat_to_biunit_func_name, *unit_nodes_arrays, + equilat_to_unit_func_name, *unit_nodes_arrays, nout=dim, verbose=False) unit_nodes = np.array([a.reshape(-1) for a in unit_nodes_arrays]) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index ea349cdf4..1dd7f419f 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -363,7 +363,7 @@ def face_vertex_indices(self): for face in mp.faces_for_shape(self._modepy_shape)) def vertex_unit_coordinates(self): - return mp.biunit_vertices_for_shape(self._modepy_shape).T + return mp.unit_vertices_for_shape(self._modepy_shape).T # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cc3080dfe..ea44e3e83 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -896,7 +896,7 @@ def test_sanity_single_element(actx_factory, dim, mesh_order, group_cls, shape = group_cls._modepy_shape_cls(dim) space = mp.space_for_shape(shape, mesh_order) - vertices = mp.biunit_vertices_for_shape(shape) + vertices = mp.unit_vertices_for_shape(shape) nodes = mp.edge_clustered_nodes_for_space(space, shape).reshape(dim, 1, -1) vertex_indices = np.arange(shape.nvertices, dtype=np.int32).reshape(1, -1) From d775213a9c391a1cda8c9e26e77af4362c396821 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Dec 2020 14:20:00 -0600 Subject: [PATCH 37/51] use more modern modepy in poly_element --- meshmode/discretization/poly_element.py | 130 +++++++++++++----------- 1 file changed, 71 insertions(+), 59 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 538ff1985..18a1f8804 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -104,11 +104,22 @@ def diff_matrices(self): # {{{ concrete element groups for simplices class SimplexElementGroupBase(ElementGroupBase): + @property + @memoize_method + def _shape(self): + return mp.Simplex(self.dim) + + @property + @memoize_method + def _space(self): + return mp.PN(self.dim, self.order) + @memoize_method def from_mesh_interp_matrix(self): meg = self.mesh_el_group + meg_space = mp.PN(meg.dim, meg.order) return mp.resampling_matrix( - mp.simplex_best_available_basis(meg.dim, meg.order), + mp.basis_for_space(meg_space, self._shape).functions, self.unit_nodes, meg.unit_nodes) @@ -118,28 +129,19 @@ class PolynomialSimplexElementGroupBase(PolynomialElementGroupBase, def is_orthogonal_basis(self): return self.dim <= 3 + @property @memoize_method - def _mode_ids_and_basis(self): - # for now, see https://gitlab.tiker.net/inducer/modepy/-/merge_requests/14 - import modepy.modes as modes - if self.dim <= 3: - return modes.simplex_onb_with_mode_ids(self.dim, self.order) - else: - return modes.simplex_monomial_basis_with_mode_ids(self.dim, self.order) + def _basis(self): + return mp.basis_for_space(self._space, self._shape) def basis(self): - mode_ids, basis = self._mode_ids_and_basis() - return basis + return self._basis.functions def mode_ids(self): - mode_ids, basis = self._mode_ids_and_basis() - return mode_ids + return self._basis.mode_ids def grad_basis(self): - if self.dim <= 3: - return mp.grad_simplex_onb(self.dim, self.order) - else: - return mp.grad_simplex_monomial_basis(self.dim, self.order) + return self._basis.gradients class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): @@ -351,77 +353,88 @@ def unit_nodes(self): # {{{ concrete element groups for tensor product elements -class _TensorProductElementGroupBase(PolynomialElementGroupBase): +class HypercubeElementGroupBase(ElementGroupBase): + @property + @memoize_method + def _shape(self): + return mp.Hypercube(self.dim) + + @property + @memoize_method + def _space(self): + return mp.QN(self.dim, self.order) + + @memoize_method + def from_mesh_interp_matrix(self): + meg = self.mesh_el_group + meg_space = mp.QN(meg.dim, meg.order) + return mp.resampling_matrix( + mp.basis_for_space(meg_space, self._shape).functions, + self.unit_nodes, + meg.unit_nodes) + + +class TensorProductElementGroupBase(PolynomialElementGroupBase, + HypercubeElementGroupBase): def __init__(self, mesh_el_group, order, index, *, - basis_1d=None, grad_basis_1d=None, - unit_nodes_1d=None, quad_weights_1d=None): + basis, unit_nodes_1d, quad_weights_1d=None): super().__init__(mesh_el_group, order, index) - self._basis_1d = basis_1d - self._grad_basis_1d = grad_basis_1d + self._basis = basis self._unit_nodes_1d = unit_nodes_1d self._quad_weights_1d = quad_weights_1d def is_orthogonal_basis(self): - return True + try: + self._basis.orthonormality_weight() + return True + except mp.BasisNotOrthonormal: + return False + + def mode_ids(self): + return self._basis.mode_ids def basis(self): - from modepy.modes import tensor_product_basis - return tensor_product_basis(self.dim, self._basis_1d) + return self._basis.functions def grad_basis(self): - from modepy.modes import grad_tensor_product_basis - return grad_tensor_product_basis(self.dim, - self._basis_1d, self._grad_basis_1d) + return self._basis.gradients @property @memoize_method def unit_nodes(self): - from modepy.nodes import tensor_product_nodes - return tensor_product_nodes(self.dim, self._unit_nodes_1d) - - @memoize_method - def from_mesh_interp_matrix(self): - from modepy.modes import legendre_tensor_product_basis - meg = self.mesh_el_group - return mp.resampling_matrix( - legendre_tensor_product_basis(self.dim, meg.order), - self.unit_nodes, - meg.unit_nodes) + # FIXME: allow non-homogeneous tensor products + return mp.tensor_product_nodes(self.dim, self._unit_nodes_1d) @property @memoize_method def weights(self): if self._quad_weights_1d is None: - import modepy as mp - mm = mp.mass_matrix(self._basis_1d, self._unit_nodes_1d) - weights = mm @ np.ones(len(self._unit_nodes_1d)) + return np.dot( + self.mass_matrix(), + np.ones(len(self.basis()))) else: - weights = self._quad_weights_1d - - from itertools import product - return np.fromiter( - (np.prod(w) for w in product(weights, repeat=self.dim)), - dtype=np.float, - count=(self.order + 1)**self.dim) + # FIXME: allow non-homogeneous tensor products + return np.prod( + mp.tensor_product_nodes(self.dim, self._quad_weights_1d), + axis=0) -class _LegendreTensorProductElementGroup(_TensorProductElementGroupBase): +class LegendreTensorProductElementGroup(TensorProductElementGroupBase): def __init__(self, mesh_el_group, order, index, *, unit_nodes_1d=None, quad_weights_1d=None): - from modepy.modes import jacobi, grad_jacobi - from functools import partial + + basis = mp.orthonormal_basis_for_space( + mp.QN(mesh_el_group.dim, order), + mp.Hypercube(mesh_el_group.dim)) super().__init__(mesh_el_group, order, index, - basis_1d=tuple( - partial(jacobi, 0, 0, i) for i in range(order + 1)), - grad_basis_1d=tuple( - partial(grad_jacobi, 0, 0, i) for i in range(order + 1)), + basis=basis, unit_nodes_1d=unit_nodes_1d, quad_weights_1d=quad_weights_1d) -class GaussLegendreTensorProductElementGroup(_LegendreTensorProductElementGroup): +class GaussLegendreTensorProductElementGroup(LegendreTensorProductElementGroup): """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. @@ -429,7 +442,6 @@ class GaussLegendreTensorProductElementGroup(_LegendreTensorProductElementGroup) No interpolation nodes are present on the boundary of the hypercube. """ def __init__(self, mesh_el_group, order, index): - import modepy as mp quad = mp.LegendreGaussQuadrature(order) super().__init__(mesh_el_group, order, index, @@ -438,7 +450,7 @@ def __init__(self, mesh_el_group, order, index): class LegendreGaussLobattoTensorProductElementGroup( - _LegendreTensorProductElementGroup): + LegendreTensorProductElementGroup): """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. @@ -453,7 +465,7 @@ def __init__(self, mesh_el_group, order, index): unit_nodes_1d=legendre_gauss_lobatto_nodes(order)) -class EquidistantTensorProductElementGroup(_LegendreTensorProductElementGroup): +class EquidistantTensorProductElementGroup(LegendreTensorProductElementGroup): """Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. From 760d92b205d7cb8e881fb3e8b760be44887252e9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Dec 2020 14:43:19 -0600 Subject: [PATCH 38/51] update modern modepy in face restriction --- meshmode/discretization/connection/face.py | 120 ++++----------------- meshmode/mesh/__init__.py | 8 +- test/test_meshmode.py | 2 +- 3 files changed, 29 insertions(+), 101 deletions(-) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 86a2b2906..b0c855531 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -74,10 +74,7 @@ def _build_boundary_connection(actx, vol_discr, bdry_discr, connection_data, bdry_grp = bdry_discr.groups[ibdry_grp] data = connection_data[igrp, face_id] - bdry_unit_nodes_01 = (bdry_grp.unit_nodes + 1)*0.5 - result_unit_nodes = ( - np.dot(data.face_basis.T, bdry_unit_nodes_01).T - + data.face_offset).T + result_unit_nodes = data.face.map_to_volume(bdry_grp.unit_nodes) batches.append( InterpolationBatch( @@ -210,7 +207,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # }}} - from meshmode.mesh import Mesh, SimplexElementGroup, TensorProductElementGroup + from meshmode.mesh import Mesh, _ModepyElementGroup bdry_mesh_groups = [] connection_data = {} @@ -222,7 +219,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, mgrp = grp.mesh_el_group - if not isinstance(mgrp, (SimplexElementGroup, TensorProductElementGroup)): + if not isinstance(mgrp, _ModepyElementGroup): raise NotImplementedError("can only take boundary of " "meshes based on SimplexElementGroup and " "TensorProductElementGroup") @@ -260,61 +257,42 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # }}} - grp_face_vertex_indices = mgrp.face_vertex_indices() - grp_vertex_unit_coordinates = mgrp.vertex_unit_coordinates() - batch_base = 0 - # group by face_id + # group by face_index - for face_id in range(mgrp.nfaces): - batch_boundary_el_numbers_in_grp = np.array( - [ - ibface_el - for ibface_el, ibface_face in group_boundary_faces - if ibface_face == face_id], - dtype=np.intp) + for face in mgrp.faces: + batch_boundary_el_numbers_in_grp = np.array([ + ibface_el + for ibface_el, ibface_face in group_boundary_faces + if ibface_face == face.face_index + ], dtype=np.intp) # {{{ preallocate arrays for mesh group nbatch_elements = len(batch_boundary_el_numbers_in_grp) - if per_face_groups or face_id == 0: + if per_face_groups or face.face_index == 0: if per_face_groups: ngroup_bdry_elements = nbatch_elements else: ngroup_bdry_elements = len(group_boundary_faces) # make up some not-terrible nodes for the boundary Mesh - if isinstance(mgrp, SimplexElementGroup): - bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order) - vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order) - - nbdry_el_vertices = (mgrp.dim-1) + 1 + space = mp.space_for_shape(face, mgrp.order) + bdry_unit_nodes = mp.edge_clustered_nodes_for_space(space, face) - elif isinstance(mgrp, TensorProductElementGroup): - bdry_unit_nodes = \ - mp.legendre_gauss_lobatto_tensor_product_nodes( - mgrp.dim-1, mgrp.order) - vol_basis = mp.legendre_tensor_product_basis( - mgrp.dim, mgrp.order) - - nbdry_el_vertices = 2**(mgrp.dim-1) - - else: - assert False + vol_basis = mp.basis_for_space( + mgrp._modepy_space, mgrp._modepy_shape).functions vertex_indices = np.empty( - (ngroup_bdry_elements, nbdry_el_vertices), + (ngroup_bdry_elements, face.nvertices), mgrp.vertex_indices.dtype) - bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5 - - nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1] + nbdry_unit_nodes = bdry_unit_nodes.shape[-1] nodes = np.empty( (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes), dtype=np.float64) - # }}} new_el_numbers = batch_base + np.arange(nbatch_elements) @@ -323,60 +301,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # {{{ no per-element axes in these computations - # Find boundary vertex indices - loc_face_vertices = list(grp_face_vertex_indices[face_id]) - - # Find unit nodes for boundary element - face_vertex_unit_coordinates = \ - grp_vertex_unit_coordinates[loc_face_vertices] - - # Find face_basis, face_offset such that - # face_basis.T @ [e_1 e_2] + face_offset == [r_1 r_2]. - # - # (Notation assumes that the volume is 3D and the face is 2D. - # Code does not.) - - face_offset = face_vertex_unit_coordinates[0] - face_basis = face_vertex_unit_coordinates[1:] - face_offset - - if isinstance(mgrp, SimplexElementGroup): - # no further action required - pass - - elif isinstance(mgrp, TensorProductElementGroup): - if mgrp.dim <= 2: - # faces are simplices, no action required - pass - - elif mgrp.dim == 3: - # The 3D vertices that generate the face basis for the - # volume cube are [a, b, c]. - - # ^ [1, 0] - # | - # b----c - # | | - # | | - # 0----a--> [0,1] - - # First, drop 'c': He's linearly dependent on a, b. - face_basis = face_basis[:-1] - - # assert that all basis vectors are axis-aligned - for bvec in face_basis: - assert sum(abs(entry) > 1e-13 for entry in bvec) == 1 - - else: - raise NotImplementedError( - f"face restriction of a {mgrp.dim}D hypercube") - - else: - assert False - - face_unit_nodes = ( - np.dot(face_basis.T, bdry_unit_nodes_01).T - + face_offset).T - + face_unit_nodes = face.map_to_volume(bdry_unit_nodes) resampling_mat = mp.resampling_matrix( vol_basis, face_unit_nodes, mgrp.unit_nodes) @@ -387,7 +312,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # Find vertex_indices glob_face_vertices = mgrp.vertex_indices[ - batch_boundary_el_numbers_in_grp][:, loc_face_vertices] + batch_boundary_el_numbers_in_grp][:, face.volume_vertex_indices] vertex_indices[new_el_numbers] = \ vol_to_bdry_vertices[glob_face_vertices] @@ -399,14 +324,13 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # }}} - connection_data[igrp, face_id] = _ConnectionBatchData( + connection_data[igrp, face.face_index] = _ConnectionBatchData( group_source_element_indices=batch_boundary_el_numbers_in_grp, group_target_element_indices=new_el_numbers, - face_basis=face_basis, - face_offset=face_offset, + face=face, ) - is_last_face = face_id + 1 == mgrp.nfaces + is_last_face = face.face_index + 1 == mgrp.nfaces if per_face_groups or is_last_face: bdry_mesh_group = type(mgrp)( diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 1dd7f419f..d7dc8f281 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -358,9 +358,13 @@ def __setstate__(self, valuedict): self._modepy_shape = self._modepy_shape_cls(self.dim) self._modepy_space = mp.space_for_shape(self._modepy_shape, self.order) + @property + @memoize_method + def faces(self): + return mp.faces_for_shape(self._modepy_shape) + def face_vertex_indices(self): - return tuple(face.volume_vertex_indices - for face in mp.faces_for_shape(self._modepy_shape)) + return tuple(face.volume_vertex_indices for face in self.faces) def vertex_unit_coordinates(self): return mp.unit_vertices_for_shape(self._modepy_shape).T diff --git a/test/test_meshmode.py b/test/test_meshmode.py index ea44e3e83..888070337 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -718,7 +718,7 @@ def f(x): print(eoc_rec) assert ( eoc_rec.order_estimate() >= order-0.5 - or eoc_rec.max_error() < 1.5e-13) + or eoc_rec.max_error() < 1.7e-13) # }}} From 4ea8c458e701013fe0a56f0f0f1745dd68aad12c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 3 Dec 2020 10:48:21 -0600 Subject: [PATCH 39/51] fix MeshElementGroup subclass pickling --- meshmode/mesh/__init__.py | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index d7dc8f281..1f64dc238 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -199,7 +199,7 @@ class MeshElementGroup(Record): def __init__(self, order, vertex_indices, nodes, element_nr_base=None, node_nr_base=None, - unit_nodes=None, dim=None): + unit_nodes=None, dim=None, **kwargs): """ :arg order: the maximum total degree used for interpolation. :arg nodes: ``[ambient_dim, nelements, nunit_nodes]`` @@ -212,12 +212,13 @@ def __init__(self, order, vertex_indices, nodes, automatically assigned. """ - Record.__init__(self, + super().__init__( order=order, vertex_indices=vertex_indices, nodes=nodes, unit_nodes=unit_nodes, - element_nr_base=element_nr_base, node_nr_base=node_nr_base) + element_nr_base=element_nr_base, node_nr_base=node_nr_base, + **kwargs) def get_copy_kwargs(self, **kwargs): if "element_nr_base" not in kwargs: @@ -296,7 +297,7 @@ def __ne__(self, other): class _ModepyElementGroup(MeshElementGroup): def __init__(self, order, vertex_indices, nodes, element_nr_base=None, node_nr_base=None, - unit_nodes=None, dim=None): + unit_nodes=None, dim=None, **kwargs): """ :arg order: the maximum total degree used for interpolation. :arg nodes: ``[ambient_dim, nelements, nunit_nodes]`` @@ -325,8 +326,8 @@ def __init__(self, order, vertex_indices, nodes, raise TypeError("'dim' must be passed if 'unit_nodes' is not passed") # dim is now usable - shape = self._modepy_shape = self._modepy_shape_cls(dim) - space = self._modepy_space = mp.space_for_shape(self._modepy_shape, order) + shape = self._modepy_shape_cls(dim) + space = mp.space_for_shape(shape, order) if unit_nodes is None: unit_nodes = mp.edge_clustered_nodes_for_space(space, shape) @@ -349,14 +350,12 @@ def __init__(self, order, vertex_indices, nodes, f" got {vertex_indices.shape[-1]}") super().__init__(order, vertex_indices, nodes, - element_nr_base, node_nr_base, unit_nodes, dim) - - def __setstate__(self, valuedict): - super().__setstate__(valuedict) - - import modepy as mp - self._modepy_shape = self._modepy_shape_cls(self.dim) - self._modepy_space = mp.space_for_shape(self._modepy_shape, self.order) + element_nr_base=element_nr_base, + node_nr_base=node_nr_base, + unit_nodes=unit_nodes, + dim=dim, + _modepy_shape=shape, + _modepy_space=space) @property @memoize_method From 63f23418eab0c81a10f2031ce6cdbe6724f26ce4 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 3 Dec 2020 13:37:32 -0600 Subject: [PATCH 40/51] use modern modepy in visualization --- meshmode/discretization/visualization.py | 111 +++++++++++++++-------- 1 file changed, 71 insertions(+), 40 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 2426b25cd..e5e19b327 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -20,10 +20,14 @@ THE SOFTWARE. """ +from functools import singledispatch + import numpy as np + from pytools import memoize_method, Record from meshmode.dof_array import flatten, thaw +from modepy.shapes import Shape, Simplex, Hypercube __doc__ = """ @@ -117,6 +121,60 @@ def primitive_element_size(self): # }}} +# {{{ vtk submeshes + +@singledispatch +def vtk_submesh_for_shape(shape: Shape, node_tuples): + raise NotImplementedError(type(shape).__name__) + + +@vtk_submesh_for_shape.register(Simplex) +def _(shape: Simplex, node_tuples): + import modepy as mp + return submesh_for_shape(shape, node_tuples) + + +@vtk_submesh_for_shape.register(Hypercube) +def vtk_submesh_for_shape(shape: Hypercube, node_tuples): + node_tuple_to_index = {nt: i for i, nt in enumerate(node_tuples)} + + # NOTE: this can't use mp.submesh_for_shape because VTK vertex order is + # counterclockwise instead of z order + el_offsets = { + 1: [(0,), (1,)], + 2: [(0, 0), (1, 0), (1, 1), (0, 1)], + 3: [ + (0, 0, 0), + (1, 0, 0), + (1, 1, 0), + (0, 1, 0), + (0, 0, 1), + (1, 0, 1), + (1, 1, 1), + (0, 1, 1), + ] + }[shape.dim] + + def element_from_origin(origin): + return tuple([ + ]) + + from pytools import add_tuples + elements = [] + for origin in node_tuples: + try: + elements.append(tuple( + node_tuple_to_index[add_tuples(origin, offset)] + for offset in el_offsets + )) + except KeyError: + pass + + return elements + +# }}} + + # {{{ vtk connectivity class VTKConnectivity: @@ -155,51 +213,24 @@ def tensor_cell_types(self): } def connectivity_for_element_group(self, grp): - from pytools import ( - generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, - generate_nonnegative_integer_tuples_below as gnitb) - from meshmode.mesh import TensorProductElementGroup, SimplexElementGroup + import modepy as mp + from meshmode.mesh import _ModepyElementGroup - if isinstance(grp.mesh_el_group, SimplexElementGroup): - node_tuples = list(gnitstam(grp.order, grp.dim)) + if isinstance(grp.mesh_el_group, _ModepyElementGroup): + shape = grp.mesh_el_group._modepy_shape + space = type(grp.mesh_el_group._modepy_space)(grp.dim, grp.order) + node_tuples = mp.node_tuples_for_space(space) - from modepy.tools import simplex_submesh el_connectivity = np.array( - simplex_submesh(node_tuples), + vtk_submesh_for_shape(shape, node_tuples), dtype=np.intp) - vtk_cell_type = self.simplex_cell_types[grp.dim] - - elif isinstance(grp.mesh_el_group, TensorProductElementGroup): - node_tuples = list(gnitb(grp.order+1, grp.dim)) - node_tuple_to_index = { - nt: i for i, nt in enumerate(node_tuples)} - - def add_tuple(a, b): - return tuple(ai+bi for ai, bi in zip(a, b)) - - el_offsets = { - 1: [(0,), (1,)], - 2: [(0, 0), (1, 0), (1, 1), (0, 1)], - 3: [ - (0, 0, 0), - (1, 0, 0), - (1, 1, 0), - (0, 1, 0), - (0, 0, 1), - (1, 0, 1), - (1, 1, 1), - (0, 1, 1), - ] - }[grp.dim] - - el_connectivity = np.array([ - [ - node_tuple_to_index[add_tuple(origin, offset)] - for offset in el_offsets] - for origin in gnitb(grp.order, grp.dim)]) - - vtk_cell_type = self.tensor_cell_types[grp.dim] + if isinstance(shape, Simplex): + vtk_cell_type = self.simplex_cell_types[shape.dim] + elif isinstance(shape, Hypercube): + vtk_cell_type = self.tensor_cell_types[shape.dim] + else: + raise TypeError(f"unsupported shape: {type(shape)}") else: raise NotImplementedError("visualization for element groups " From 0f37a6ada9262cccf1b201a59f01959228c80a59 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 3 Dec 2020 13:38:08 -0600 Subject: [PATCH 41/51] use modern modepy in make_group_from_vertices --- meshmode/mesh/generation.py | 28 ++++++++++++---------------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 15bad928d..71056f209 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -333,13 +333,11 @@ def make_group_from_vertices(vertices, vertex_indices, order, # dim, nunit_nodes if unit_nodes is None: - if dim <= 3: - unit_nodes = mp.warp_and_blend_nodes(dim, order) - else: - unit_nodes = mp.equidistant_nodes(dim, order) + shape = mp.Simplex(dim) + space = mp.space_for_shape(shape, order) + unit_nodes = mp.edge_clustered_nodes_for_space(space, shape) unit_nodes_01 = 0.5 + 0.5*unit_nodes - nodes = np.einsum( "si,des->dei", unit_nodes_01, spanning_vectors) + el_origins @@ -352,24 +350,22 @@ def make_group_from_vertices(vertices, vertex_indices, order, raise ValueError("invalid number of vertices for tensor-product " "elements, must be power of two") + shape = mp.Hypercube(dim) + space = mp.space_for_shape(shape, order) + if unit_nodes is None: - from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes - unit_nodes = mp.tensor_product_nodes(dim, - legendre_gauss_lobatto_nodes(order)) + unit_nodes = mp.edge_clustered_nodes_for_space(space, shape) # shape: (dim, nnodes) unit_nodes_01 = 0.5 + 0.5*unit_nodes _, nnodes = unit_nodes.shape - from pytools import generate_nonnegative_integer_tuples_below as gnitb - # Modepy's convention is 'first coordinate varies fastest'. Adapt to it - # for the vertices. It doesn't matter for the basis. - id_tuples = [id_tuple[::-1] for id_tuple in gnitb(2, dim)] - assert len(id_tuples) == nvertices + node_tuples = mp.node_tuples_for_space(type(space)(dim, 1)) + assert len(node_tuples) == nvertices vdm = np.empty((nvertices, nvertices)) - for i, vertex_tuple in enumerate(id_tuples): - for j, func_tuple in enumerate(id_tuples): + for i, vertex_tuple in enumerate(node_tuples): + for j, func_tuple in enumerate(node_tuples): vertex_ref = np.array(vertex_tuple, dtype=np.float64) vdm[i, j] = np.prod(vertex_ref**func_tuple) @@ -379,7 +375,7 @@ def make_group_from_vertices(vertices, vertex_indices, order, coeffs[d] = la.solve(vdm, el_vertices[d].T).T vdm_nodes = np.zeros((nnodes, nvertices)) - for j, func_tuple in enumerate(id_tuples): + for j, func_tuple in enumerate(node_tuples): vdm_nodes[:, j] = np.prod( unit_nodes_01 ** np.array(func_tuple).reshape(-1, 1), axis=0) From 4465d55e9750efe6cafbd705f0aeff02313ddf60 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 3 Dec 2020 13:38:27 -0600 Subject: [PATCH 42/51] update MeshElementGroup docs --- meshmode/discretization/visualization.py | 2 +- meshmode/mesh/__init__.py | 27 +++++++++++------------- 2 files changed, 13 insertions(+), 16 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index e5e19b327..71b94f150 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -131,7 +131,7 @@ def vtk_submesh_for_shape(shape: Shape, node_tuples): @vtk_submesh_for_shape.register(Simplex) def _(shape: Simplex, node_tuples): import modepy as mp - return submesh_for_shape(shape, node_tuples) + return mp.submesh_for_shape(shape, node_tuples) @vtk_submesh_for_shape.register(Hypercube) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 1f64dc238..328eefdc1 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -202,11 +202,10 @@ def __init__(self, order, vertex_indices, nodes, unit_nodes=None, dim=None, **kwargs): """ :arg order: the maximum total degree used for interpolation. - :arg nodes: ``[ambient_dim, nelements, nunit_nodes]`` + :arg nodes: ``(ambient_dim, nelements, nunit_nodes)`` The nodes are assumed to be mapped versions of *unit_nodes*. - :arg unit_nodes: ``[dim, nunit_nodes]`` - The unit nodes of which *nodes* is a mapped - version. + :arg unit_nodes: ``(dim, nunit_nodes)`` The unit nodes of which *nodes* + is a mapped version. Do not supply *element_nr_base* and *node_nr_base*, they will be automatically assigned. @@ -257,20 +256,20 @@ def nunit_nodes(self): @property def is_affine(self): - raise NotImplementedError() + raise NotImplementedError def face_vertex_indices(self): """Return a tuple of tuples indicating which vertices (in mathematically positive ordering) make up each face of an element in this group. """ - raise NotImplementedError() + raise NotImplementedError def vertex_unit_coordinates(self): """Return an array of shape ``(nfaces, dim)`` with the unit coordinates of each vertex. """ - raise NotImplementedError() + raise NotImplementedError @property def nfaces(self): @@ -300,15 +299,13 @@ def __init__(self, order, vertex_indices, nodes, unit_nodes=None, dim=None, **kwargs): """ :arg order: the maximum total degree used for interpolation. - :arg nodes: ``[ambient_dim, nelements, nunit_nodes]`` + :arg nodes: ``(ambient_dim, nelements, nunit_nodes)`` The nodes are assumed to be mapped versions of *unit_nodes*. - :arg unit_nodes: ``[dim, nunit_nodes]`` - The unit nodes of which *nodes* is a mapped - version. If unspecified, the nodes from - :func:`modepy.warp_and_blend_nodes` for *dim* - are assumed. These must be in unit coordinates - as defined in :mod:`modepy`. - :arg dim: only used if *unit_nodes* is None, to get + :arg unit_nodes: ``(dim, nunit_nodes)`` The unit nodes of which + *nodes* is a mapped version. If unspecified, the nodes from + :func:`modepy.edge_clustered_nodes_for_space` are assumed. + These must be in unit coordinates as defined in :mod:`modepy`. + :arg dim: only used if *unit_nodes* is *None*, to get the default unit nodes. Do not supply *element_nr_base* and *node_nr_base*, they will be From f99d6ebb16f37f6b352a3405bd2858187c7c6eb2 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 3 Dec 2020 15:03:47 -0600 Subject: [PATCH 43/51] use new tensor_product_nodes in poly_element --- meshmode/discretization/poly_element.py | 60 +++++++++++++++---------- 1 file changed, 36 insertions(+), 24 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 18a1f8804..470c5db5e 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -376,13 +376,20 @@ def from_mesh_interp_matrix(self): class TensorProductElementGroupBase(PolynomialElementGroupBase, HypercubeElementGroupBase): - def __init__(self, mesh_el_group, order, index, *, - basis, unit_nodes_1d, quad_weights_1d=None): + def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes_1d): + """ + :arg basis: a :class:`modepy.Basis`. + :arg unit_nodes_1d: a :class:`tuple` of one-dimensional unit nodes, + one of each dimension of *mesh_el_group* + """ super().__init__(mesh_el_group, order, index) + if len(unit_nodes_1d) != mesh_el_group.dim: + raise ValueError("insufficient nodes for element group dimension. " + f"expected {mesh_le_group.dim}, got {len(unit_nodes_1d)}.") + self._basis = basis self._unit_nodes_1d = unit_nodes_1d - self._quad_weights_1d = quad_weights_1d def is_orthogonal_basis(self): try: @@ -403,26 +410,18 @@ def grad_basis(self): @property @memoize_method def unit_nodes(self): - # FIXME: allow non-homogeneous tensor products - return mp.tensor_product_nodes(self.dim, self._unit_nodes_1d) + return mp.tensor_product_nodes(self._unit_nodes_1d) @property @memoize_method def weights(self): - if self._quad_weights_1d is None: - return np.dot( - self.mass_matrix(), - np.ones(len(self.basis()))) - else: - # FIXME: allow non-homogeneous tensor products - return np.prod( - mp.tensor_product_nodes(self.dim, self._quad_weights_1d), - axis=0) + return np.dot( + self.mass_matrix(), + np.ones(len(self.basis()))) class LegendreTensorProductElementGroup(TensorProductElementGroupBase): - def __init__(self, mesh_el_group, order, index, *, - unit_nodes_1d=None, quad_weights_1d=None): + def __init__(self, mesh_el_group, order, index, *, unit_nodes_1d): basis = mp.orthonormal_basis_for_space( mp.QN(mesh_el_group.dim, order), @@ -430,8 +429,7 @@ def __init__(self, mesh_el_group, order, index, *, super().__init__(mesh_el_group, order, index, basis=basis, - unit_nodes_1d=unit_nodes_1d, - quad_weights_1d=quad_weights_1d) + unit_nodes_1d=unit_nodes_1d) class GaussLegendreTensorProductElementGroup(LegendreTensorProductElementGroup): @@ -442,11 +440,21 @@ class GaussLegendreTensorProductElementGroup(LegendreTensorProductElementGroup): No interpolation nodes are present on the boundary of the hypercube. """ def __init__(self, mesh_el_group, order, index): - quad = mp.LegendreGaussQuadrature(order) - super().__init__(mesh_el_group, order, index, - unit_nodes_1d=quad.nodes, - quad_weights_1d=quad.weights) + unit_nodes_1d=[None] * mesh_el_group.dim) + + @property + @memoize_method + def _quadrature_rule(self): + return mp.LegendreGaussTensorProductQuadrature(self.order, self.dim) + + @property + def unit_nodes(self): + return self._quadrature_rule.nodes + + @property + def weights(self): + return self._quadrature_rule.weights class LegendreGaussLobattoTensorProductElementGroup( @@ -461,8 +469,10 @@ class LegendreGaussLobattoTensorProductElementGroup( def __init__(self, mesh_el_group, order, index): from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes + nodes_1d = legendre_gauss_lobatto_nodes(order) + super().__init__(mesh_el_group, order, index, - unit_nodes_1d=legendre_gauss_lobatto_nodes(order)) + unit_nodes_1d=[nodes_1d] * mesh_el_group.dim) class EquidistantTensorProductElementGroup(LegendreTensorProductElementGroup): @@ -476,8 +486,10 @@ class EquidistantTensorProductElementGroup(LegendreTensorProductElementGroup): def __init__(self, mesh_el_group, order, index): from modepy.nodes import equidistant_nodes + nodes_1d = equidistant_nodes(1, order)[0] + super().__init__(mesh_el_group, order, index, - unit_nodes_1d=equidistant_nodes(1, order)[0]) + unit_nodes_1d=[nodes_1d] * mesh_el_group.dim) # }}} From bef4afb3c823085153aa2eb40620ee431dc682a3 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 3 Dec 2020 15:19:18 -0600 Subject: [PATCH 44/51] fix shadowing --- meshmode/discretization/poly_element.py | 2 +- meshmode/discretization/visualization.py | 6 +----- test/test_meshmode.py | 16 ++++++++-------- 3 files changed, 10 insertions(+), 14 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 470c5db5e..a7e128060 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -386,7 +386,7 @@ def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes_1d): if len(unit_nodes_1d) != mesh_el_group.dim: raise ValueError("insufficient nodes for element group dimension. " - f"expected {mesh_le_group.dim}, got {len(unit_nodes_1d)}.") + f"expected {mesh_el_group.dim}, got {len(unit_nodes_1d)}.") self._basis = basis self._unit_nodes_1d = unit_nodes_1d diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 71b94f150..1b2643ace 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -135,7 +135,7 @@ def _(shape: Simplex, node_tuples): @vtk_submesh_for_shape.register(Hypercube) -def vtk_submesh_for_shape(shape: Hypercube, node_tuples): +def _(shape: Hypercube, node_tuples): node_tuple_to_index = {nt: i for i, nt in enumerate(node_tuples)} # NOTE: this can't use mp.submesh_for_shape because VTK vertex order is @@ -155,10 +155,6 @@ def vtk_submesh_for_shape(shape: Hypercube, node_tuples): ] }[shape.dim] - def element_from_origin(origin): - return tuple([ - ]) - from pytools import add_tuples elements = [] for origin in node_tuples: diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 888070337..3e7328512 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -144,20 +144,20 @@ def Get_size(self): # noqa: N802 assert(filecmp.cmp("ref-"+pvtu_filename, pvtu_filename)) -@pytest.mark.parametrize(("dim", "group_factory"), [ +@pytest.mark.parametrize(("dim", "group_cls"), [ (1, SimplexElementGroup), (2, SimplexElementGroup), (3, SimplexElementGroup), (2, TensorProductElementGroup), (3, TensorProductElementGroup), ]) -def test_visualizers(actx_factory, dim, group_factory): +def test_visualizers(actx_factory, dim, group_cls): actx = actx_factory() nelements = 64 target_order = 4 - is_simplex = issubclass(group_factory, SimplexElementGroup) + is_simplex = issubclass(group_cls, SimplexElementGroup) if dim == 1: mesh = mgen.make_curve_mesh( mgen.NArmedStarfish(5, 0.25), @@ -169,7 +169,7 @@ def test_visualizers(actx_factory, dim, group_factory): else: mesh = mgen.generate_regular_rect_mesh( a=(0,)*dim, b=(1,)*dim, n=(5,)*dim, - group_factory=group_factory, + group_cls=group_cls, order=target_order) elif dim == 3: if is_simplex: @@ -177,18 +177,18 @@ def test_visualizers(actx_factory, dim, group_factory): else: mesh = mgen.generate_regular_rect_mesh( a=(0,)*dim, b=(1,)*dim, n=(5,)*dim, - group_factory=group_factory, + group_cls=group_cls, order=target_order) else: raise ValueError("unknown dimensionality") if is_simplex: - discr_group_factory = InterpolatoryQuadratureSimplexGroupFactory + group_factory = InterpolatoryQuadratureSimplexGroupFactory else: - discr_group_factory = LegendreGaussLobattoTensorProductGroupFactory + group_factory = LegendreGaussLobattoTensorProductGroupFactory from meshmode.discretization import Discretization - discr = Discretization(actx, mesh, discr_group_factory(target_order)) + discr = Discretization(actx, mesh, group_factory(target_order)) nodes = thaw(actx, discr.nodes()) f = actx.np.sqrt(sum(nodes**2)) + 1j*nodes[0] From ec213e6ee23ebd8d1d5f8544d399c1d844b69227 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Thu, 10 Dec 2020 18:33:41 -0600 Subject: [PATCH 45/51] Point modepy back to master in req.txt --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index c65ea9f71..7027fed52 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,7 +3,7 @@ recursivenodes git+https://github.com/inducer/pytools.git#egg=pytools git+https://github.com/inducer/gmsh_interop.git#egg=gmsh_interop git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile -git+https://github.com/alexfikl/modepy.git@tensor-face-mass-matrix#egg=modepy +git+https://github.com/inducer/modepy.git#egg=modepy git+https://github.com/inducer/pyopencl.git#egg=pyopencl git+https://github.com/inducer/islpy.git#egg=islpy From 79fcade771777a244b4b1deb425214f3b358afbd Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 21 Dec 2020 13:57:31 -0600 Subject: [PATCH 46/51] handle all cases in if condition MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- test/test_meshmode.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 3e7328512..2054ddd07 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -332,6 +332,8 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False): num_on_bdy = dim * (nelem-1)**(dim-1) elif group_cls is SimplexElementGroup: num_on_bdy = dim * (dim-1) * (nelem-1)**(dim-1) + else: + assert False assert not is_boundary_tag_empty(mesh, "btag_test_1") assert not is_boundary_tag_empty(mesh, "btag_test_2") From adbc99025732c663cd474017a1843b68fd150101 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 21 Dec 2020 14:13:37 -0600 Subject: [PATCH 47/51] use stricter orthogonality check in tensor products MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/discretization/poly_element.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index a7e128060..ef5e93c35 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -393,8 +393,7 @@ def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes_1d): def is_orthogonal_basis(self): try: - self._basis.orthonormality_weight() - return True + return self._basis.orthonormality_weight() == 1 except mp.BasisNotOrthonormal: return False From e401f6132360e3a7a0a80e4986d39bc81546480e Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 21 Dec 2020 14:16:11 -0600 Subject: [PATCH 48/51] add comment for tensor product is_orthogonal check --- meshmode/discretization/poly_element.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index ef5e93c35..d4c1f63c0 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -393,6 +393,8 @@ def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes_1d): def is_orthogonal_basis(self): try: + # NOTE: meshmode kind of assumes that the basis is also orthonormal, + # which is why this check is stricter than expected return self._basis.orthonormality_weight() == 1 except mp.BasisNotOrthonormal: return False From d9fbabb631d02743d9ec68fb7436de20a5471da4 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 21 Dec 2020 14:33:27 -0600 Subject: [PATCH 49/51] take full unit_nodes in tensor product element groups --- meshmode/discretization/poly_element.py | 58 +++++++++++++------------ 1 file changed, 31 insertions(+), 27 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index d4c1f63c0..17ed7a93b 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -376,20 +376,24 @@ def from_mesh_interp_matrix(self): class TensorProductElementGroupBase(PolynomialElementGroupBase, HypercubeElementGroupBase): - def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes_1d): + def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes): """ - :arg basis: a :class:`modepy.Basis`. - :arg unit_nodes_1d: a :class:`tuple` of one-dimensional unit nodes, - one of each dimension of *mesh_el_group* + :arg basis: a :class:`modepy.TensorProductBasis`. + :arg unit_nodes: unit nodes for the tensor product, obtained by + using :func:`modepy.tensor_product_nodes`, for example. """ super().__init__(mesh_el_group, order, index) - if len(unit_nodes_1d) != mesh_el_group.dim: - raise ValueError("insufficient nodes for element group dimension. " - f"expected {mesh_el_group.dim}, got {len(unit_nodes_1d)}.") + if basis._dim != mesh_el_group.dim: + raise ValueError("basis dimension does not match element group: " + f"expected {mesh_el_group.dim}, got {basis._dim}.") + + if unit_nodes.shape[0] != mesh_el_group.dim: + raise ValueError("unit node dimension does not match element group: " + f"expected {mesh_el_group.dim}, got {unit_nodes.shape[0]}.") self._basis = basis - self._unit_nodes_1d = unit_nodes_1d + self._unit_nodes = unit_nodes def is_orthogonal_basis(self): try: @@ -399,19 +403,21 @@ def is_orthogonal_basis(self): except mp.BasisNotOrthonormal: return False + @memoize_method def mode_ids(self): return self._basis.mode_ids + @memoize_method def basis(self): return self._basis.functions + @memoize_method def grad_basis(self): return self._basis.gradients @property - @memoize_method def unit_nodes(self): - return mp.tensor_product_nodes(self._unit_nodes_1d) + return self._unit_nodes @property @memoize_method @@ -422,15 +428,14 @@ def weights(self): class LegendreTensorProductElementGroup(TensorProductElementGroupBase): - def __init__(self, mesh_el_group, order, index, *, unit_nodes_1d): - + def __init__(self, mesh_el_group, order, index, *, unit_nodes): basis = mp.orthonormal_basis_for_space( mp.QN(mesh_el_group.dim, order), mp.Hypercube(mesh_el_group.dim)) super().__init__(mesh_el_group, order, index, basis=basis, - unit_nodes_1d=unit_nodes_1d) + unit_nodes=unit_nodes) class GaussLegendreTensorProductElementGroup(LegendreTensorProductElementGroup): @@ -440,18 +445,13 @@ class GaussLegendreTensorProductElementGroup(LegendreTensorProductElementGroup): No interpolation nodes are present on the boundary of the hypercube. """ - def __init__(self, mesh_el_group, order, index): - super().__init__(mesh_el_group, order, index, - unit_nodes_1d=[None] * mesh_el_group.dim) - @property - @memoize_method - def _quadrature_rule(self): - return mp.LegendreGaussTensorProductQuadrature(self.order, self.dim) + def __init__(self, mesh_el_group, order, index): + self._quadrature_rule = mp.LegendreGaussTensorProductQuadrature( + order, mesh_el_group.dim) - @property - def unit_nodes(self): - return self._quadrature_rule.nodes + super().__init__(mesh_el_group, order, index, + unit_nodes=self._quadrature_rule.nodes) @property def weights(self): @@ -470,10 +470,12 @@ class LegendreGaussLobattoTensorProductElementGroup( def __init__(self, mesh_el_group, order, index): from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes - nodes_1d = legendre_gauss_lobatto_nodes(order) + unit_nodes_1d = legendre_gauss_lobatto_nodes(order) super().__init__(mesh_el_group, order, index, - unit_nodes_1d=[nodes_1d] * mesh_el_group.dim) + unit_nodes=mp.tensor_product_nodes( + [unit_nodes_1d] * mesh_el_group.dim) + ) class EquidistantTensorProductElementGroup(LegendreTensorProductElementGroup): @@ -487,10 +489,12 @@ class EquidistantTensorProductElementGroup(LegendreTensorProductElementGroup): def __init__(self, mesh_el_group, order, index): from modepy.nodes import equidistant_nodes - nodes_1d = equidistant_nodes(1, order)[0] + unit_nodes_1d = equidistant_nodes(1, order)[0] super().__init__(mesh_el_group, order, index, - unit_nodes_1d=[nodes_1d] * mesh_el_group.dim) + unit_nodes=mp.tensor_product_nodes( + [unit_nodes_1d] * mesh_el_group.dim) + ) # }}} From f2b33e389dc105dbe325247dd7f1ec3302386178 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 21 Dec 2020 15:11:17 -0600 Subject: [PATCH 50/51] rename some variables --- meshmode/discretization/connection/face.py | 2 +- meshmode/mesh/__init__.py | 4 ++-- meshmode/mesh/generation.py | 10 +++++----- meshmode/mesh/io.py | 3 ++- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index b0c855531..e3494c019 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -261,7 +261,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag, # group by face_index - for face in mgrp.faces: + for face in mgrp._modepy_faces: batch_boundary_el_numbers_in_grp = np.array([ ibface_el for ibface_el, ibface_face in group_boundary_faces diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 328eefdc1..ea298b185 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -356,11 +356,11 @@ def __init__(self, order, vertex_indices, nodes, @property @memoize_method - def faces(self): + def _modepy_faces(self): return mp.faces_for_shape(self._modepy_shape) def face_vertex_indices(self): - return tuple(face.volume_vertex_indices for face in self.faces) + return tuple(face.volume_vertex_indices for face in self._modepy_faces) def vertex_unit_coordinates(self): return mp.unit_vertices_for_shape(self._modepy_shape).T diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 71056f209..27ea045ac 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -360,12 +360,12 @@ def make_group_from_vertices(vertices, vertex_indices, order, unit_nodes_01 = 0.5 + 0.5*unit_nodes _, nnodes = unit_nodes.shape - node_tuples = mp.node_tuples_for_space(type(space)(dim, 1)) - assert len(node_tuples) == nvertices + vertex_tuples = mp.node_tuples_for_space(type(space)(dim, 1)) + assert len(vertex_tuples) == nvertices vdm = np.empty((nvertices, nvertices)) - for i, vertex_tuple in enumerate(node_tuples): - for j, func_tuple in enumerate(node_tuples): + for i, vertex_tuple in enumerate(vertex_tuples): + for j, func_tuple in enumerate(vertex_tuples): vertex_ref = np.array(vertex_tuple, dtype=np.float64) vdm[i, j] = np.prod(vertex_ref**func_tuple) @@ -375,7 +375,7 @@ def make_group_from_vertices(vertices, vertex_indices, order, coeffs[d] = la.solve(vdm, el_vertices[d].T).T vdm_nodes = np.zeros((nnodes, nvertices)) - for j, func_tuple in enumerate(node_tuples): + for j, func_tuple in enumerate(vertex_tuples): vdm_nodes[:, j] = np.prod( unit_nodes_01 ** np.array(func_tuple).reshape(-1, 1), axis=0) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index e8a24d941..d0ac94f7d 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -233,7 +233,8 @@ def get_mesh(self): unit_nodes=unit_nodes ) else: - pass + # NOTE: already checked above + assert False groups.append(group) From 0c10f6be58b0df1bec1f409e5e58e10e911fa4da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Wed, 23 Dec 2020 16:46:29 -0600 Subject: [PATCH 51/51] Tweak is_orthogonal_basis comment --- meshmode/discretization/poly_element.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 17ed7a93b..5714ca8df 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -397,8 +397,8 @@ def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes): def is_orthogonal_basis(self): try: - # NOTE: meshmode kind of assumes that the basis is also orthonormal, - # which is why this check is stricter than expected + # NOTE: meshmode kind of assumes that the basis is orthonormal + # with weight 1, which is why this check is stricter than expected. return self._basis.orthonormality_weight() == 1 except mp.BasisNotOrthonormal: return False