diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index cef3d66785d24df7f8b764139baa810a39fbe473..5f636e78165f46b32ecb9c41e2883927ece81f82 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -26,7 +26,9 @@ THE SOFTWARE. import numpy as np #import numpy.linalg as la from pytools import memoize_method -from meshmode.mesh import SimplexElementGroup as _MeshSimplexElementGroup +from meshmode.mesh import ( + SimplexElementGroup as _MeshSimplexElementGroup, + TensorProductElementGroup as _MeshTensorProductElementGroup) import modepy as mp @@ -38,7 +40,8 @@ Group types .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup .. autoclass:: PolynomialWarpAndBlendElementGroup -.. autoclass:: PolynomialEquidistantElementGroup +.. autoclass:: PolynomialEquidistantSimplexElementGroup +.. autoclass:: LegendreGaussLobattoTensorProductElementGroup Group factories ^^^^^^^^^^^^^^^ @@ -46,7 +49,8 @@ Group factories .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory .. autoclass:: PolynomialWarpAndBlendGroupFactory -.. autoclass:: PolynomialEquidistantGroupFactory +.. autoclass:: PolynomialEquidistantSimplexGroupFactory +.. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ # FIXME Most of the loopy kernels will break as soon as we start using multiple @@ -61,24 +65,9 @@ Group factories from meshmode.discretization import ElementGroupBase -# {{{ concrete element groups - -class PolynomialSimplexElementGroupBase(ElementGroupBase): - def is_orthogonal_basis(self): - return self.dim <= 3 - - def basis(self): - if self.dim <= 3: - return mp.simplex_onb(self.dim, self.order) - else: - return mp.simplex_monomial_basis(self.dim, self.order) - - 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) +# {{{ base class for poynomial elements +class PolynomialElementGroupBase(ElementGroupBase): @memoize_method def mass_matrix(self): assert self.is_orthogonal_basis() @@ -88,14 +77,6 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase): self.basis(), self.unit_nodes) - @memoize_method - def from_mesh_interp_matrix(self): - meg = self.mesh_el_group - return mp.resampling_matrix( - mp.simplex_best_available_basis(meg.dim, meg.order), - self.unit_nodes, - meg.unit_nodes) - @memoize_method def diff_matrices(self): if len(self.basis()) != self.unit_nodes.shape[1]: @@ -116,6 +97,35 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase): else: return result +# }}} + + +# {{{ concrete element groups for simplices + +class PolynomialSimplexElementGroupBase(PolynomialElementGroupBase): + def is_orthogonal_basis(self): + return self.dim <= 3 + + def basis(self): + if self.dim <= 3: + return mp.simplex_onb(self.dim, self.order) + else: + return mp.simplex_monomial_basis(self.dim, self.order) + + 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) + + @memoize_method + def from_mesh_interp_matrix(self): + meg = self.mesh_el_group + return mp.resampling_matrix( + mp.simplex_best_available_basis(meg.dim, meg.order), + self.unit_nodes, + meg.unit_nodes) + class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): """Elemental discretization supplying a high-order quadrature rule @@ -189,7 +199,7 @@ class _MassMatrixQuadratureElementGroup(PolynomialSimplexElementGroupBase): @memoize_method def weights(self): return np.dot( - mp.mass_matrix(self.basis(), self.unit_nodes), + self.mass_matrix(), np.ones(len(self.basis()))) @@ -210,7 +220,7 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): return result -class PolynomialEquidistantElementGroup(_MassMatrixQuadratureElementGroup): +class PolynomialEquidistantSimplexElementGroup(_MassMatrixQuadratureElementGroup): """Elemental discretization with a number of nodes matching the number of polynomials in :math:`P^k`, hence usable for differentiation and interpolation. Interpolation nodes are present on the boundary of the @@ -231,13 +241,54 @@ class PolynomialEquidistantElementGroup(_MassMatrixQuadratureElementGroup): # }}} +# {{{ concrete element groups for tensor product elements + +class LegendreGaussLobattoTensorProductElementGroup(PolynomialElementGroupBase): + def is_orthogonal_basis(self): + return True + + def basis(self): + from modepy.modes import tensor_product_basis, jacobi + from functools import partial + return tensor_product_basis( + self.dim, tuple( + partial(jacobi, 0, 0, i) + for i in range(self.order+1))) + + def grad_basis(self): + raise NotImplementedError() + + @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)) + + @memoize_method + def from_mesh_interp_matrix(self): + from modepy.modes import tensor_product_basis, jacobi + from functools import partial + meg = self.mesh_el_group + + basis = tensor_product_basis( + self.dim, tuple( + partial(jacobi, 0, 0, i) + for i in range(meg.order+1))) + + return mp.resampling_matrix(basis, self.unit_nodes, meg.unit_nodes) + +# }}} + + # {{{ group factories class ElementGroupFactory(object): pass -class OrderBasedGroupFactory(ElementGroupFactory): +class HomogeneousOrderBasedGroupFactory(ElementGroupFactory): def __init__(self, order): self.order = order @@ -249,28 +300,58 @@ class OrderBasedGroupFactory(ElementGroupFactory): return self.group_class(mesh_el_group, self.order, node_nr_base) -class InterpolatoryQuadratureSimplexGroupFactory(OrderBasedGroupFactory): +class OrderAndTypeBasedGroupFactory(ElementGroupFactory): + def __init__(self, order, simplex_group_class, tensor_product_group_class): + self.order = order + self.simplex_group_class = simplex_group_class + self.tensor_product_group_class = tensor_product_group_class + + def __call__(self, mesh_el_group, node_nr_base): + if isinstance(mesh_el_group, _MeshSimplexElementGroup): + group_class = self.simplex_group_class + elif isinstance(mesh_el_group, _MeshTensorProductElementGroup): + group_class = self.tensor_product_group_class + else: + raise TypeError("only mesh element groups of type '%s' and '%s' " + "are supported" % ( + _MeshSimplexElementGroup.__name__, + _MeshTensorProductElementGroup.__name__)) + + return group_class(mesh_el_group, self.order, node_nr_base) + + +# {{{ group factories for simplices + +class InterpolatoryQuadratureSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = InterpolatoryQuadratureSimplexElementGroup -class QuadratureSimplexGroupFactory(OrderBasedGroupFactory): +class QuadratureSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = QuadratureSimplexElementGroup -class PolynomialWarpAndBlendGroupFactory(OrderBasedGroupFactory): +class PolynomialWarpAndBlendGroupFactory(HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = PolynomialWarpAndBlendElementGroup -class PolynomialEquidistantGroupFactory(OrderBasedGroupFactory): +class PolynomialEquidistantSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): """ .. versionadded:: 2016.1 """ mesh_group_class = _MeshSimplexElementGroup - group_class = PolynomialEquidistantElementGroup + group_class = PolynomialEquidistantSimplexElementGroup + +# }}} + + +class LegendreGaussLobattoTensorProductGroupFactory( + HomogeneousOrderBasedGroupFactory): + mesh_group_class = _MeshTensorProductElementGroup + group_class = LegendreGaussLobattoTensorProductElementGroup # }}} diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index d68ec0c58078e68417721c8205fc3f24021d0bad..1021d8f4f7556e5306fdd74703dcb94d531b4c5c 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -1,6 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" @@ -24,8 +22,9 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range import numpy as np -from pytools import memoize_method +from pytools import memoize_method, Record import pyopencl as cl __doc__ = """ @@ -71,17 +70,50 @@ def separate_by_real_and_imag(data, real_only): yield (name, field) +class _VisConnectivityGroup(Record): + """ + .. attribute:: vis_connectivity + + an array of shape ``(group.nelements,nsubelements,primitive_element_size)`` + + .. attribute:: vtk_cell_type + + .. attribute:: subelement_nr_base + """ + + @property + def nsubelements(self): + return self.nelements * self.nsubelements_per_element + + @property + def nelements(self): + return self.vis_connectivity.shape[0] + + @property + def nsubelements_per_element(self): + return self.vis_connectivity.shape[1] + + @property + def primitive_element_size(self): + return self.vis_connectivity.shape[2] + + class Visualizer(object): """ .. automethod:: show_scalar_in_mayavi .. automethod:: write_vtk_file """ - def __init__(self, connection): + def __init__(self, connection, element_shrink_factor=None): self.connection = connection self.discr = connection.from_discr self.vis_discr = connection.to_discr + if element_shrink_factor is None: + element_shrink_factor = 1 + + self.element_shrink_factor = element_shrink_factor + def _resample_and_get(self, queue, vec): from pytools.obj_array import with_object_array_or_scalar @@ -90,42 +122,101 @@ class Visualizer(object): return with_object_array_or_scalar(resample_and_get_one, vec) + # {{{ vis sub-element connectivity + @memoize_method def _vis_connectivity(self): """ - :return: an array of shape - ``(vis_discr.nelements,nsubelements,primitive_element_size)`` + :return: a list of :class:`_VisConnectivityGroup` instances. """ # Assume that we're using modepy's default node ordering. - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam, single_valued - vis_order = single_valued( - group.order for group in self.vis_discr.groups) - node_tuples = list(gnitstam(vis_order, self.vis_discr.dim)) + 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 + + result = [] - from modepy.tools import submesh - el_connectivity = np.array( - submesh(node_tuples), - dtype=np.intp) + from pyvisfile.vtk import ( + VTK_LINE, VTK_TRIANGLE, VTK_TETRA, + VTK_QUAD, VTK_HEXAHEDRON) - nelements = sum(group.nelements for group in self.vis_discr.groups) - vis_connectivity = np.empty( - (nelements,) + el_connectivity.shape, dtype=np.intp) + subel_nr_base = 0 - el_nr_base = 0 for group in self.vis_discr.groups: + if isinstance(group.mesh_el_group, SimplexElementGroup): + node_tuples = list(gnitstam(group.order, group.dim)) + + from modepy.tools import submesh + el_connectivity = np.array( + submesh(node_tuples), + dtype=np.intp) + vtk_cell_type = { + 1: VTK_LINE, + 2: VTK_TRIANGLE, + 3: VTK_TETRA, + }[group.dim] + + elif isinstance(group.mesh_el_group, TensorProductElementGroup): + node_tuples = list(gnitb(group.order+1, group.dim)) + node_tuple_to_index = dict( + (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), + ] + }[group.dim] + + el_connectivity = np.array([ + [ + node_tuple_to_index[add_tuple(origin, offset)] + for offset in el_offsets] + for origin in gnitb(group.order, group.dim)]) + + vtk_cell_type = { + 1: VTK_LINE, + 2: VTK_QUAD, + 3: VTK_HEXAHEDRON, + }[group.dim] + + else: + raise NotImplementedError("visualization for element groups " + "of type '%s'" % type(group.mesh_el_group).__name__) + assert len(node_tuples) == group.nunit_nodes - vis_connectivity[el_nr_base:el_nr_base+group.nelements] = ( - np.arange( - el_nr_base*group.nunit_nodes, - (el_nr_base+group.nelements)*group.nunit_nodes, - group.nunit_nodes + vis_connectivity = ( + group.node_nr_base + np.arange( + 0, group.nelements*group.nunit_nodes, group.nunit_nodes )[:, np.newaxis, np.newaxis] - + el_connectivity) - el_nr_base += group.nelements + + el_connectivity).astype(np.intp) + + vgrp = _VisConnectivityGroup( + vis_connectivity=vis_connectivity, + vtk_cell_type=vtk_cell_type, + subelement_nr_base=subel_nr_base) + result.append(vgrp) + + subel_nr_base += vgrp.nsubelements + + return result + + # }}} - return vis_connectivity + # {{{ mayavi def show_scalar_in_mayavi(self, field, **kwargs): import mayavi.mlab as mlab @@ -140,7 +231,7 @@ class Visualizer(object): assert nodes.shape[0] == self.vis_discr.ambient_dim #mlab.points3d(nodes[0], nodes[1], 0*nodes[0]) - vis_connectivity = self._vis_connectivity() + vis_connectivity, = self._vis_connectivity() if self.vis_discr.dim == 1: nodes = list(nodes) @@ -153,6 +244,7 @@ class Visualizer(object): # http://docs.enthought.com/mayavi/mayavi/auto/example_plotting_many_lines.html # noqa src = mlab.pipeline.scalar_scatter(*args) + src.mlab_source.dataset.lines = vis_connectivity.reshape(-1, 2) lines = mlab.pipeline.stripper(src) mlab.pipeline.surface(lines, **kwargs) @@ -175,21 +267,17 @@ class Visualizer(object): if do_show: mlab.show() + # }}} + + # {{{ vtk + def write_vtk_file(self, file_name, names_and_fields, compressor=None, real_only=False): from pyvisfile.vtk import ( UnstructuredGrid, DataArray, AppendedDataXMLGenerator, - VTK_LINE, VTK_TRIANGLE, VTK_TETRA, VF_LIST_OF_COMPONENTS) - el_types = { - 1: VTK_LINE, - 2: VTK_TRIANGLE, - 3: VTK_TETRA, - } - - el_type = el_types[self.vis_discr.dim] with cl.CommandQueue(self.vis_discr.cl_context) as queue: nodes = self.vis_discr.nodes().with_queue(queue).get() @@ -198,20 +286,41 @@ class Visualizer(object): (name, self._resample_and_get(queue, fld)) for name, fld in names_and_fields] - connectivity = self._vis_connectivity() + vc_groups = self._vis_connectivity() + + # {{{ create cell_types + + nsubelements = sum(vgrp.nsubelements for vgrp in vc_groups) + cell_types = np.empty(nsubelements, dtype=np.uint8) + cell_types.fill(255) + for vgrp in vc_groups: + cell_types[ + vgrp.subelement_nr_base: + vgrp.subelement_nr_base + vgrp.nsubelements] = \ + vgrp.vtk_cell_type + assert (cell_types < 255).all() + + # }}} - nprimitive_elements = ( - connectivity.shape[0] - * connectivity.shape[1]) + nodes = nodes + if self.element_shrink_factor != 1: + for vgrp in self.vis_discr.groups: + nodes_view = vgrp.view(nodes) + el_centers = np.mean(nodes_view, axis=-1) + nodes_view[:] = ( + (self.element_shrink_factor * nodes_view) + + (1-self.element_shrink_factor) + * el_centers[:, :, np.newaxis]) grid = UnstructuredGrid( (self.vis_discr.nnodes, DataArray("points", nodes.reshape(self.vis_discr.ambient_dim, -1), vector_format=VF_LIST_OF_COMPONENTS)), - cells=connectivity.reshape(-1), - cell_types=np.asarray([el_type] * nprimitive_elements, - dtype=np.uint8)) + cells=np.hstack([ + vgrp.vis_connectivity.reshape(-1) + for vgrp in vc_groups]), + cell_types=cell_types) # for name, field in separate_by_real_and_imag(cell_data, real_only): # grid.add_celldata(DataArray(name, field, @@ -230,18 +339,26 @@ class Visualizer(object): AppendedDataXMLGenerator(compressor)(grid).write(outf) -def make_visualizer(queue, discr, vis_order): +def make_visualizer(queue, discr, vis_order, element_shrink_factor=None): from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory + from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendElementGroup, + LegendreGaussLobattoTensorProductElementGroup, + OrderAndTypeBasedGroupFactory) vis_discr = Discretization( discr.cl_context, discr.mesh, - PolynomialWarpAndBlendGroupFactory(vis_order), + OrderAndTypeBasedGroupFactory( + vis_order, + simplex_group_class=PolynomialWarpAndBlendElementGroup, + tensor_product_group_class=( + LegendreGaussLobattoTensorProductElementGroup)), real_dtype=discr.real_dtype) from meshmode.discretization.connection import \ make_same_mesh_connection - return Visualizer(make_same_mesh_connection(vis_discr, discr)) + return Visualizer( + make_same_mesh_connection(vis_discr, discr), + element_shrink_factor=element_shrink_factor) # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index a793743ce6ff50cc67aea4b31225a1b9e4df4e3d..5178cf406512547e0bd1424ba983e0759a9eb617 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -36,7 +36,7 @@ from pyopencl.tools import ( # noqa from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, - PolynomialEquidistantGroupFactory, + PolynomialEquidistantSimplexGroupFactory, ) from meshmode.mesh import BTAG_ALL from meshmode.discretization.connection import \ @@ -454,32 +454,49 @@ def test_element_orientation(): def test_merge_and_map(ctx_getter, visualize=False): from meshmode.mesh.io import generate_gmsh, FileSource + from meshmode.mesh.generation import generate_box_mesh + from meshmode.mesh import TensorProductElementGroup + from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory) mesh_order = 3 - mesh = generate_gmsh( - FileSource("blob-2d.step"), 2, order=mesh_order, - force_ambient_dim=2, - other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"] - ) + if 1: + mesh = generate_gmsh( + FileSource("blob-2d.step"), 2, order=mesh_order, + force_ambient_dim=2, + other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"] + ) + + discr_grp_factory = PolynomialWarpAndBlendGroupFactory(3) + else: + mesh = generate_box_mesh( + ( + np.linspace(0, 1, 4), + np.linspace(0, 1, 4), + np.linspace(0, 1, 4), + ), + 10, group_factory=TensorProductElementGroup) + + discr_grp_factory = LegendreGaussLobattoTensorProductGroupFactory(3) from meshmode.mesh.processing import merge_disjoint_meshes, affine_map - mesh2 = affine_map(mesh, A=np.eye(2), b=np.array([5, 0])) + mesh2 = affine_map(mesh, + A=np.eye(mesh.ambient_dim), + b=np.array([5, 0, 0])[:mesh.ambient_dim]) mesh3 = merge_disjoint_meshes((mesh2, mesh)) if visualize: from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) - discr = Discretization(cl_ctx, mesh3, - PolynomialWarpAndBlendGroupFactory(3)) + discr = Discretization(cl_ctx, mesh3, discr_grp_factory) from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, 1) + vis = make_visualizer(queue, discr, 3, element_shrink_factor=0.8) vis.write_vtk_file("merged.vtu", []) # }}} @@ -595,9 +612,9 @@ def test_sanity_qhull_nd(ctx_getter, dim, order): from meshmode.discretization import Discretization low_discr = Discretization(ctx, mesh, - PolynomialEquidistantGroupFactory(order)) + PolynomialEquidistantSimplexGroupFactory(order)) high_discr = Discretization(ctx, mesh, - PolynomialEquidistantGroupFactory(order+1)) + PolynomialEquidistantSimplexGroupFactory(order+1)) from meshmode.discretization.connection import make_same_mesh_connection cnx = make_same_mesh_connection(high_discr, low_discr) diff --git a/test/test_refinement.py b/test/test_refinement.py index 8098766d446874b08e676d304d5ac73b88e1163a..558b14088dc012a001cd736d5bc2a0fd854862ed 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -38,7 +38,7 @@ from meshmode.mesh.refinement import Refiner from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, - PolynomialEquidistantGroupFactory, + PolynomialEquidistantSimplexGroupFactory, ) import logging @@ -144,7 +144,7 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, - PolynomialEquidistantGroupFactory + PolynomialEquidistantSimplexGroupFactory ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ ("circle", 1, [20, 30, 40]),