diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml new file mode 100644 index 0000000000000000000000000000000000000000..609ef1425d0a0df37ace94b784f8ba346efe7b15 --- /dev/null +++ b/.gitlab-ci.yml @@ -0,0 +1,36 @@ +Python 2.7 AMD CPU: + script: + - export PY_EXE=python2.7 + - export PYOPENCL_TEST=amd:pu + - export EXTRA_INSTALL="numpy mako" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh + - ". ./build-and-test-py-project.sh" + tags: + - python2.7 + - amd-cl-cpu + except: + - tags +Python 2.7 POCL: + script: + - export PY_EXE=python2.7 + - export PYOPENCL_TEST=portable + - export EXTRA_INSTALL="numpy mako" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh + - ". ./build-and-test-py-project.sh" + tags: + - python2.7 + - pocl + except: + - tags +Python 3.4 POCL: + script: + - export PY_EXE=python3.4 + - export PYOPENCL_TEST=portable + - export EXTRA_INSTALL="numpy mako" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh + - ". ./build-and-test-py-project.sh" + tags: + - python3.4 + - pocl + except: + - tags diff --git a/README.rst b/README.rst index 929ac20b35e2b6103772839755b69b76f44f1e3f..3ca2d6c34f44e78a074468ac94b3547b3c603742 100644 --- a/README.rst +++ b/README.rst @@ -1,4 +1,7 @@ meshmode ======== +* `Source code on Github `_ +* `Documentation `_ + .. TODO diff --git a/doc/_static/akdoc.css b/doc/_static/akdoc.css deleted file mode 100644 index d8b61e3ff7a358e5d5c0f132b5040ec7c3f43e42..0000000000000000000000000000000000000000 --- a/doc/_static/akdoc.css +++ /dev/null @@ -1,59 +0,0 @@ -pre { - line-height: 110%; -} - -.footer { - background-color: #eee; -} - -body > div.container { - margin-top:10px; -} - -dd { - margin-left: 40px; -} - -tt.descname { - font-size: 100%; -} - -code { - color: rgb(51,51,51); -} - -h1 { - padding-bottom:7px; - border-bottom: 1px solid #ccc; -} - -h2 { - padding-bottom:5px; - border-bottom: 1px solid #ccc; -} - -h3 { - padding-bottom:5px; - border-bottom: 1px solid #ccc; -} - -.rubric { - font-size: 120%; - padding-bottom:1px; - border-bottom: 1px solid #ccc; -} - -.headerlink { - padding-left: 1ex; - padding-right: 1ex; -} - -a.headerlink:hover { - text-decoration: none; -} - -blockquote p { - font-size: 100%; - font-weight: normal; - line-height: normal; -}; diff --git a/doc/_templates/layout.html b/doc/_templates/layout.html deleted file mode 100644 index 400e7ec1d49677537aff6bf744e2803ef5c01e9e..0000000000000000000000000000000000000000 --- a/doc/_templates/layout.html +++ /dev/null @@ -1,2 +0,0 @@ -{% extends "!layout.html" %} -{% set bootswatch_css_custom = ['_static/akdoc.css']%} diff --git a/doc/conf.py b/doc/conf.py index f43c6d3596d715838132b0dfc6bfbe6d29eaa386..1e5a0227f1671b4097dfd3b92417b3c7729ac195 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -103,26 +103,24 @@ pygments_style = 'sphinx' # -- Options for HTML output ---------------------------------------------- -try: - import sphinx_bootstrap_theme -except: - from warnings import warn - warn("I would like to use the sphinx bootstrap theme, but can't find it.\n" - "'pip install sphinx_bootstrap_theme' to fix.") -else: - # Activate the theme. - html_theme = 'bootstrap' - html_theme_path = sphinx_bootstrap_theme.get_html_theme_path() - - # Theme options are theme-specific and customize the look and feel of a theme - # further. For a list of options available for each theme, see the - # documentation. - html_theme_options = { - "navbar_fixed_top": "true", - "navbar_site_name": "Contents", - 'bootstrap_version': '3', - 'source_link_position': 'footer', +html_theme = "alabaster" + +html_theme_options = { + "extra_nav_links": { + "🚀 Github": "https://github.com/inducer/meshmode", + "💾 Download Releases": "https://pypi.python.org/pypi/meshmode", } + } + +html_sidebars = { + '**': [ + 'about.html', + 'navigation.html', + 'relations.html', + 'searchbox.html', + ] +} + # Add any paths that contain custom themes here, relative to this directory. #html_theme_path = [] diff --git a/examples/multiple-meshes.py b/examples/multiple-meshes.py new file mode 100644 index 0000000000000000000000000000000000000000..84918e1779a34932ed6269a8884e30e63bb846d0 --- /dev/null +++ b/examples/multiple-meshes.py @@ -0,0 +1,26 @@ +from __future__ import division + +import numpy as np # noqa + + +order = 4 + + +def main(): + from meshmode.mesh.generation import ( # noqa + make_curve_mesh, starfish) + mesh1 = make_curve_mesh(starfish, np.linspace(0, 1, 20), 4) + + from meshmode.mesh.processing import affine_map, merge_disjoint_meshes + mesh2 = affine_map(mesh1, b=np.array([2, 3])) + + mesh = merge_disjoint_meshes((mesh1, mesh2)) + + from meshmode.mesh.visualization import draw_2d_mesh + draw_2d_mesh(mesh, set_bounding_box=True) + + import matplotlib.pyplot as pt + pt.show() + +if __name__ == "__main__": + main() diff --git a/meshmode/__init__.py b/meshmode/__init__.py index e1a4c5e6aedbffd76ceed3f426ee347380d13cce..1c38d87894841834007489ac2a4eea5175ba1428 100644 --- a/meshmode/__init__.py +++ b/meshmode/__init__.py @@ -25,7 +25,7 @@ THE SOFTWARE. __doc__ = """ .. exception:: Error -.. exception:: ConnectivityUnavailable +.. exception:: DataUnavailable """ @@ -33,5 +33,5 @@ class Error(RuntimeError): pass -class ConnectivityUnavailable(Error): +class DataUnavailable(Error): pass diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py index 6088753222fcdfde70e51b0cca237f2af2dd9f78..e1ed150f132a516d409757567d4b5219db669235 100644 --- a/meshmode/discretization/__init__.py +++ b/meshmode/discretization/__init__.py @@ -23,9 +23,10 @@ THE SOFTWARE. """ import numpy as np -from pytools import memoize_method, memoize_method_nested +from pytools import memoize_method, memoize_in import loopy as lp import pyopencl as cl +import pyopencl.array # noqa __doc__ = """ .. autoclass:: ElementGroupBase @@ -37,11 +38,50 @@ __doc__ = """ # {{{ element group base class ElementGroupBase(object): - """Container for the :class:`QBXGrid` data corresponding to + """Container for the :class:`Discretization` data corresponding to one :class:`meshmode.mesh.MeshElementGroup`. .. attribute :: mesh_el_group + .. attribute :: order .. attribute :: node_nr_base + + .. autoattribute:: nelements + .. autoattribute:: nunit_nodes + .. autoattribute:: nnodes + .. autoattribute:: dim + .. automethod:: view + + .. method:: unit_nodes() + + Returns a :class:`numpy.ndarray` of shape ``(dim, nunit_nodes)`` + of reference coordinates of interpolation nodes. + + .. method:: basis() + + Returns a :class:`list` of basis functions that take arrays + of shape ``(dim, n)`` and return an array of shape (n,)`` + (which performs evaluation of the basis function). + + .. method:: grad_basis() + + :returns: a :class:`tuple` of functions, each of which + accepts arrays of shape *(dims, npts)* + and returns a :class:`tuple` of length *dims* containing + the derivatives along each axis as an array of size *npts*. + 'Scalar' evaluation, by passing just one vector of length *dims*, + is also supported. + + .. method:: diff_matrices() + + Return a :attr:`dim`-long :class:`tuple` of matrices of + shape ``(nunit_nodes, nunit_nodes)``, each of which, + when applied to an array of nodal values, take derivatives + in the reference (r,s,t) directions. + + .. method:: weights() + + Returns an array of length :attr:`nunit_nodes` containing + quadrature weights. """ def __init__(self, mesh_el_group, order, node_nr_base): @@ -79,6 +119,11 @@ class ElementGroupBase(object): (-1, -1)) def view(self, global_array): + """Return a view of *global_array* of shape ``(..., nelements, + nunit_nodes)`` where *global_array* is of shape ``(..., nnodes)``, + where *nnodes* is the global (per-discretization) node count. + """ + return global_array[ ..., self.node_nr_base:self.node_nr_base + self.nnodes] \ .reshape( @@ -127,7 +172,9 @@ class Discretization(object): .. attribute :: groups - .. method:: empty(dtype, queue=None, extra_dims=None) + .. automethod:: empty + + .. automethod:: zeros .. method:: nodes() @@ -141,11 +188,6 @@ class Discretization(object): """ def __init__(self, cl_ctx, mesh, group_factory, real_dtype=np.float64): - """ - :arg order: A polynomial-order-like parameter passed unmodified to - :attr:`group_class`. See subclasses for more precise definition. - """ - self.cl_context = cl_ctx self.mesh = mesh @@ -170,7 +212,20 @@ class Discretization(object): def ambient_dim(self): return self.mesh.ambient_dim - def empty(self, dtype, queue=None, extra_dims=None): + def empty(self, queue=None, dtype=None, extra_dims=None, allocator=None): + """Return an empty DOF vector. + + :arg dtype: type special value 'c' will result in a + vector of dtype :attr:`self.complex_dtype`. If + *None* (the default), a real vector will be returned. + """ + if dtype is None: + dtype = self.real_dtype + elif dtype == "c": + dtype = self.complex_dtype + else: + dtype = np.dtype(dtype) + if queue is None: first_arg = self.cl_context else: @@ -180,11 +235,14 @@ class Discretization(object): if extra_dims is not None: shape = extra_dims + shape - return cl.array.empty(first_arg, shape, dtype=dtype) + return cl.array.empty(first_arg, shape, dtype=dtype, allocator=allocator) + + def zeros(self, queue, dtype=None, extra_dims=None, allocator=None): + return self.empty(queue, dtype=dtype, extra_dims=extra_dims, + allocator=allocator).fill(0) - def num_reference_derivative( - self, queue, ref_axes, vec): - @memoize_method_nested + def num_reference_derivative(self, queue, ref_axes, vec): + @memoize_in(self, "reference_derivative_knl") def knl(): knl = lp.make_kernel( """{[k,i,j]: @@ -196,9 +254,12 @@ class Discretization(object): knl = lp.split_iname(knl, "i", 16, inner_tag="l.0") return lp.tag_inames(knl, dict(k="g.0")) - result = self.empty(vec.dtype) + result = self.empty(dtype=vec.dtype) for grp in self.groups: + if grp.nelements == 0: + continue + mat = None for ref_axis in ref_axes: next_mat = grp.diff_matrices()[ref_axis] @@ -212,24 +273,28 @@ class Discretization(object): return result def quad_weights(self, queue): - @memoize_method_nested + @memoize_in(self, "quad_weights_knl") def knl(): knl = lp.make_kernel( "{[k,i]: 0<=kdei", - resampling_mat, - mgrp.nodes[:, batch_boundary_el_numbers_in_grp, :]) - - # }}} - - 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, - ) - - batch_base += len(batch_boundary_el_numbers_in_grp) - - bdry_mesh_group = SimplexElementGroup( - mgrp.order, vertex_indices, nodes, unit_nodes=bdry_unit_nodes) - bdry_mesh_groups.append(bdry_mesh_group) - - bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups) - - from meshmode.discretization import Discretization - bdry_discr = Discretization( - discr.cl_context, bdry_mesh, group_factory) - - connection = _build_boundary_connection( - queue, discr, bdry_discr, connection_data) - - logger.info("building boundary connection: done") - - return bdry_mesh, bdry_discr, connection - -# }}} - - -# {{{ refinement connection - -def make_refinement_connection(refiner, coarse_discr): - pass - -# }}} - -# vim: foldmethod=marker diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..269ed5d14c61e3b33ef3d7e6d4354b105e9b5f50 --- /dev/null +++ b/meshmode/discretization/connection/__init__.py @@ -0,0 +1,444 @@ +from __future__ import division, print_function, absolute_import + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from six.moves import range, zip + +import numpy as np +import pyopencl as cl +import pyopencl.array # noqa +from pytools import memoize_method, memoize_in + +from meshmode.discretization.connection.same_mesh import \ + make_same_mesh_connection +from meshmode.discretization.connection.face import ( + FRESTR_INTERIOR_FACES, FRESTR_ALL_FACES, + make_face_restriction, make_face_to_all_faces_embedding) +from meshmode.discretization.connection.opposite_face import \ + make_opposite_face_connection + +import logging +logger = logging.getLogger(__name__) + + +__all__ = [ + "DiscretizationConnection", + "make_same_mesh_connection", + "FRESTR_INTERIOR_FACES", "FRESTR_ALL_FACES", + "make_face_restriction", + "make_face_to_all_faces_embedding", + "make_opposite_face_connection" + ] + +__doc__ = """ +.. autoclass:: DiscretizationConnection + +.. autofunction:: make_same_mesh_connection + +.. autofunction:: FRESTR_INTERIOR_FACES +.. autofunction:: FRESTR_ALL_FACES +.. autofunction:: make_face_restriction +.. autofunction:: make_face_to_all_faces_embedding + +.. autofunction:: make_opposite_face_connection + +Implementation details +^^^^^^^^^^^^^^^^^^^^^^ + +.. autoclass:: InterpolationBatch + +.. autoclass:: DiscretizationConnectionElementGroup +""" + + +# {{{ interpolation batch + +class InterpolationBatch(object): + """One interpolation batch captures how a batch of elements *within* an + element group should be an interpolated. Note that while it's possible that + an interpolation batch takes care of interpolating an entire element group + from source to target, that's not *necessarily* the case. Consider the case + of extracting boundary values of a discretization. For, say, a triangle, at + least three different interpolation batches are needed to cover boundary + edges that fall onto each of the three edges of the unit triangle. + + .. attribute:: from_group_index + + An integer indicating from which element group in the *from* discretization + the data should be interpolated. + + .. attribute:: from_element_indices + + ``element_id_t [nelements]``. (a :class:`pyopencl.array.Array`) + This contains the (group-local) element index (relative to + :attr:`from_group_index` from which this "*to*" element's data will be + interpolated. + + .. attribute:: to_element_indices + + ``element_id_t [nelements]``. (a :class:`pyopencl.array.Array`) + This contains the (group-local) element index to which this "*to*" + element's data will be interpolated. + + .. attribute:: result_unit_nodes + + A :class:`numpy.ndarray` of shape + ``(from_group.dim,to_group.nelements,to_group.nunit_nodes)`` + storing the coordinates of the nodes (in unit coordinates + of the *from* reference element) from which the node + locations of this element should be interpolated. + + .. autoattribute:: nelements + + .. attribute:: to_element_face + + *int* or *None*. (a :class:`pyopencl.array.Array` if existent) If this + interpolation batch targets interpolation *to* a face, then this number + captures the face number (on all elements referenced by + :attr:`from_element_indices` to which this batch interpolates. (Since + there is a fixed set of "from" unit nodes per batch, one batch will + always go to a single face index.) + """ + + def __init__(self, from_group_index, from_element_indices, + to_element_indices, result_unit_nodes, to_element_face): + self.from_group_index = from_group_index + self.from_element_indices = from_element_indices + self.to_element_indices = to_element_indices + self.result_unit_nodes = result_unit_nodes + self.to_element_face = to_element_face + + @property + def nelements(self): + return len(self.from_element_indices) + +# }}} + + +# {{{ connection element group + +class DiscretizationConnectionElementGroup(object): + """ + .. attribute:: batches + + A list of :class:`InterpolationBatch` instances. + """ + def __init__(self, batches): + self.batches = batches + +# }}} + + +# {{{ connection class + +class DiscretizationConnection(object): + """Data supporting an interpolation-like operation that takes in data on + one discretization and returns it on another. Implemented applications + include: + + * upsampling/downsampling on the same mesh + * restricition to the boundary + * interpolation to a refined/coarsened mesh + * interpolation onto opposing faces + + .. attribute:: from_discr + + .. attribute:: to_discr + + .. attribute:: groups + + a list of :class:`DiscretizationConnectionElementGroup` + instances, with a one-to-one correspondence to the groups in + :attr:`to_discr`. + + .. attribute:: is_surjective + + A :class:`bool` indicating whether every output degree + of freedom is set by the connection. + + .. automethod:: __call__ + + .. automethod:: full_resample_matrix + + """ + + def __init__(self, from_discr, to_discr, groups, is_surjective): + if from_discr.cl_context != to_discr.cl_context: + raise ValueError("from_discr and to_discr must live in the " + "same OpenCL context") + + self.cl_context = from_discr.cl_context + + if from_discr.mesh.vertex_id_dtype != to_discr.mesh.vertex_id_dtype: + raise ValueError("from_discr and to_discr must agree on the " + "vertex_id_dtype") + + if from_discr.mesh.element_id_dtype != to_discr.mesh.element_id_dtype: + raise ValueError("from_discr and to_discr must agree on the " + "element_id_dtype") + + self.cl_context = from_discr.cl_context + + self.from_discr = from_discr + self.to_discr = to_discr + self.groups = groups + + self.is_surjective = is_surjective + + @memoize_method + def _resample_matrix(self, to_group_index, ibatch_index): + import modepy as mp + ibatch = self.groups[to_group_index].batches[ibatch_index] + from_grp = self.from_discr.groups[ibatch.from_group_index] + + result = mp.resampling_matrix( + from_grp.basis(), + ibatch.result_unit_nodes, from_grp.unit_nodes) + + with cl.CommandQueue(self.cl_context) as queue: + return cl.array.to_device(queue, result).with_queue(None) + + @memoize_method + def _resample_point_pick_indices(self, to_group_index, ibatch_index, + tol_multiplier=None): + """If :meth:`_resample_matrix` *R* is a row subset of a permutation matrix *P*, + return the index subset I so that, loosely, ``x[I] == R @ x``. + + Will return *None* if no such index array exists, or a + :class:`pyopencl.array.Array` containing the index subset. + """ + + with cl.CommandQueue(self.cl_context) as queue: + mat = self._resample_matrix(to_group_index, ibatch_index).get( + queue=queue) + + nrows, ncols = mat.shape + result = np.zeros(nrows, dtype=self.to_discr.mesh.element_id_dtype) + + if tol_multiplier is None: + tol_multiplier = 50 + + tol = np.finfo(mat.dtype).eps * tol_multiplier + + for irow in range(nrows): + one_indices, = np.where(np.abs(mat[irow] - 1) < tol) + zero_indices, = np.where(np.abs(mat[irow]) < tol) + + if len(one_indices) != 1: + return None + if len(zero_indices) != ncols - 1: + return None + + one_index, = one_indices + result[irow] = one_index + + with cl.CommandQueue(self.cl_context) as queue: + return cl.array.to_device(queue, result).with_queue(None) + + def full_resample_matrix(self, queue): + """Build a dense matrix representing this discretization connection. + + .. warning:: + + On average, this will be exceedingly expensive (:math:`O(N^2)` in + the number *N* of discretization points) in terms of memory usage + and thus not what you'd typically want. + """ + + @memoize_in(self, "oversample_mat_knl") + def knl(): + import loopy as lp + knl = lp.make_kernel( + """{[k,i,j]: + 0<=k not looking for those + continue + + group_boundary_faces.extend( + zip(fagrp.elements, fagrp.element_faces)) + + elif boundary_tag is FRESTR_ALL_FACES: + group_boundary_faces.extend( + (iel, iface) + for iface in range(grp.mesh_el_group.nfaces) + for iel in range(grp.nelements) + ) + + else: + bdry_grp = fagrp_map.get(None) + if bdry_grp is not None: + nb_el_bits = -bdry_grp.neighbors + face_relevant_flags = (nb_el_bits & btag_bit) != 0 + + group_boundary_faces.extend( + zip( + bdry_grp.elements[face_relevant_flags], + bdry_grp.element_faces[face_relevant_flags])) + + # }}} + + grp_face_vertex_indices = mgrp.face_vertex_indices() + grp_vertex_unit_coordinates = mgrp.vertex_unit_coordinates() + + batch_base = 0 + + # group by face_id + + 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) + + # {{{ 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: + ngroup_bdry_elements = nbatch_elements + else: + ngroup_bdry_elements = len(group_boundary_faces) + + vertex_indices = np.empty( + (ngroup_bdry_elements, mgrp.dim+1-1), + 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), + dtype=np.float64) + + # }}} + + new_el_numbers = batch_base + np.arange(nbatch_elements) + if not per_face_groups: + batch_base += nbatch_elements + + # {{{ 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 A, b such that A [e_1 e_2] + b = [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_vertex_unit_coordinates[1:] + - face_vertex_unit_coordinates[0]).T + + face_unit_nodes = (np.dot(A, bdry_unit_nodes_01).T + b).T + + resampling_mat = mp.resampling_matrix( + vol_basis, + face_unit_nodes, mgrp.unit_nodes) + + # }}} + + # {{{ build information for mesh element group + + # Find vertex_indices + glob_face_vertices = mgrp.vertex_indices[ + batch_boundary_el_numbers_in_grp][:, loc_face_vertices] + vertex_indices[new_el_numbers] = \ + vol_to_bdry_vertices[glob_face_vertices] + + # Find nodes + nodes[:, new_el_numbers, :] = np.einsum( + "ij,dej->dei", + resampling_mat, + mgrp.nodes[:, batch_boundary_el_numbers_in_grp, :]) + + # }}} + + 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, + ) + + is_last_face = face_id + 1 == mgrp.nfaces + + if per_face_groups or is_last_face: + bdry_mesh_group = SimplexElementGroup( + mgrp.order, vertex_indices, nodes, + unit_nodes=bdry_unit_nodes) + bdry_mesh_groups.append(bdry_mesh_group) + + bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups) + + from meshmode.discretization import Discretization + bdry_discr = Discretization( + discr.cl_context, bdry_mesh, group_factory) + + with cl.CommandQueue(discr.cl_context) as queue: + connection = _build_boundary_connection( + queue, discr, bdry_discr, connection_data, + per_face_groups) + + logger.info("building face restriction: done") + + return connection + +# }}} + + +# {{{ face -> all_faces connection + +def make_face_to_all_faces_embedding(faces_connection, all_faces_discr): + """Return a + :class:`meshmode.discretization.connection.DiscretizationConnection` + connecting a discretization containing some faces of a discretization + to one containing all faces. + + :arg faces_connection: must be the (connection) result of calling + :func:`meshmode.discretization.connection.make_face_restriction` + with + :class:`meshmode.discretization.connection.FRESTR_INTERIOR_FACES` + or a boundary tag. + :arg all_faces_discr: must be the (discretization) result of calling + :func:`meshmode.discretization.connection.make_face_restriction` + with + :class:`meshmode.discretization.connection.FRESTR_ALL_FACES` + for the same volume discretization as the one from which + *faces_discr* was obtained. + """ + + vol_discr = faces_connection.from_discr + faces_discr = faces_connection.to_discr + + per_face_groups = ( + len(vol_discr.groups) != len(faces_discr.groups)) + + if len(faces_discr.groups) != len(all_faces_discr.groups): + raise ValueError("faces_discr and all_faces_discr must have the " + "same number of groups") + if len(faces_connection.groups) != len(all_faces_discr.groups): + raise ValueError("faces_connection and all_faces_discr must have the " + "same number of groups") + + from meshmode.discretization.connection import ( + DiscretizationConnection, + DiscretizationConnectionElementGroup, + InterpolationBatch) + + i_faces_grp = 0 + + with cl.CommandQueue(vol_discr.cl_context) as queue: + groups = [] + for ivol_grp, vol_grp in enumerate(vol_discr.groups): + batches = [] + + nfaces = vol_grp.mesh_el_group.nfaces + for iface in range(nfaces): + faces_grp = faces_discr.groups[i_faces_grp] + all_faces_grp = all_faces_discr.groups[i_faces_grp] + + if per_face_groups: + assert len(faces_connection.groups[i_faces_grp].batches) == 1 + else: + assert (len(faces_connection.groups[i_faces_grp].batches) + == nfaces) + + assert np.array_equal( + faces_grp.unit_nodes, all_faces_grp.unit_nodes) + + # {{{ find src_batch + + src_batches = faces_connection.groups[i_faces_grp].batches + if per_face_groups: + src_batch, = src_batches + else: + src_batch = src_batches[iface] + del src_batches + + # }}} + + if per_face_groups: + to_element_indices = src_batch.from_element_indices + else: + assert all_faces_grp.nelements == nfaces * vol_grp.nelements + + to_element_indices = ( + vol_grp.nelements*iface + + src_batch.from_element_indices.with_queue(queue) + ).with_queue(None) + + batches.append( + InterpolationBatch( + from_group_index=i_faces_grp, + from_element_indices=src_batch.to_element_indices, + to_element_indices=to_element_indices, + result_unit_nodes=all_faces_grp.unit_nodes, + to_element_face=None)) + + is_last_face = iface + 1 == nfaces + if per_face_groups or is_last_face: + groups.append( + DiscretizationConnectionElementGroup(batches=batches)) + batches = [] + + i_faces_grp += 1 + + return DiscretizationConnection( + faces_connection.to_discr, + all_faces_discr, + groups, + is_surjective=False) + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py new file mode 100644 index 0000000000000000000000000000000000000000..70728522b2c06f119a5b20598616e2c464eea249 --- /dev/null +++ b/meshmode/discretization/connection/opposite_face.py @@ -0,0 +1,395 @@ +from __future__ import division, print_function, absolute_import + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +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 +import numpy.linalg as la +import pyopencl as cl +import pyopencl.array # noqa + +import logging +logger = logging.getLogger(__name__) + + +# {{{ opposite-face connection + +def _make_cross_face_batches( + queue, vol_discr, bdry_discr, + i_tgt_grp, i_src_grp, + i_face_tgt, + adj_grp, + vbc_tgt_grp_face_batch, src_grp_el_lookup): + + # {{{ index wrangling + + # Assert that the adjacency group and the restriction + # interpolation batch and the adjacency group have the same + # element ordering. + + adj_grp_tgt_flags = adj_grp.element_faces == i_face_tgt + + assert ( + np.array_equal( + adj_grp.elements[adj_grp_tgt_flags], + vbc_tgt_grp_face_batch.from_element_indices + .get(queue=queue))) + + # find to_element_indices + + to_bdry_element_indices = ( + vbc_tgt_grp_face_batch.to_element_indices + .get(queue=queue)) + + # find from_element_indices + + from_vol_element_indices = adj_grp.neighbors[adj_grp_tgt_flags] + from_element_faces = adj_grp.neighbor_faces[adj_grp_tgt_flags] + + from_bdry_element_indices = src_grp_el_lookup[ + from_vol_element_indices, from_element_faces] + + # }}} + + # {{{ visualization (for debugging) + + if 0: + print("TVE", adj_grp.elements[adj_grp_tgt_flags]) + print("TBE", to_bdry_element_indices) + print("FVE", from_vol_element_indices) + from meshmode.mesh.visualization import draw_2d_mesh + import matplotlib.pyplot as pt + draw_2d_mesh(vol_discr.mesh, draw_element_numbers=True, + set_bounding_box=True, + draw_vertex_numbers=False, + draw_face_numbers=True, + fill=None) + pt.figure() + + draw_2d_mesh(bdry_discr.mesh, draw_element_numbers=True, + set_bounding_box=True, + draw_vertex_numbers=False, + draw_face_numbers=True, + fill=None) + + pt.show() + # }}} + + # {{{ invert face map (using Gauss-Newton) + + to_bdry_nodes = ( + # FIXME: This should view-then-transfer (but PyOpenCL doesn't do + # non-contiguous transfers for now). + bdry_discr.groups[i_tgt_grp].view( + bdry_discr.nodes().get(queue=queue)) + [:, to_bdry_element_indices]) + + tol = 1e4 * np.finfo(to_bdry_nodes.dtype).eps + + from_mesh_grp = bdry_discr.mesh.groups[i_src_grp] + from_grp = bdry_discr.groups[i_src_grp] + + dim = from_grp.dim + ambient_dim, nelements, nto_unit_nodes = to_bdry_nodes.shape + + initial_guess = np.mean(from_mesh_grp.vertex_unit_coordinates(), axis=0) + from_unit_nodes = np.empty((dim, nelements, nto_unit_nodes)) + from_unit_nodes[:] = initial_guess.reshape(-1, 1, 1) + + import modepy as mp + from_vdm = mp.vandermonde(from_grp.basis(), from_grp.unit_nodes) + from_inv_t_vdm = la.inv(from_vdm.T) + from_nfuncs = len(from_grp.basis()) + + # (ambient_dim, nelements, nfrom_unit_nodes) + from_bdry_nodes = ( + # FIXME: This should view-then-transfer (but PyOpenCL doesn't do + # non-contiguous transfers for now). + bdry_discr.groups[i_src_grp].view( + bdry_discr.nodes().get(queue=queue)) + [:, from_bdry_element_indices]) + + def apply_map(unit_nodes): + # unit_nodes: (dim, nelements, nto_unit_nodes) + + # basis_at_unit_nodes + basis_at_unit_nodes = np.empty((from_nfuncs, nelements, nto_unit_nodes)) + + for i, f in enumerate(from_grp.basis()): + basis_at_unit_nodes[i] = ( + f(unit_nodes.reshape(dim, -1)) + .reshape(nelements, nto_unit_nodes)) + + intp_coeffs = np.einsum("fj,jet->fet", from_inv_t_vdm, basis_at_unit_nodes) + + # If we're interpolating 1, we had better get 1 back. + one_deviation = np.abs(np.sum(intp_coeffs, axis=0) - 1) + assert (one_deviation < tol).all(), np.max(one_deviation) + + return np.einsum("fet,aef->aet", intp_coeffs, from_bdry_nodes) + + def get_map_jacobian(unit_nodes): + # unit_nodes: (dim, nelements, nto_unit_nodes) + + # basis_at_unit_nodes + dbasis_at_unit_nodes = np.empty( + (dim, from_nfuncs, nelements, nto_unit_nodes)) + + for i, df in enumerate(from_grp.grad_basis()): + df_result = df(unit_nodes.reshape(dim, -1)) + + for rst_axis, df_r in enumerate(df_result): + dbasis_at_unit_nodes[rst_axis, i] = ( + df_r.reshape(nelements, nto_unit_nodes)) + + dintp_coeffs = np.einsum( + "fj,rjet->rfet", from_inv_t_vdm, dbasis_at_unit_nodes) + + return np.einsum("rfet,aef->raet", dintp_coeffs, from_bdry_nodes) + + # {{{ test map applier and jacobian + + if 0: + u = from_unit_nodes + f = apply_map(u) + for h in [1e-1, 1e-2]: + du = h*np.random.randn(*u.shape) + + f_2 = apply_map(u+du) + + jf = get_map_jacobian(u) + + f2_2 = f + np.einsum("raet,ret->aet", jf, du) + + print(h, la.norm((f_2-f2_2).ravel())) + + # }}} + + # {{{ visualize initial guess + + if 0: + import matplotlib.pyplot as pt + guess = apply_map(from_unit_nodes) + goals = to_bdry_nodes + + from meshmode.discretization.visualization import draw_curve + draw_curve(bdry_discr) + + pt.plot(guess[0].reshape(-1), guess[1].reshape(-1), "or") + pt.plot(goals[0].reshape(-1), goals[1].reshape(-1), "og") + pt.plot(from_bdry_nodes[0].reshape(-1), from_bdry_nodes[1].reshape(-1), "o", + color="purple") + pt.show() + + # }}} + + logger.info("make_opposite_face_connection: begin gauss-newton") + + niter = 0 + while True: + resid = apply_map(from_unit_nodes) - to_bdry_nodes + + df = get_map_jacobian(from_unit_nodes) + df_inv_resid = np.empty_like(from_unit_nodes) + + # For the 1D/2D accelerated versions, we'll use the normal + # equations and Cramer's rule. If you're looking for high-end + # numerics, look no further than meshmode. + + if dim == 1: + # A is df.T + ata = np.einsum("iket,jket->ijet", df, df) + atb = np.einsum("iket,ket->iet", df, resid) + + df_inv_resid = atb / ata[0, 0] + + elif dim == 2: + # A is df.T + ata = np.einsum("iket,jket->ijet", df, df) + atb = np.einsum("iket,ket->iet", df, resid) + + det = ata[0, 0]*ata[1, 1] - ata[0, 1]*ata[1, 0] + + df_inv_resid = np.empty_like(from_unit_nodes) + df_inv_resid[0] = 1/det * (ata[1, 1] * atb[0] - ata[1, 0]*atb[1]) + df_inv_resid[1] = 1/det * (-ata[0, 1] * atb[0] + ata[0, 0]*atb[1]) + + else: + # The boundary of a 3D mesh is 2D, so that's the + # highest-dimensional case we genuinely care about. + # + # This stinks, performance-wise, because it's not vectorized. + # But we'll only hit it for boundaries of 4+D meshes, in which + # case... good luck. :) + for e in range(nelements): + for t in range(nto_unit_nodes): + df_inv_resid[:, e, t], _, _, _ = \ + la.lstsq(df[:, :, e, t].T, resid[:, e, t]) + + from_unit_nodes = from_unit_nodes - df_inv_resid + + max_resid = np.max(np.abs(resid)) + logger.debug("gauss-newton residual: %g" % max_resid) + + if max_resid < tol: + logger.info("make_opposite_face_connection: gauss-newton: done, " + "final residual: %g" % max_resid) + break + + niter += 1 + if niter > 10: + raise RuntimeError("Gauss-Newton (for finding opposite-face reference " + "coordinates) did not converge") + + # }}} + + # {{{ find groups of from_unit_nodes + + def to_dev(ary): + return cl.array.to_device(queue, ary, array_queue=None) + + done_elements = np.zeros(nelements, dtype=np.bool) + while True: + todo_elements, = np.where(~done_elements) + if not len(todo_elements): + return + + template_unit_nodes = from_unit_nodes[:, todo_elements[0], :] + + unit_node_dist = np.max(np.max(np.abs( + from_unit_nodes[:, todo_elements, :] + - + template_unit_nodes.reshape(dim, 1, -1)), + axis=2), axis=0) + + close_els = todo_elements[unit_node_dist < tol] + done_elements[close_els] = True + + unit_node_dist = np.max(np.max(np.abs( + from_unit_nodes[:, todo_elements, :] + - + template_unit_nodes.reshape(dim, 1, -1)), + axis=2), axis=0) + + from meshmode.discretization.connection import InterpolationBatch + yield InterpolationBatch( + from_group_index=i_src_grp, + from_element_indices=to_dev(from_bdry_element_indices[close_els]), + to_element_indices=to_dev(to_bdry_element_indices[close_els]), + result_unit_nodes=template_unit_nodes, + to_element_face=None) + + # }}} + + +def _find_ibatch_for_face(vbc_tgt_grp_batches, iface): + vbc_tgt_grp_face_batches = [ + batch + for batch in vbc_tgt_grp_batches + if batch.to_element_face == iface] + + assert len(vbc_tgt_grp_face_batches) == 1 + + vbc_tgt_grp_face_batch, = vbc_tgt_grp_face_batches + + return vbc_tgt_grp_face_batch + + +def _make_el_lookup_table(queue, connection, igrp): + from_nelements = connection.from_discr.groups[igrp].nelements + from_nfaces = connection.from_discr.mesh.groups[igrp].nfaces + + iel_lookup = np.empty((from_nelements, from_nfaces), + dtype=connection.from_discr.mesh.element_id_dtype) + iel_lookup.fill(-1) + + for ibatch, batch in enumerate(connection.groups[igrp].batches): + from_element_indices = batch.from_element_indices.get(queue=queue) + iel_lookup[from_element_indices, batch.to_element_face] = \ + batch.to_element_indices.get(queue=queue) + + return iel_lookup + + +def make_opposite_face_connection(volume_to_bdry_conn): + """Given a boundary restriction connection *volume_to_bdry_conn*, + return a :class:`DiscretizationConnection` that performs data + exchange across opposite faces. + """ + + vol_discr = volume_to_bdry_conn.from_discr + vol_mesh = vol_discr.mesh + bdry_discr = volume_to_bdry_conn.to_discr + + # make sure we were handed a volume-to-boundary connection + for i_tgrp, conn_grp in enumerate(volume_to_bdry_conn.groups): + for batch in conn_grp.batches: + assert batch.from_group_index == i_tgrp + assert batch.to_element_face is not None + + ngrps = len(volume_to_bdry_conn.groups) + assert ngrps == len(vol_discr.groups) + assert ngrps == len(bdry_discr.groups) + + # One interpolation batch in this connection corresponds + # to a key (i_tgt_grp,) (i_src_grp, i_face_tgt,) + + with cl.CommandQueue(vol_discr.cl_context) as queue: + # a list of batches for each group + groups = [[] for i_tgt_grp in range(ngrps)] + + for i_src_grp in range(ngrps): + src_grp_el_lookup = _make_el_lookup_table( + queue, volume_to_bdry_conn, i_src_grp) + + for i_tgt_grp in range(ngrps): + vbc_tgt_grp_batches = volume_to_bdry_conn.groups[i_tgt_grp].batches + + adj_grp = vol_mesh.facial_adjacency_groups[i_tgt_grp][i_src_grp] + + for i_face_tgt in range(vol_mesh.groups[i_tgt_grp].nfaces): + vbc_tgt_grp_face_batch = _find_ibatch_for_face( + vbc_tgt_grp_batches, i_face_tgt) + + groups[i_tgt_grp].extend( + _make_cross_face_batches( + queue, vol_discr, bdry_discr, + i_tgt_grp, i_src_grp, + i_face_tgt, + adj_grp, + vbc_tgt_grp_face_batch, src_grp_el_lookup)) + + from meshmode.discretization.connection import ( + DiscretizationConnection, DiscretizationConnectionElementGroup) + return DiscretizationConnection( + from_discr=bdry_discr, + to_discr=bdry_discr, + groups=[ + DiscretizationConnectionElementGroup(batches=batches) + for batches in groups], + is_surjective=True) + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/discretization/connection/same_mesh.py b/meshmode/discretization/connection/same_mesh.py new file mode 100644 index 0000000000000000000000000000000000000000..3b73e34fefd5115e01b97728ab414cc8c77fd178 --- /dev/null +++ b/meshmode/discretization/connection/same_mesh.py @@ -0,0 +1,65 @@ +from __future__ import division, print_function, absolute_import + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np +import pyopencl as cl +import pyopencl.array # noqa + + +# {{{ same-mesh constructor + +def make_same_mesh_connection(to_discr, from_discr): + from meshmode.discretization.connection import ( + InterpolationBatch, DiscretizationConnectionElementGroup, + DiscretizationConnection) + + if from_discr.mesh is not to_discr.mesh: + raise ValueError("from_discr and to_discr must be based on " + "the same mesh") + + assert to_discr.cl_context == from_discr.cl_context + + with cl.CommandQueue(to_discr.cl_context) as queue: + groups = [] + for igrp, (fgrp, tgrp) in enumerate(zip(from_discr.groups, to_discr.groups)): + all_elements = cl.array.arange(queue, + fgrp.nelements, + dtype=np.intp).with_queue(None) + ibatch = InterpolationBatch( + from_group_index=igrp, + from_element_indices=all_elements, + to_element_indices=all_elements, + result_unit_nodes=tgrp.unit_nodes, + to_element_face=None) + + groups.append( + DiscretizationConnectionElementGroup([ibatch])) + + return DiscretizationConnection( + from_discr, to_discr, groups, + is_surjective=True) + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index d59346691c93f7e967a1174b0a4622fdfff05426..a5319bc84f47ab6689f502629ff6009893033832 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -38,6 +38,7 @@ Group types .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup .. autoclass:: PolynomialWarpAndBlendElementGroup +.. autoclass:: PolynomialEquidistantElementGroup Group factories ^^^^^^^^^^^^^^^ @@ -45,6 +46,7 @@ Group factories .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory .. autoclass:: PolynomialWarpAndBlendGroupFactory +.. autoclass:: PolynomialEquidistantGroupFactory """ # FIXME Most of the loopy kernels will break as soon as we start using multiple @@ -63,13 +65,16 @@ from meshmode.discretization import ElementGroupBase class PolynomialSimplexElementGroupBase(ElementGroupBase): def basis(self): - return mp.simplex_onb(self.dim, self.order) + return mp.simplex_best_available_basis(self.dim, self.order) + + def grad_basis(self): + return mp.grad_simplex_best_available_basis(self.dim, self.order) @memoize_method def from_mesh_interp_matrix(self): meg = self.mesh_el_group return mp.resampling_matrix( - mp.simplex_onb(meg.dim, meg.order), + mp.simplex_best_available_basis(meg.dim, meg.order), self.unit_nodes, meg.unit_nodes) @@ -77,7 +82,7 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase): def diff_matrices(self): result = mp.differentiation_matrices( self.basis(), - mp.grad_simplex_onb(self.dim, self.order), + self.grad_basis(), self.unit_nodes) if not isinstance(result, tuple): @@ -85,13 +90,6 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase): else: return result - @memoize_method - def resampling_matrix(self): - meg = self.mesh_el_group - return mp.resampling_matrix( - mp.simplex_onb(self.dim, meg.order), - self.unit_nodes, meg.unit_nodes) - class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): """Elemental discretization supplying a high-order quadrature rule @@ -154,11 +152,20 @@ class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): return self._quadrature_rule().weights -class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase): +class _MassMatrixQuadratureElementGroup(PolynomialSimplexElementGroupBase): + @property + @memoize_method + def weights(self): + return np.dot( + mp.mass_matrix(self.basis(), self.unit_nodes), + np.ones(len(self.basis()))) + + +class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): """Elemental discretization with a number of nodes matching the number of polynomials in :math:`P^k`, hence usable for differentiation and - interpolation. Interpolation nodes are present on the boundary of the - simplex. + interpolation. Interpolation nodes edge-clustered for avoidance of Runge + phenomena. Nodes are present on the boundary of the simplex. """ @property @memoize_method @@ -170,12 +177,24 @@ class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase): assert dim2 == dim return result + +class PolynomialEquidistantElementGroup(_MassMatrixQuadratureElementGroup): + """Elemental discretization with a number of nodes matching the number of + polynomials in :math:`P^k`, hence usable for differentiation and + interpolation. Interpolation nodes are present on the boundary of the + simplex. + + .. versionadded:: 2016.1 + """ @property @memoize_method - def weights(self): - return np.dot( - mp.mass_matrix(self.basis(), self.unit_nodes), - np.ones(len(self.basis()))) + def unit_nodes(self): + dim = self.mesh_el_group.dim + result = mp.equidistant_nodes(dim, self.order) + + dim2, nunit_nodes = result.shape + assert dim2 == dim + return result # }}} @@ -212,6 +231,15 @@ class PolynomialWarpAndBlendGroupFactory(OrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = PolynomialWarpAndBlendElementGroup + +class PolynomialEquidistantGroupFactory(OrderBasedGroupFactory): + """ + .. versionadded:: 2016.1 + """ + + mesh_group_class = _MeshSimplexElementGroup + group_class = PolynomialEquidistantElementGroup + # }}} diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 3c2565fdb5a7ac206d76b0c2fb5ed44853939e54..d68ec0c58078e68417721c8205fc3f24021d0bad 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -34,7 +34,7 @@ __doc__ = """ .. autoclass:: Visualizer -.. autofunction:: write_mesh_connectivity_vtk_file +.. autofunction:: write_nodal_adjacency_vtk_file """ @@ -77,10 +77,10 @@ class Visualizer(object): .. automethod:: write_vtk_file """ - def __init__(self, discr, vis_discr, connection): - self.discr = discr - self.vis_discr = vis_discr + def __init__(self, connection): self.connection = connection + self.discr = connection.from_discr + self.vis_discr = connection.to_discr def _resample_and_get(self, queue, vec): from pytools.obj_array import with_object_array_or_scalar @@ -240,16 +240,36 @@ def make_visualizer(queue, discr, vis_order): real_dtype=discr.real_dtype) from meshmode.discretization.connection import \ make_same_mesh_connection - cnx = make_same_mesh_connection(queue, vis_discr, discr) - return Visualizer(discr, vis_discr, cnx) + return Visualizer(make_same_mesh_connection(vis_discr, discr)) # }}} -# {{{ connectivity +# {{{ draw_curve -def write_mesh_connectivity_vtk_file(file_name, mesh, compressor=None,): +def draw_curve(discr): + mesh = discr.mesh + + import matplotlib.pyplot as pt + pt.plot(mesh.vertices[0], mesh.vertices[1], "o") + + color = pt.cm.rainbow(np.linspace(0, 1, len(discr.groups))) + with cl.CommandQueue(discr.cl_context) as queue: + for i, group in enumerate(discr.groups): + group_nodes = group.view(discr.nodes()).get(queue=queue) + pt.plot( + group_nodes[0].T, + group_nodes[1].T, "-x", + label="Group %d" % i, + color=color[i]) + +# }}} + + +# {{{ adjacency + +def write_nodal_adjacency_vtk_file(file_name, mesh, compressor=None,): from pyvisfile.vtk import ( UnstructuredGrid, DataArray, AppendedDataXMLGenerator, @@ -266,16 +286,16 @@ def write_mesh_connectivity_vtk_file(file_name, mesh, compressor=None,): np.sum(mesh.vertices[:, grp.vertex_indices], axis=-1) / grp.vertex_indices.shape[-1]) - cnx = mesh.element_connectivity + adj = mesh.nodal_adjacency - nconnections = len(cnx.neighbors) + nconnections = len(adj.neighbors) connections = np.empty((nconnections, 2), dtype=np.int32) - nb_starts = cnx.neighbors_starts + nb_starts = adj.neighbors_starts for iel in range(mesh.nelements): connections[nb_starts[iel]:nb_starts[iel+1], 0] = iel - connections[:, 1] = cnx.neighbors + connections[:, 1] = adj.neighbors grid = UnstructuredGrid( (mesh.nelements, diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 377f86c98b37d72de34cc8f8a1bbbc0f7ffa53b6..5c2b5c1cb6ec7be972ed1ca15e0eb2db1c4b852c 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -24,6 +24,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import six import numpy as np import modepy as mp import numpy.linalg as la @@ -39,13 +40,59 @@ __doc__ = """ :members: :undoc-members: -.. autoclass:: ElementConnectivity +.. autoclass:: NodalAdjacency +.. autoclass:: FacialAdjacencyGroup .. autofunction:: as_python +.. autofunction:: check_bc_coverage +.. autofunction:: is_boundary_tag_empty +Predefined Boundary tags +------------------------ + +.. autoclass:: BTAG_NONE +.. autoclass:: BTAG_ALL +.. autoclass:: BTAG_REALLY_ALL +.. autoclass:: BTAG_NO_BOUNDARY """ +# {{{ element tags + +class BTAG_NONE(object): + """A boundary tag representing an empty boundary or volume.""" + pass + + +class BTAG_ALL(object): + """A boundary tag representing the entire boundary or volume. + + In the case of the boundary, TAG_ALL does not include rank boundaries, + or, more generally, anything tagged with TAG_NO_BOUNDARY.""" + pass + + +class BTAG_REALLY_ALL(object): + """A boundary tag representing the entire boundary. + + Unlike :class:`TAG_ALL`, this includes rank boundaries, + or, more generally, everything tagged with :class:`TAG_NO_BOUNDARY`.""" + pass + + +class BTAG_NO_BOUNDARY(object): + """A boundary tag indicating that this edge should not fall under + :class:`TAG_ALL`. Among other things, this is used to keep rank boundaries + out of :class:`BTAG_ALL`. + """ + pass + + +SYSTEM_TAGS = set([BTAG_NONE, BTAG_ALL, BTAG_REALLY_ALL, BTAG_NO_BOUNDARY]) + +# }}} + + # {{{ element group class MeshElementGroup(Record): @@ -81,6 +128,11 @@ class MeshElementGroup(Record): *Not* the ambient dimension, see :attr:`Mesh.ambient_dim` for that. + .. automethod:: face_vertex_indices + .. automethod:: vertex_unit_coordinates + + .. attribute:: nfaces + .. automethod:: __eq__ .. automethod:: __ne__ """ @@ -139,6 +191,23 @@ class MeshElementGroup(Record): def nunit_nodes(self): return self.unit_nodes.shape[-1] + 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() + + def vertex_unit_coordinates(self): + """Return an array of shape ``(nfaces, dim)`` with the unit + coordinates of each vertex. + """ + raise NotImplementedError() + + @property + def nfaces(self): + return len(self.face_vertex_indices()) + def __eq__(self, other): return ( type(self) == type(other) @@ -182,7 +251,10 @@ class SimplexElementGroup(MeshElementGroup): raise TypeError("'dim' must be passed " "if 'unit_nodes' is not passed") - unit_nodes = mp.warp_and_blend_nodes(dim, order) + if dim <= 3: + unit_nodes = mp.warp_and_blend_nodes(dim, order) + else: + unit_nodes = mp.equidistant_nodes(dim, order) dims = unit_nodes.shape[0] @@ -191,12 +263,15 @@ class SimplexElementGroup(MeshElementGroup): "element. expected: %d, got: %d" % (dims+1, vertex_indices.shape[-1])) - MeshElementGroup.__init__(self, order, vertex_indices, nodes, + super(SimplexElementGroup, self).__init__(order, vertex_indices, nodes, element_nr_base, node_nr_base, unit_nodes, dim) def face_vertex_indices(self): if self.dim == 1: - return ((0, 1),) + return ( + (0,), + (1,), + ) elif self.dim == 2: return ( (0, 1), @@ -214,36 +289,21 @@ class SimplexElementGroup(MeshElementGroup): raise NotImplementedError("dim=%d" % self.dim) def vertex_unit_coordinates(self): - if self.dim == 1: - return np.array([ - [-1], [1] - ], dtype=np.float64) - elif self.dim == 2: - return np.array([ - [-1, -1], - [1, -1], - [-1, 1], - ], dtype=np.float64) - elif self.dim == 3: - return np.array([ - [-1, -1, -1], - [1, -1, -1], - [-1, 1, -1], - [-1, -1, 1], - ], dtype=np.float64) - else: - raise NotImplementedError("dim=%d" % self.dim) + from modepy.tools import unit_vertices + return unit_vertices(self.dim) # }}} -# {{{ mesh +# {{{ nodal adjacency + +class NodalAdjacency(Record): + """Describes nodal element adjacency information, i.e. information about + elements that touch in at least one point. -class ElementConnectivity(Record): - """ .. attribute:: neighbors_starts - ``element_id_t [nelments+1]`` + ``element_id_t [nelements+1]`` Use together with :attr:`neighbors`. ``neighbors_starts[iel]`` and ``neighbors_starts[iel+1]`` together indicate a ranges of element indices @@ -269,6 +329,77 @@ class ElementConnectivity(Record): def __ne__(self, other): return not self.__eq__(other) +# }}} + + +# {{{ facial adjacency + +class FacialAdjacencyGroup(Record): + """Describes facial element adjacency information for one + :class:`MeshElementGroup`, i.e. information about elements that share (part + of) a face. + + .. attribute:: igroup + + .. attribute:: ineighbor_group + + ID of neighboring group, or *None* for boundary faces. If identical + to :attr:`igroup`, then this contains the self-connectivity in this + group. + + .. attribute:: elements + + ``element_id_t [nfagrp_elements]``. ``elements[i]`` + Group-local element numbers. + + .. attribute:: element_faces + + ``face_id_t [nfagrp_elements]``. ``element_faces[i]`` + indicate what face index of the opposite element indicated in + ``neighbors[iel_grp][iface]`` touches face number *iface* of element + number *iel_grp* in this element group. + + .. attribute:: neighbors + + ``element_id_t [nfagrp_elements]``. ``neighbors[i]`` + gives the element number within :attr:`ineighbor_group` of the element + opposite ``elements[i]``. + + If this number is negative, then this indicates that this is a + boundary face, and the bits set in ``-neighbors[i]`` + should be interpreted according to :attr:`Mesh.boundary_tags`. + + .. attribute:: neighbor_faces + + ``face_id_t [nfagrp_elements]``. ``neighbor_faces[i]`` indicate what + face index of the opposite element indicated in ``neighbors[i]`` has + facial contact with face number ``element_faces[i]`` of element number + ``elements[i]`` in element group *igroup*. + + Zero if ``neighbors[i]`` is negative. + + .. automethod:: __eq__ + .. automethod:: __ne__ + """ + + def __eq__(self, other): + return ( + type(self) == type(other) + and self.igroup == other.igroup + and self.ineighbor_group == other.ineighbor_group + and np.array_equal(self.elements, other.elements) + and np.array_equal(self.element_faces, other.element_faces) + and np.array_equal(self.neighbors, other.neighbors) + and np.array_equal(self.neighbor_faces, other.neighbor_faces) + ) + + def __ne__(self, other): + return not self.__eq__(other) + +# }}} + + +# {{{ mesh class Mesh(Record): """ @@ -281,12 +412,36 @@ class Mesh(Record): A list of :class:`MeshElementGroup` instances. - .. attribute:: element_connectivity + .. attribute:: nodal_adjacency + + An instance of :class:`NodalAdjacency`. + + Referencing this attribute may raise + :exc:`meshmode.DataUnavailable`. + + .. attribute:: facial_adjacency_groups + + A list of mappings from neighbor group IDs to instances of + :class:`FacialAdjacencyGroup`. - An instance of :class:`ElementConnectivity`. + ``facial_adjacency_groups[igrp][ineighbor_group]`` gives + the set of facial adjacency relations between group *igrp* + and *ineighbor_group*. *ineighbor_group* and *igrp* may be + identical, or *ineighbor_group* may be *None*, in which case + a group containing boundary faces is returned. Referencing this attribute may raise - :exc:`meshmode.ConnectivityUnavailable`. + :exc:`meshmode.DataUnavailable`. + + .. attribute:: boundary_tags + + A tuple of boundary tag identifiers. :class:`BTAG_ALL` and + :class:`BTAG_REALLY_ALL` are guranateed to exist. + + .. attribute:: btag_to_index + + A mapping that maps boundary tag identifiers to their + corresponding index. .. attribute:: vertex_id_dtype @@ -296,8 +451,12 @@ class Mesh(Record): .. automethod:: __ne__ """ + face_id_dtype = np.int8 + def __init__(self, vertices, groups, skip_tests=False, - element_connectivity=False, + nodal_adjacency=False, + facial_adjacency_groups=False, + boundary_tags=None, vertex_id_dtype=np.int32, element_id_dtype=np.int32): """ @@ -305,7 +464,7 @@ class Mesh(Record): :arg skip_tests: Skip mesh tests, in case you want to load a broken mesh anyhow and then fix it inside of this data structure. - :arg element_connectivity: One of three options: + :arg nodal_adjacency: One of three options: *None*, in which case this information will be deduced from vertex adjacency. *False*, in which case this information will be marked unavailable (such as if there are @@ -313,8 +472,16 @@ class Mesh(Record): the full picture), and references to :attr:`element_neighbors_starts` and :attr:`element_neighbors` will result in exceptions. Lastly, a tuple - *(element_neighbors_starts, element_neighbors)*, representing the - correspondingly-named attributes. + :class:`NodalAdjacency` object. + :arg facial_adjacency_groups: One of three options: + *None*, in which case this information + will be deduced from vertex adjacency. *False*, in which case + this information will be marked unavailable (such as if there are + hanging nodes in the geometry, so that vertex adjacency does not convey + the full picture), and references to + :attr:`element_neighbors_starts` and :attr:`element_neighbors` + will result in exceptions. Lastly, a data structure as described in + :attr:`facial_adjacency_groups` may be passed. """ el_nr = 0 node_nr = 0 @@ -326,18 +493,48 @@ class Mesh(Record): el_nr += ng.nelements node_nr += ng.nnodes - if element_connectivity is not False and element_connectivity is not None: - nb_starts, nbs = element_connectivity - element_connectivity = ElementConnectivity( - neighbors_starts=nb_starts, - neighbors=nbs) + # {{{ boundary tags + + if boundary_tags is None: + boundary_tags = [] + else: + boundary_tags = boundary_tags[:] + + if BTAG_NONE in boundary_tags: + raise ValueError("BTAG_NONE is not allowed to be part of " + "boundary_tags") + if BTAG_ALL not in boundary_tags: + boundary_tags.append(BTAG_ALL) + if BTAG_REALLY_ALL not in boundary_tags: + boundary_tags.append(BTAG_REALLY_ALL) + + max_boundary_tag_count = int( + np.log(np.iinfo(element_id_dtype).max)/np.log(2)) + if len(boundary_tags) > max_boundary_tag_count: + raise ValueError("too few bits in element_id_dtype to represent all " + "boundary tags") + + btag_to_index = dict( + (btag, i) for i, btag in enumerate(boundary_tags)) + + # }}} + + if nodal_adjacency is not False and nodal_adjacency is not None: + if not isinstance(nodal_adjacency, NodalAdjacency): + nb_starts, nbs = nodal_adjacency + nodal_adjacency = NodalAdjacency( + neighbors_starts=nb_starts, + neighbors=nbs) - del nb_starts - del nbs + del nb_starts + del nbs Record.__init__( self, vertices=vertices, groups=new_groups, - _element_connectivity=element_connectivity, + _nodal_adjacency=nodal_adjacency, + _facial_adjacency_groups=facial_adjacency_groups, + boundary_tags=boundary_tags, + btag_to_index=btag_to_index, vertex_id_dtype=np.dtype(vertex_id_dtype), element_id_dtype=np.dtype(element_id_dtype), ) @@ -347,6 +544,32 @@ class Mesh(Record): for g in self.groups: assert g.vertex_indices.dtype == self.vertex_id_dtype + if nodal_adjacency: + assert nodal_adjacency.neighbor_starts.shape == (self.nelements+1,) + assert len(nodal_adjacency.neighbors.shape) == 1 + + assert nodal_adjacency.neighbor_starts.dtype == self.element_id_dtype + assert nodal_adjacency.neighbors.dtype == self.element_id_dtype + + if facial_adjacency_groups: + assert len(facial_adjacency_groups) == len(self.groups) + for fagrp_map in facial_adjacency_groups: + for fagrp in six.itervalues(fagrp_map): + grp = self.groups[fagrp.igroup] + + fvi = grp.face_vertex_indices() + assert fagrp.neighbors.dtype == self.element_id_dtype + assert fagrp.neighbors.shape == ( + grp.nelements, len(fvi)) + assert fagrp.neighbors.dtype == self.face_id_dtype + assert fagrp.neighbor_faces.shape == ( + grp.nelements, len(fvi)) + + is_bdry = fagrp.neighbors < 0 + assert ((1 << btag_to_index[BTAG_REALLY_ALL]) + & fagrp.neighbors[is_bdry]).all(), \ + "boundary faces without BTAG_REALLY_ALL found" + from meshmode.mesh.processing import \ test_volume_mesh_element_orientations @@ -364,19 +587,47 @@ class Mesh(Record): from pytools import single_valued return single_valued(grp.dim for grp in self.groups) + @property + def nvertices(self): + return self.vertices.shape[-1] + @property def nelements(self): return sum(grp.nelements for grp in self.groups) @property - def element_connectivity(self): - if self._element_connectivity is False: - from meshmode import ConnectivityUnavailable - raise ConnectivityUnavailable() - elif self._element_connectivity is None: - self._element_connectivity = _compute_connectivity_from_vertices(self) + def nodal_adjacency(self): + if self._nodal_adjacency is False: + from meshmode import DataUnavailable + raise DataUnavailable("nodal_adjacency") + elif self._nodal_adjacency is None: + self._nodal_adjacency = _compute_nodal_adjacency_from_vertices(self) + + return self._nodal_adjacency + + def nodal_adjacency_init_arg(self): + """Returns an 'nodal_adjacency' argument that can be + passed to a Mesh constructor. + """ - return self._element_connectivity + return self._nodal_adjacency + + @property + def facial_adjacency_groups(self): + if self._facial_adjacency_groups is False: + from meshmode import DataUnavailable + raise DataUnavailable("facial_adjacency_groups") + elif self._facial_adjacency_groups is None: + self._facial_adjacency_groups = \ + _compute_facial_adjacency_from_vertices(self) + + return self._facial_adjacency_groups + + def boundary_tag_bit(self, boundary_tag): + try: + return 1 << self.btag_to_index[boundary_tag] + except KeyError: + return 0 def __eq__(self, other): return ( @@ -385,8 +636,11 @@ class Mesh(Record): and self.groups == other.groups and self.vertex_id_dtype == other.vertex_id_dtype and self.element_id_dtype == other.element_id_dtype - and (self._element_connectivity - == other._element_connectivity)) + and (self._nodal_adjacency + == other._nodal_adjacency) + and (self._facial_adjacency_groups + == other._facial_adjacency_groups) + and self.boundary_tags == other.boundary_tags) def __ne__(self, other): return not self.__eq__(other) @@ -402,10 +656,13 @@ class Mesh(Record): # {{{ node-vertex consistency test def _test_node_vertex_consistency_simplex(mesh, mgrp): - from modepy.tools import UNIT_VERTICES + if mgrp.nelements == 0: + return True + resampling_mat = mp.resampling_matrix( - mp.simplex_onb(mgrp.dim, mgrp.order), - UNIT_VERTICES[mgrp.dim].T.copy(), mgrp.unit_nodes) + mp.simplex_best_available_basis(mgrp.dim, mgrp.order), + mgrp.vertex_unit_coordinates().T, + mgrp.unit_nodes) # dim, nelments, nnvertices map_vertices = np.einsum( @@ -448,9 +705,9 @@ def _test_node_vertex_consistency(mesh): # }}} -# {{{ vertex-based connectivity +# {{{ vertex-based nodal adjacency -def _compute_connectivity_from_vertices(mesh): +def _compute_nodal_adjacency_from_vertices(mesh): # FIXME Native code would make this faster _, nvertices = mesh.vertices.shape @@ -483,13 +740,122 @@ def _compute_connectivity_from_vertices(mesh): assert neighbors_starts[-1] == len(neighbors) - return ElementConnectivity( + return NodalAdjacency( neighbors_starts=neighbors_starts, neighbors=neighbors) # }}} +# {{{ vertex-based facial adjacency + +def _compute_facial_adjacency_from_vertices(mesh): + # FIXME Native code would make this faster + + # create face_map, which is a mapping of + # (vertices on a face) -> + # [(igrp, iel_grp, face_idx) for elements bordering that face] + face_map = {} + for igrp, grp in enumerate(mesh.groups): + for fid, face_vertex_indices in enumerate(grp.face_vertex_indices()): + all_fvi = grp.vertex_indices[:, face_vertex_indices] + + for iel_grp, fvi in enumerate(all_fvi): + face_map.setdefault( + frozenset(fvi), []).append((igrp, iel_grp, fid)) + + # maps tuples (igrp, ineighbor_group) to number of elements + group_count = {} + for face_tuples in six.itervalues(face_map): + if len(face_tuples) == 2: + (igrp, _, _), (inb_grp, _, _) = face_tuples + group_count[igrp, inb_grp] = group_count.get((igrp, inb_grp), 0) + 1 + group_count[inb_grp, igrp] = group_count.get((inb_grp, igrp), 0) + 1 + elif len(face_tuples) == 1: + (igrp, _, _), = face_tuples + group_count[igrp, None] = group_count.get((igrp, None), 0) + 1 + else: + raise RuntimeError("unexpected number of adjacent faces") + + # {{{ build facial_adjacency_groups data structure, still empty + + facial_adjacency_groups = [] + for igroup in range(len(mesh.groups)): + grp_map = {} + facial_adjacency_groups.append(grp_map) + + bdry_count = group_count.get((igroup, None)) + if bdry_count is not None: + elements = np.empty(bdry_count, dtype=mesh.element_id_dtype) + element_faces = np.empty(bdry_count, dtype=mesh.face_id_dtype) + neighbors = np.empty(bdry_count, dtype=mesh.element_id_dtype) + neighbor_faces = np.zeros(bdry_count, dtype=mesh.face_id_dtype) + + neighbors.fill(-( + mesh.boundary_tag_bit(BTAG_ALL) + | mesh.boundary_tag_bit(BTAG_REALLY_ALL))) + + grp_map[None] = FacialAdjacencyGroup( + igroup=igroup, + ineighbor_group=None, + elements=elements, + element_faces=element_faces, + neighbors=neighbors, + neighbor_faces=neighbor_faces) + + for ineighbor_group in range(len(mesh.groups)): + nb_count = group_count.get((igroup, ineighbor_group)) + if nb_count is not None: + elements = np.empty(nb_count, dtype=mesh.element_id_dtype) + element_faces = np.empty(nb_count, dtype=mesh.face_id_dtype) + neighbors = np.empty(nb_count, dtype=mesh.element_id_dtype) + neighbor_faces = np.empty(nb_count, dtype=mesh.face_id_dtype) + + grp_map[ineighbor_group] = FacialAdjacencyGroup( + igroup=igroup, + ineighbor_group=ineighbor_group, + elements=elements, + element_faces=element_faces, + neighbors=neighbors, + neighbor_faces=neighbor_faces) + + # }}} + + # maps tuples (igrp, ineighbor_group) to number of elements filled in group + fill_count = {} + for face_tuples in six.itervalues(face_map): + if len(face_tuples) == 2: + for (igroup, iel, iface), (inb_group, inb_el, inb_face) in [ + (face_tuples[0], face_tuples[1]), + (face_tuples[1], face_tuples[0]), + ]: + idx = fill_count.get((igrp, inb_grp), 0) + fill_count[igrp, inb_grp] = idx + 1 + + fagrp = facial_adjacency_groups[igroup][inb_grp] + fagrp.elements[idx] = iel + fagrp.element_faces[idx] = iface + fagrp.neighbors[idx] = inb_el + fagrp.neighbor_faces[idx] = inb_face + + elif len(face_tuples) == 1: + (igroup, iel, iface), = face_tuples + + idx = fill_count.get((igrp, None), 0) + fill_count[igrp, None] = idx + 1 + + fagrp = facial_adjacency_groups[igroup][None] + fagrp.elements[idx] = iel + fagrp.element_faces[idx] = iface + + else: + raise RuntimeError("unexpected number of adjacent faces") + + return facial_adjacency_groups + +# }}} + + # {{{ as_python def _numpy_array_as_python(array): @@ -505,11 +871,15 @@ def as_python(mesh, function_name="make_mesh"): from pytools.py_codegen import PythonCodeGenerator, Indentation cg = PythonCodeGenerator() - cg("# generated by meshmode.mesh.as_python") - cg("") - cg("import numpy as np") - cg("from meshmode.mesh import Mesh, MeshElementGroup") - cg("") + cg(""" + # generated by meshmode.mesh.as_python + + import numpy as np + from meshmode.mesh import ( + Mesh, MeshElementGroup, FacialAdjacencyGroup, + BTAG_ALL, BTAG_REALLY_ALL) + + """) cg("def %s():" % function_name) with Indentation(cg): @@ -530,25 +900,136 @@ def as_python(mesh, function_name="make_mesh"): cg(" unit_nodes=%s))" % _numpy_array_as_python(group.unit_nodes)) + # {{{ facial adjacency groups + + def fagrp_params_str(fagrp): + params = { + "igroup": fagrp.igroup, + "ineighbor_group": repr(fagrp.ineighbor_group), + "elements": _numpy_array_as_python(fagrp.elements), + "element_faces": _numpy_array_as_python(fagrp.element_faces), + "neighbors": _numpy_array_as_python(fagrp.neighbors), + "neighbor_faces": _numpy_array_as_python(fagrp.neighbor_faces), + } + return ",\n ".join("%s=%s" % (k, v) for k, v in six.iteritems(params)) + + if mesh._facial_adjacency_groups: + cg("facial_adjacency_groups = []") + + for igrp, fagrp_map in enumerate(mesh.facial_adjacency_groups): + cg("facial_adjacency_groups.append({%s})" % ",\n ".join( + "%r: FacialAdjacencyGroup(%s)" % ( + inb_grp, fagrp_params_str(fagrp)) + for inb_grp, fagrp in six.iteritems(fagrp_map))) + + else: + cg("facial_adjacency_groups = %r" % mesh._facial_adjacency_groups) + + # }}} + + # {{{ boundary tags + + def strify_boundary_tag(btag): + if isinstance(btag, type): + return btag.__name__ + else: + return repr(btag) + + btags_str = ", ".join( + strify_boundary_tag(btag) for btag in mesh.boundary_tags) + + # }}} + cg("return Mesh(vertices, groups, skip_tests=True,") cg(" vertex_id_dtype=np.%s," % mesh.vertex_id_dtype.name) cg(" element_id_dtype=np.%s," % mesh.element_id_dtype.name) - if isinstance(mesh._element_connectivity, ElementConnectivity): + if isinstance(mesh._nodal_adjacency, NodalAdjacency): el_con_str = "(%s, %s)" % ( _numpy_array_as_python( - mesh._element_connectivity.neighbors_starts), + mesh._nodal_adjacency.neighbors_starts), _numpy_array_as_python( - mesh._element_connectivity.neighbors), + mesh._nodal_adjacency.neighbors), ) else: - el_con_str = repr(mesh._element_connectivity) + el_con_str = repr(mesh._nodal_adjacency) - cg(" element_connectivity=%s)" % el_con_str) + cg(" nodal_adjacency=%s," % el_con_str) + cg(" facial_adjacency_groups=facial_adjacency_groups,") + cg(" boundary_tags=[%s])" % btags_str) + + # FIXME: Handle facial adjacency, boundary tags return cg.get() # }}} +# {{{ check_bc_coverage + +def check_bc_coverage(mesh, boundary_tags, incomplete_ok=False): + """Verify boundary condition coverage. + + Given a list of boundary tags as *boundary_tags*, this function verifies + that + + 1. the union of all these boundaries gives the complete boundary, + 1. all these boundaries are disjoint. + + :arg incomplete_ok: Do not report an error if some faces are not covered + by the boundary conditions. + """ + + for igrp, fagrp_map in enumerate(mesh.facial_adjacency_groups): + bdry_grp = fagrp_map.get(None) + if bdry_grp is None: + continue + + nb_elements = bdry_grp.neighbors + assert (nb_elements < 0).all() + + nb_el_bits = -nb_elements + + seen = np.zeros_like(nb_el_bits, dtype=np.bool) + + for btag in boundary_tags: + tag_bit = mesh.boundary_tag_bit(btag) + tag_set = (nb_el_bits & tag_bit) != 0 + + if (seen & tag_set).any(): + raise RuntimeError("faces with multiple boundary conditions found") + + seen = seen | tag_set + + if not incomplete_ok and not seen.all(): + raise RuntimeError("found faces without boundary conditions") + +# }}} + + +# {{{ is_boundary_tag_empty + +def is_boundary_tag_empty(mesh, boundary_tag): + """Return *True* if the corresponding boundary tag does not occur as part of + *mesh*. + """ + + btag_bit = mesh.boundary_tag_bit(boundary_tag) + if not btag_bit: + return True + + for igrp in range(len(mesh.groups)): + bdry_fagrp = mesh.facial_adjacency_groups[igrp].get(None, None) + if bdry_fagrp is None: + continue + + neg = bdry_fagrp.neighbors < 0 + if (-bdry_fagrp.neighbors[neg] & btag_bit).any(): + return False + + return True + +# }}} + + # vim: foldmethod=marker diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 9a72c5377c28ac876ff13d1391393860c7dc7986..65761586e01c8b0aa60f82d6f8b4bc00941eeed7 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.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) 2013 Andreas Kloeckner" @@ -24,6 +22,7 @@ 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 import numpy.linalg as la @@ -59,6 +58,7 @@ Volumes .. autofunction:: generate_box_mesh .. autofunction:: generate_regular_rect_mesh +.. autofunction:: generate_warped_rect_mesh """ @@ -212,8 +212,10 @@ def make_curve_mesh(curve_f, element_boundaries, order): nodes=nodes, unit_nodes=unodes) - return Mesh(vertices=vertices, groups=[egroup], - element_connectivity=None) + return Mesh( + vertices=vertices, groups=[egroup], + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} @@ -232,7 +234,11 @@ def make_group_from_vertices(vertices, vertex_indices, order): dim = nspan_vectors # dim, nunit_nodes - unit_nodes = mp.warp_and_blend_nodes(dim, order) + if dim <= 3: + unit_nodes = mp.warp_and_blend_nodes(dim, order) + else: + unit_nodes = mp.equidistant_nodes(dim, order) + unit_nodes_01 = 0.5 + 0.5*unit_nodes nodes = np.einsum( @@ -285,8 +291,10 @@ def generate_icosahedron(r, order): grp = make_group_from_vertices(vertices, vertex_indices, order) from meshmode.mesh import Mesh - return Mesh(vertices, [grp], - element_connectivity=None) + return Mesh( + vertices, [grp], + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} @@ -302,8 +310,10 @@ def generate_icosphere(r, order): nodes=grp.nodes * r / np.sqrt(np.sum(grp.nodes**2, axis=0))) from meshmode.mesh import Mesh - return Mesh(mesh.vertices, [grp], - element_connectivity=None) + return Mesh( + mesh.vertices, [grp], + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} @@ -365,7 +375,11 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner, nodes[2] = b*np.sin(minor_theta) from meshmode.mesh import Mesh - return (Mesh(vertices, [grp.copy(nodes=nodes)], element_connectivity=None), + return ( + Mesh( + vertices, [grp.copy(nodes=nodes)], + nodal_adjacency=None, + facial_adjacency_groups=None), [idx(i, 0) for i in range(n_outer)], [idx(0, j) for j in range(n_inner)]) @@ -470,7 +484,8 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64): from meshmode.mesh import Mesh return Mesh(vertices, [grp], - element_connectivity=None) + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} @@ -495,4 +510,35 @@ 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): + """Generate a mesh of a warped line/square/cube. Mainly useful for testing + functionality with curvilinear meshes. + """ + + assert dim in [1, 2, 3] + mesh = generate_regular_rect_mesh( + a=(-0.5,)*dim, b=(0.5,)*dim, + n=(n,)*dim, order=order) + + def m(x): + result = np.empty_like(x) + result[0] = ( + 1.5*x[0] + np.cos(x[0]) + + 0.1*np.sin(10*x[1])) + result[1] = ( + 0.05*np.cos(10*x[0]) + + 1.3*x[1] + np.sin(x[1])) + if len(x) == 3: + result[2] = x[2] + np.sin(x[0]) + return result + + from meshmode.mesh.processing import map_mesh + return map_mesh(mesh, m) + +# }}} + + # vim: fdm=marker diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 292ee55077a178c70694bf8fd719b2fe42eef143..e27f8de352016e6f70c86be905fd4e247213d72b 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -1,8 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -import six -from six.moves import range -from six.moves import zip +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" @@ -26,6 +22,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import six +from six.moves import range, zip import numpy as np from meshpy.gmsh_reader import ( # noqa @@ -39,6 +37,8 @@ __doc__ = """ .. autofunction:: read_gmsh .. autofunction:: generate_gmsh +.. autofunction:: from_meshpy +.. autofunction:: from_vertices_and_simplices """ @@ -171,11 +171,16 @@ class GmshMeshReceiver(GmshMeshReceiverBase): groups.append(group) - return Mesh(vertices, groups, element_connectivity=None) + return Mesh( + vertices, groups, + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} +# {{{ gmsh + def read_gmsh(filename, force_ambient_dim=None): """Read a gmsh mesh file from *filename* and return a :class:`meshmode.mesh.Mesh`. @@ -213,7 +218,7 @@ def generate_gmsh(source, dimensions, order=None, other_options=[], mesh = recv.get_mesh() if force_ambient_dim is None: - AXIS_NAMES = "xyz" + AXIS_NAMES = "xyz" # noqa dim = mesh.vertices.shape[0] for idim in range(dim): @@ -226,3 +231,66 @@ def generate_gmsh(source, dimensions, order=None, other_options=[], break return mesh + +# }}} + + +# {{{ meshpy + +def from_meshpy(mesh_info, order=1): + """Imports a mesh from a :mod:`meshpy` *mesh_info* data structure, + which may be generated by either :mod:`meshpy.triangle` or + :mod:`meshpy_tet`. + """ + from meshmode.mesh import Mesh + from meshmode.mesh.generation import make_group_from_vertices + + vertices = np.array(mesh_info.points).T + elements = np.array(mesh_info.elements, np.int32) + + grp = make_group_from_vertices(vertices, elements, order) + + # FIXME: Should transfer boundary/volume markers + + return Mesh( + vertices=vertices, groups=[grp], + nodal_adjacency=None, + facial_adjacency_groups=None) + +# }}} + + +# {{{ from_vertices_and_simplices + +def from_vertices_and_simplices(vertices, simplices, order=1, fix_orientation=False): + """Imports a mesh from a numpy array of vertices and an array + of simplices. + + :arg vertices: + An array of vertex coordinates with shape + *(ambient_dim, nvertices)* + :arg simplices: + An array *(nelements, nvertices)* of (mesh-wide) + vertex indices. + """ + from meshmode.mesh import Mesh + from meshmode.mesh.generation import make_group_from_vertices + + grp = make_group_from_vertices(vertices, simplices, order) + + if fix_orientation: + from meshmode.mesh.processing import ( + find_volume_mesh_element_group_orientation, + flip_simplex_element_group) + orient = find_volume_mesh_element_group_orientation(vertices, grp) + grp = flip_simplex_element_group(vertices, grp, orient < 0) + + return Mesh( + vertices=vertices, groups=[grp], + nodal_adjacency=None, + facial_adjacency_groups=None) + +# }}} + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 94cd8b0951f5eb937c556bc239a2ebd7110bfb04..06631f9dedd6fdb1bcd3112a8c8157d4c72cb7f7 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -1,7 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range -from functools import reduce +from __future__ import division, absolute_import, print_function __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" @@ -25,6 +22,9 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range +from functools import reduce + import numpy as np import numpy.linalg as la import modepy as mp @@ -34,14 +34,15 @@ __doc__ = """ .. autofunction:: find_volume_mesh_element_orientations .. autofunction:: perform_flips .. autofunction:: find_bounding_box -.. autofunction:: merge_dijsoint_meshes +.. autofunction:: merge_disjoint_meshes +.. autofunction:: map_mesh .. autofunction:: affine_map """ # {{{ orientations -def find_volume_mesh_element_group_orientation(mesh, grp): +def find_volume_mesh_element_group_orientation(vertices, grp): """Return a positive floating point number for each positively oriented element, and a negative floating point number for each negatively oriented element. @@ -56,11 +57,11 @@ def find_volume_mesh_element_group_orientation(mesh, grp): "exclusively SimplexElementGroup-based meshes") # (ambient_dim, nelements, nvertices) - vertices = mesh.vertices[:, grp.vertex_indices] + my_vertices = vertices[:, grp.vertex_indices] # (ambient_dim, nelements, nspan_vectors) spanning_vectors = ( - vertices[:, :, 1:] - vertices[:, :, 0][:, :, np.newaxis]) + my_vertices[:, :, 1:] - my_vertices[:, :, 0][:, :, np.newaxis]) ambient_dim = spanning_vectors.shape[0] nspan_vectors = spanning_vectors.shape[-1] @@ -81,6 +82,10 @@ def find_volume_mesh_element_group_orientation(mesh, grp): from operator import xor outer_prod = -reduce(xor, mvs) + if grp.dim == 1: + # FIXME: This is a little weird. + outer_prod = -outer_prod + return (outer_prod.I | outer_prod).as_scalar() @@ -102,7 +107,8 @@ def find_volume_mesh_element_orientations(mesh, tolerate_unimplemented_checks=Fa if tolerate_unimplemented_checks: try: signed_area_elements = \ - find_volume_mesh_element_group_orientation(mesh, grp) + find_volume_mesh_element_group_orientation( + mesh.vertices, grp) except NotImplementedError: result_grp_view[:] = float("nan") else: @@ -110,7 +116,8 @@ def find_volume_mesh_element_orientations(mesh, tolerate_unimplemented_checks=Fa result_grp_view[:] = signed_area_elements else: signed_area_elements = \ - find_volume_mesh_element_group_orientation(mesh, grp) + find_volume_mesh_element_group_orientation( + mesh.vertices, grp) assert not np.isnan(signed_area_elements).any() result_grp_view[:] = signed_area_elements @@ -158,7 +165,7 @@ def flip_simplex_element_group(vertices, grp, grp_flip_flags): flipped_unit_nodes = barycentric_to_unit(flipped_bary_unit_nodes) flip_matrix = mp.resampling_matrix( - mp.simplex_onb(grp.dim, grp.order), + mp.simplex_best_available_basis(grp.dim, grp.order), flipped_unit_nodes, grp.unit_nodes) flip_matrix[np.abs(flip_matrix) < 1e-15] = 0 @@ -220,7 +227,7 @@ def find_bounding_box(mesh): # {{{ merging -def merge_dijsoint_meshes(meshes, skip_tests=False): +def merge_disjoint_meshes(meshes, skip_tests=False): if not meshes: raise ValueError("must pass at least one mesh") @@ -271,27 +278,48 @@ def merge_dijsoint_meshes(meshes, skip_tests=False): # }}} -# {{{ affine map +# {{{ map -def affine_map(mesh, A=None, b=None): - """Apply the affine map *f(x)=Ax+b* to the geometry of *mesh*.""" +def map_mesh(mesh, f): # noqa + """Apply the map *f* to the mesh. *f* needs to accept and return arrays of + shape ``(ambient_dim, npoints)``.""" - vertices = np.einsum("ij,jv->iv", A, mesh.vertices) + b[:, np.newaxis] + vertices = f(mesh.vertices) # {{{ assemble new groups list new_groups = [] for group in mesh.groups: - nodes = ( - np.einsum("ij,jen->ien", A, group.nodes) - + b[:, np.newaxis, np.newaxis]) - new_groups.append(group.copy(nodes=nodes)) + mapped_nodes = f(group.nodes.reshape(mesh.ambient_dim, -1)) + new_groups.append(group.copy( + nodes=mapped_nodes.reshape(*group.nodes.shape))) # }}} from meshmode.mesh import Mesh - return Mesh(vertices, new_groups, skip_tests=True) + return Mesh(vertices, new_groups, skip_tests=True, + nodal_adjacency=mesh.nodal_adjacency_init_arg(), + facial_adjacency_groups=mesh._facial_adjacency_groups) + +# }}} + + +# {{{ affine map + +def affine_map(mesh, A=None, b=None): # noqa + """Apply the affine map *f(x)=Ax+b* to the geometry of *mesh*.""" + + if A is None: + A = np.eye(mesh.ambient_dim) # noqa + + if b is None: + b = np.zeros(A.shape[0]) + + def f(x): + return np.einsum("ij,jv->iv", A, x) + b[:, np.newaxis] + + return map_mesh(mesh, f) # }}} diff --git a/meshmode/mesh/refinement.py b/meshmode/mesh/refinement.py index 23fe893dc9b3383ced0088981a9e409b78965815..0a6518f0ec040c3b1a0adc8f18900a65ff197caf 100644 --- a/meshmode/mesh/refinement.py +++ b/meshmode/mesh/refinement.py @@ -22,6 +22,8 @@ THE SOFTWARE. import numpy as np +import itertools +from six.moves import range class TreeRayNode(object): @@ -39,7 +41,7 @@ class TreeRayNode(object): List containing elements indices of elements adjacent to this ray segment. """ - def __init__(self, left_vertex, right_vertex, adjacent_elements = []): + def __init__(self, left_vertex, right_vertex, adjacent_elements=[]): import copy self.left = None self.right = None @@ -50,16 +52,15 @@ class TreeRayNode(object): self.adjacent_elements = copy.deepcopy(adjacent_elements) self.adjacent_add_diff = [] + class Refiner(object): def __init__(self, mesh): - #print 'herlkjjlkjasdf' - from llist import dllist, dllistnode - from meshmode.mesh.tesselate import tesselatetet, tesselatetri + from meshmode.mesh.tesselate import tesselatetet, tesselatetri self.lazy = False self.seen_tuple = {} tri_node_tuples, tri_result = tesselatetri() tet_node_tuples, tet_result = tesselatetet() - print tri_node_tuples, tri_result + print(tri_node_tuples, tri_result) #quadrilateral_node_tuples = [ #print tri_result, tet_result self.simplex_node_tuples = [None, None, tri_node_tuples, tet_node_tuples] @@ -67,38 +68,36 @@ class Refiner(object): #print tri_node_tuples, tri_result #self.simplex_node_tuples, self.simplex_result = tesselatetet() self.last_mesh = mesh - + # {{{ initialization - # mapping: (vertex_1, vertex_2) -> TreeRayNode + # mapping: (vertex_1, vertex_2) -> TreeRayNode # where vertex_i represents a vertex number # # Assumption: vertex_1 < vertex_2 self.pair_map = {} nvertices = len(mesh.vertices[0]) - + #array containing element whose edge lies on corresponding vertex self.hanging_vertex_element = [] - import six - for i in six.moves.range(nvertices): + for i in range(nvertices): self.hanging_vertex_element.append([]) - import six for grp in mesh.groups: iel_base = grp.element_nr_base - for iel_grp in six.moves.range(grp.nelements): + for iel_grp in range(grp.nelements): vert_indices = grp.vertex_indices[iel_grp] - for i in six.moves.range(len(vert_indices)): - for j in six.moves.range(i+1, len(vert_indices)): - - #minimum and maximum of the two indices for storing + for i in range(len(vert_indices)): + for j in range(i+1, len(vert_indices)): + + #minimum and maximum of the two indices for storing #data in vertex_pair min_index = min(vert_indices[i], vert_indices[j]) max_index = max(vert_indices[i], vert_indices[j]) vertex_pair = (min_index, max_index) - #print vertex_pair + #print(vertex_pair) if vertex_pair not in self.pair_map: self.pair_map[vertex_pair] = TreeRayNode(min_index, max_index) self.pair_map[vertex_pair].adjacent_elements.append(iel_base+iel_grp) @@ -107,39 +106,39 @@ class Refiner(object): adjacent_elements.append(iel_base+iel_grp)) # }}} - #print vert_indices + #print(vert_indices) #generate reference tuples self.index_to_node_tuple = [] self.index_to_midpoint_tuple = [] - for d in six.moves.range(len(vert_indices)): + for d in range(len(vert_indices)): dim = d + 1 cur_index_to_node_tuple = [()] * dim - for i in six.moves.range(0, dim-1): + for i in range(0, dim-1): cur_index_to_node_tuple[0] = cur_index_to_node_tuple[0] + (0,) - for i in six.moves.range(1, dim): - for j in six.moves.range(1, dim): + for i in range(1, dim): + for j in range(1, dim): if i == j: cur_index_to_node_tuple[i] = cur_index_to_node_tuple[i] + (2,) else: cur_index_to_node_tuple[i] = cur_index_to_node_tuple[i] + (0,) cur_index_to_midpoint_tuple = [()] * (int((dim * (dim - 1)) / 2)) curind = 0 - for ind1 in six.moves.range(0, len(cur_index_to_node_tuple)): - for ind2 in six.moves.range(ind1+1, len(cur_index_to_node_tuple)): + for ind1 in range(0, len(cur_index_to_node_tuple)): + for ind2 in range(ind1+1, len(cur_index_to_node_tuple)): i = cur_index_to_node_tuple[ind1] j = cur_index_to_node_tuple[ind2] - #print i, j - for k in six.moves.range(0, dim-1): + #print(i, j) + for k in range(0, dim-1): cur = int((i[k] + j[k]) / 2) cur_index_to_midpoint_tuple[curind] = cur_index_to_midpoint_tuple[curind] + (cur,) curind += 1 self.index_to_node_tuple.append(cur_index_to_node_tuple) self.index_to_midpoint_tuple.append(cur_index_to_midpoint_tuple) ''' - self.ray_vertices = np.empty([len(mesh.groups[0].vertex_indices) * - len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1) / 2, 2], + self.ray_vertices = np.empty([len(mesh.groups[0].vertex_indices) * + len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1) / 2, 2], dtype=np.int32) - self.ray_elements = np.zeros([len(mesh.groups[0].vertex_indices) * + self.ray_elements = np.zeros([len(mesh.groups[0].vertex_indices) * len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1) / 2, 1, 2], dtype=np.int32) self.vertex_to_ray = [] @@ -158,7 +157,7 @@ class Refiner(object): self.vertex_to_ray[grp.vertex_indices[i][k]].append(temp2) count += 1 ind = 0 - #print self.ray_vertices + #print(self.ray_vertices) for i in self.ray_vertices: which = 0 for grp in mesh.groups: @@ -188,16 +187,16 @@ class Refiner(object): np.bool) def get_current_mesh(self): - + from meshmode.mesh import Mesh #return Mesh(vertices, [grp], element_connectivity=self.generate_connectivity(len(self.last_mesh.groups[group].vertex_indices) \ # + count*3)) - groups = [] + groups = [] grpn = 0 for grp in self.last_mesh.groups: groups.append(np.empty([len(grp.vertex_indices), len(self.last_mesh.groups[grpn].vertex_indices[0])], dtype=np.int32)) - for iel_grp in xrange(grp.nelements): + for iel_grp in range(grp.nelements): for i in range(0, len(grp.vertex_indices[iel_grp])): groups[grpn][iel_grp][i] = grp.vertex_indices[iel_grp][i] grpn += 1 @@ -209,9 +208,9 @@ class Refiner(object): self.last_mesh = Mesh(self.last_mesh.vertices, grp,\ element_connectivity=self.generate_connectivity(len(self.last_mesh.groups[0].vertex_indices),\ len(self.last_mesh.vertices[0]))) - + return self.last_mesh - + def get_leaves(self, cur_node): queue = [cur_node] @@ -236,7 +235,7 @@ class Refiner(object): queue.append(vertex.left) queue.append(vertex.right) return res - + def apply_diff(self, cur_node, new_hanging_vertex_elements): for el in cur_node.adjacent_add_diff: if el not in cur_node.adjacent_elements: @@ -311,30 +310,27 @@ class Refiner(object): #refine_flag tells you which elements to split as a numpy array of bools def refine(self, refine_flags): - import six import numpy as np - from sets import Set #vertices and groups for next generation nvertices = len(self.last_mesh.vertices[0]) groups = [] - midpoint_already = Set([]) + midpoint_already = set() grpn = 0 totalnelements = 0 - for grp in self.last_mesh.groups: iel_base = grp.element_nr_base nelements = 0 - for iel_grp in six.moves.range(grp.nelements): + for iel_grp in range(grp.nelements): nelements += 1 vertex_indices = grp.vertex_indices[iel_grp] if refine_flags[iel_base+iel_grp]: cur_dim = len(grp.vertex_indices[iel_grp])-1 nelements += len(self.simplex_result[cur_dim]) - 1 - for i in six.moves.range(len(vertex_indices)): - for j in six.moves.range(i+1, len(vertex_indices)): + for i in range(len(vertex_indices)): + for j in range(i+1, len(vertex_indices)): i_index = vertex_indices[i] j_index = vertex_indices[j] index_tuple = (i_index, j_index) if i_index < j_index else (j_index, i_index) @@ -349,11 +345,10 @@ class Refiner(object): vertices = np.empty([len(self.last_mesh.vertices), nvertices]) new_hanging_vertex_element = [ - [] for i in six.moves.range(nvertices)] + [] for i in range(nvertices)] # def remove_element_from_connectivity(vertices, new_hanging_vertex_elements, to_remove): -# #print vertices -# import six +# #print(vertices) # import itertools # if len(vertices) == 2: # min_vertex = min(vertices[0], vertices[1]) @@ -366,8 +361,8 @@ class Refiner(object): # element_rays = [] # midpoints = [] # split_possible = True -# for i in six.moves.range(len(vertices)): -# for j in six.moves.range(i+1, len(vertices)): +# for i in range(len(vertices)): +# for j in range(i+1, len(vertices)): # min_vertex = min(vertices[i], vertices[j]) # max_vertex = max(vertices[i], vertices[j]) # element_rays.append(self.pair_map[(min_vertex, max_vertex)]) @@ -385,13 +380,13 @@ class Refiner(object): # node_tuple_to_coord[node_tuple] = vertices[node_index] # for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]): # node_tuple_to_coord[midpoint_tuple] = midpoints[midpoint_index] -# for i in six.moves.range(len(self.simplex_result[cur_dim])): +# for i in range(len(self.simplex_result[cur_dim])): # next_vertices = [] -# for j in six.moves.range(len(self.simplex_result[cur_dim][i])): +# for j in range(len(self.simplex_result[cur_dim][i])): # next_vertices.append(node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]]) # all_rays_present = True -# for v1 in six.moves.range(len(next_vertices)): -# for v2 in six.moves.range(v1+1, len(next_vertices)): +# for v1 in range(len(next_vertices)): +# for v2 in range(v1+1, len(next_vertices)): # vertex_tuple = (min(next_vertices[v1], next_vertices[v2]), max(next_vertices[v1], next_vertices[v2])) # if vertex_tuple not in self.pair_map: # all_rays_present = False @@ -400,13 +395,11 @@ class Refiner(object): # else: # split_possible = False # if not split_possible: -# next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1)) +# next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1)) # for next_vertices in next_vertices_list: # remove_element_from_connectivity(next_vertices, new_hanging_vertex_elements, to_remove) def add_element_to_connectivity(vertices, new_hanging_vertex_elements, to_add): - import six - import itertools if len(vertices) == 2: min_vertex = min(vertices[0], vertices[1]) max_vertex = max(vertices[0], vertices[1]) @@ -418,8 +411,8 @@ class Refiner(object): element_rays = [] midpoints = [] split_possible = True - for i in six.moves.range(len(vertices)): - for j in six.moves.range(i+1, len(vertices)): + for i in range(len(vertices)): + for j in range(i+1, len(vertices)): min_vertex = min(vertices[i], vertices[j]) max_vertex = max(vertices[i], vertices[j]) element_rays.append(self.pair_map[(min_vertex, max_vertex)]) @@ -435,13 +428,13 @@ class Refiner(object): node_tuple_to_coord[node_tuple] = vertices[node_index] for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]): node_tuple_to_coord[midpoint_tuple] = midpoints[midpoint_index] - for i in six.moves.range(len(self.simplex_result[cur_dim])): + for i in range(len(self.simplex_result[cur_dim])): next_vertices = [] - for j in six.moves.range(len(self.simplex_result[cur_dim][i])): + for j in range(len(self.simplex_result[cur_dim][i])): next_vertices.append(node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]]) all_rays_present = True - for v1 in six.moves.range(len(next_vertices)): - for v2 in six.moves.range(v1+1, len(next_vertices)): + for v1 in range(len(next_vertices)): + for v2 in range(v1+1, len(next_vertices)): vertex_tuple = (min(next_vertices[v1], next_vertices[v2]), max(next_vertices[v1], next_vertices[v2])) if vertex_tuple not in self.pair_map: all_rays_present = False @@ -450,10 +443,9 @@ class Refiner(object): else: split_possible = False if not split_possible: - next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1)) + next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1)) for next_vertices in next_vertices_list: add_element_to_connectivity(next_vertices, new_hanging_vertex_elements, to_add) -# import six # for node in element_rays: # self.add_element_to_connectivity(node, new_hanging_vertex_elements, to_add) # leaves = self.get_subtree(node) @@ -466,8 +458,8 @@ class Refiner(object): # new_hanging_vertex_elements[leaf.right_vertex].append(to_add) # next_element_rays = [] -# for i in six.moves.range(len(element_rays)): -# for j in six.moves.range(i+1, len(element_rays)): +# for i in range(len(element_rays)): +# for j in range(i+1, len(element_rays)): # if element_rays[i].midpoint is not None and element_rays[j].midpoint is not None: # min_midpoint = min(element_rays[i].midpoint, element_rays[j].midpoint) # max_midpoint = max(element_rays[i].midpoint, element_rays[j].midpoint) @@ -503,10 +495,10 @@ class Refiner(object): def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_grp): for grp in groups: iel_base = 0 - for iel_grp in six.moves.range(nelements_in_grp): + for iel_grp in range(nelements_in_grp): vertex_indices = grp[iel_grp] - for i in six.moves.range(len(vertex_indices)): - for j in six.moves.range(i+1, len(vertex_indices)): + for i in range(len(vertex_indices)): + for j in range(i+1, len(vertex_indices)): min_index = min(vertex_indices[i], vertex_indices[j]) max_index = max(vertex_indices[i], vertex_indices[j]) cur_node = self.pair_map[(min_index, max_index)] @@ -519,16 +511,16 @@ class Refiner(object): assert((iel_base+iel_grp) in new_hanging_vertex_elements[cur_node.left_vertex]) assert((iel_base+iel_grp) in new_hanging_vertex_elements[cur_node.right_vertex]) - for i in six.moves.range(len(self.last_mesh.vertices)): - for j in six.moves.range(len(self.last_mesh.vertices[i])): + for i in range(len(self.last_mesh.vertices)): + for j in range(len(self.last_mesh.vertices[i])): vertices[i,j] = self.last_mesh.vertices[i,j] import copy if i == 0: new_hanging_vertex_element[j] = copy.deepcopy(self.hanging_vertex_element[j]) grpn = 0 for grp in self.last_mesh.groups: - for iel_grp in six.moves.range(grp.nelements): - for i in six.moves.range(len(grp.vertex_indices[iel_grp])): + for iel_grp in range(grp.nelements): + for i in range(len(grp.vertex_indices[iel_grp])): groups[grpn][iel_grp][i] = grp.vertex_indices[iel_grp][i] grpn += 1 @@ -537,15 +529,15 @@ class Refiner(object): nelements_in_grp = grp.nelements for grp in self.last_mesh.groups: iel_base = grp.element_nr_base - for iel_grp in six.moves.range(grp.nelements): + for iel_grp in range(grp.nelements): if refine_flags[iel_base+iel_grp]: midpoint_vertices = [] midpoint_tuples = [] vertex_indices = grp.vertex_indices[iel_grp] #if simplex if len(grp.vertex_indices[iel_grp]) == len(self.last_mesh.vertices)+1: - for i in six.moves.range(len(vertex_indices)): - for j in six.moves.range(i+1, len(vertex_indices)): + for i in range(len(vertex_indices)): + for j in range(i+1, len(vertex_indices)): min_index = min(vertex_indices[i], vertex_indices[j]) max_index = max(vertex_indices[i], vertex_indices[j]) cur_node = self.pair_map[(min_index, max_index)] @@ -562,7 +554,7 @@ class Refiner(object): vertex_pair2 = (max_index, vertices_index) self.pair_map[vertex_pair1] = cur_node.left self.pair_map[vertex_pair2] = cur_node.right - for k in six.moves.range(len(self.last_mesh.vertices)): + for k in range(len(self.last_mesh.vertices)): vertices[k, vertices_index] = \ (self.last_mesh.vertices[k, vertex_indices[i]] + self.last_mesh.vertices[k, vertex_indices[j]]) / 2.0 @@ -575,8 +567,8 @@ class Refiner(object): #generate new rays cur_dim = len(grp.vertex_indices[0])-1 - for i in six.moves.range(len(midpoint_vertices)): - for j in six.moves.range(i+1, len(midpoint_vertices)): + for i in range(len(midpoint_vertices)): + for j in range(i+1, len(midpoint_vertices)): min_index = min(midpoint_vertices[i], midpoint_vertices[j]) max_index = max(midpoint_vertices[i], midpoint_vertices[j]) vertex_pair = (min_index, max_index) @@ -588,8 +580,8 @@ class Refiner(object): node_tuple_to_coord[node_tuple] = grp.vertex_indices[iel_grp][node_index] for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]): node_tuple_to_coord[midpoint_tuple] = midpoint_vertices[midpoint_index] - for i in six.moves.range(len(self.simplex_result[cur_dim])): - for j in six.moves.range(len(self.simplex_result[cur_dim][i])): + for i in range(len(self.simplex_result[cur_dim])): + for j in range(len(self.simplex_result[cur_dim][i])): if i == 0: groups[grpn][iel_grp][j] = \ node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]] @@ -605,14 +597,14 @@ class Refiner(object): # node_tuple_to_coord[node_tuple] = grp.vertex_indices[iel_grp][node_index] # def generate_all_tuples(cur_list): # if len(cur_list[len(cur_list)-1]) - + #clear connectivity data for grp in self.last_mesh.groups: iel_base = grp.element_nr_base - for iel_grp in six.moves.range(grp.nelements): - for i in six.moves.range(len(grp.vertex_indices[iel_grp])): - for j in six.moves.range(i+1, len(grp.vertex_indices[iel_grp])): + for iel_grp in range(grp.nelements): + for i in range(len(grp.vertex_indices[iel_grp])): + for j in range(i+1, len(grp.vertex_indices[iel_grp])): min_vert = min(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j]) max_vert = max(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j]) vertex_pair = (min_vert, max_vert) @@ -629,9 +621,9 @@ class Refiner(object): nelements_in_grp = grp.nelements for grp in groups: - for iel_grp in six.moves.range(len(grp)): + for iel_grp in range(len(grp)): add_verts = [] - for i in six.moves.range(len(grp[iel_grp])): + for i in range(len(grp[iel_grp])): add_verts.append(grp[iel_grp][i]) add_element_to_connectivity(add_verts, new_hanging_vertex_element, iel_base+iel_grp) #assert ray connectivity @@ -652,60 +644,55 @@ class Refiner(object): return self.last_mesh def print_rays(self, ind): - import six - for i in six.moves.range(len(self.last_mesh.groups[0].vertex_indices[ind])): - for j in six.moves.range(i+1, len(self.last_mesh.groups[0].vertex_indices[ind])): + for i in range(len(self.last_mesh.groups[0].vertex_indices[ind])): + for j in range(i+1, len(self.last_mesh.groups[0].vertex_indices[ind])): mn = min(self.last_mesh.groups[0].vertex_indices[ind][i], self.last_mesh.groups[0].vertex_indices[ind][j]) mx = max(self.last_mesh.groups[0].vertex_indices[ind][j], self.last_mesh.groups[0].vertex_indices[ind][i]) vertex_pair = (mn, mx) - print 'LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex - print 'RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex - print 'ADJACENT:' + print('LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex) + print('RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex) + print('ADJACENT:') rays = self.get_leaves(self.pair_map[vertex_pair]) for k in rays: - print k.adjacent_elements + print(k.adjacent_elements) ''' def print_rays(self, groups, ind): - import six - for i in six.moves.range(len(groups[0][ind])): - for j in six.moves.range(i+1, len(groups[0][ind])): + for i in range(len(groups[0][ind])): + for j in range(i+1, len(groups[0][ind])): mn = min(groups[0][ind][i], groups[0][ind][j]) mx = max(groups[0][ind][i], groups[0][ind][j]) vertex_pair = (mn, mx) - print 'LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex - print 'RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex - print 'ADJACENT:' + print('LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex) + print('RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex) + print('ADJACENT:') rays = self.get_leaves(self.pair_map[vertex_pair]) for k in rays: - print k.adjacent_elements + print(k.adjacent_elements) ''' - def print_hanging_elements(self, ind): - import six for i in self.last_mesh.groups[0].vertex_indices[ind]: - print "IND:", i, self.hanging_vertex_element[i] + print("IND:", i, self.hanging_vertex_element[i]) def generate_connectivity(self, nelements, nvertices, groups): # medium-term FIXME: make this an incremental update # rather than build-from-scratch - import six - vertex_to_element = [[] for i in xrange(nvertices)] + vertex_to_element = [[] for i in range(nvertices)] element_index = 0 for grp in groups: - for iel_grp in xrange(len(grp)): + for iel_grp in range(len(grp)): for ivertex in grp[iel_grp]: vertex_to_element[ivertex].append(element_index) element_index += 1 - element_to_element = [set() for i in xrange(nelements)] + element_to_element = [set() for i in range(nelements)] element_index = 0 if self.lazy: for grp in groups: - for iel_grp in six.moves.range(len(grp)): - for i in six.moves.range(len(grp[iel_grp])): - for j in six.moves.range(i+1, len(grp[iel_grp])): + for iel_grp in range(len(grp)): + for i in range(len(grp[iel_grp])): + for j in range(i+1, len(grp[iel_grp])): vertex_pair = (min(grp[iel_grp][i], grp[iel_grp][j]), max(grp[iel_grp][i], grp[iel_grp][j])) #print 'iel:', iel_grp, 'pair:', vertex_pair if vertex_pair not in self.seen_tuple: @@ -714,7 +701,7 @@ class Refiner(object): else: for grp in groups: - for iel_grp in xrange(len(grp)): + for iel_grp in range(len(grp)): for ivertex in grp[iel_grp]: element_to_element[element_index].update( vertex_to_element[ivertex]) @@ -723,8 +710,8 @@ class Refiner(object): if element_index != hanging_element: element_to_element[element_index].update([hanging_element]) element_to_element[hanging_element].update([element_index]) - for i in six.moves.range(len(grp[iel_grp])): - for j in six.moves.range(i+1, len(grp[iel_grp])): + for i in range(len(grp[iel_grp])): + for j in range(i+1, len(grp[iel_grp])): vertex_pair = (min(grp[iel_grp][i], grp[iel_grp][j]), max(grp[iel_grp][i], grp[iel_grp][j])) #element_to_element[element_index].update( @@ -746,11 +733,11 @@ class Refiner(object): element_to_element[self.hanging_vertex_element[ivertex][0]].update([element_index]) ''' element_index += 1 - print len(element_to_element) + print(len(element_to_element)) for iel, neighbors in enumerate(element_to_element): if iel in neighbors: neighbors.remove(iel) - #print self.ray_elements + #print(self.ray_elements) ''' for ray in self.rays: curnode = ray.first diff --git a/meshmode/mesh/tesselate.py b/meshmode/mesh/tesselate.py index b8770ff457cf09f804cbe78a0a59184f751f45fe..4dbfefab948decd297ad00abc303d5cf878f2054 100644 --- a/meshmode/mesh/tesselate.py +++ b/meshmode/mesh/tesselate.py @@ -1,27 +1,29 @@ -''' from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ as gnitstam -node_tuples = list(gnitstam(2, 2)) -node_dict = dict( - (ituple, idx) - for idx, ituple in enumerate(node_tuples)) def add_tuples(a, b): - return tuple(ac+bc for ac, bc in zip(a, b)) - -def try_add_tri(current, d1, d2, d3): - try: - result.append(( - node_dict[add_tuples(current, d1)], - node_dict[add_tuples(current, d2)], - node_dict[add_tuples(current, d3)], - )) - except KeyError: - pass - -result = [] + return tuple(ac+bc for ac, bc in zip(a, b)) + + def tesselatetri(): + result = [] + + node_tuples = list(gnitstam(2, 2)) + node_dict = dict( + (ituple, idx) + for idx, ituple in enumerate(node_tuples)) + + def try_add_tri(current, d1, d2, d3): + try: + result.append(( + node_dict[add_tuples(current, d1)], + node_dict[add_tuples(current, d2)], + node_dict[add_tuples(current, d3)], + )) + except KeyError: + pass + if len(result) > 0: return [node_tuples, result] for current in node_tuples: @@ -32,49 +34,40 @@ def tesselatetri(): try_add_tri(current, (0, 0), (1, 0), (0, 1)) try_add_tri(current, (1, 0), (1, 1), (0, 1)) return [node_tuples, result] -#tesselatetri() -#print node_tuples -#print result -''' -from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - -node_tuples = list(gnitstam(2, 3)) -node_dict = dict( - (ituple, idx) - for idx, ituple in enumerate(node_tuples)) -def add_tuples(a, b): - return tuple(ac+bc for ac, bc in zip(a, b)) - -def try_add_tet(current, d1, d2, d3, d4): - try: - result.append(( - node_dict[add_tuples(current, d1)], - node_dict[add_tuples(current, d2)], - node_dict[add_tuples(current, d3)], - node_dict[add_tuples(current, d4)], - )) - except KeyError: - pass - -result = [] def tesselatetet(): - if len(result) > 0: - return [node_tuples, result] - for current in node_tuples: - # this is a tesselation of a cube into six tets. - # subtets that fall outside of the master tet are simply not added. - - # positively oriented - try_add_tet(current, (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)) - try_add_tet(current, (1, 0, 1), (1, 0, 0), (0, 0, 1), (0, 1, 0)) - try_add_tet(current, (1, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 1)) - - try_add_tet(current, (1, 0, 0), (0, 1, 0), (1, 0, 1), (1, 1, 0)) - try_add_tet(current, (0, 1, 1), (0, 1, 0), (1, 1, 0), (1, 0, 1)) - try_add_tet(current, (0, 1, 1), (1, 1, 1), (1, 0, 1), (1, 1, 0)) - print result - return [node_tuples, result] + node_tuples = list(gnitstam(2, 3)) + + node_dict = dict( + (ituple, idx) + for idx, ituple in enumerate(node_tuples)) + def try_add_tet(current, d1, d2, d3, d4): + try: + result.append(( + node_dict[add_tuples(current, d1)], + node_dict[add_tuples(current, d2)], + node_dict[add_tuples(current, d3)], + node_dict[add_tuples(current, d4)], + )) + except KeyError: + pass + + result = [] + + if len(result) > 0: + return [node_tuples, result] + for current in node_tuples: + # this is a tesselation of a cube into six tets. + # subtets that fall outside of the master tet are simply not added. + + # positively oriented + try_add_tet(current, (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)) + try_add_tet(current, (1, 0, 1), (1, 0, 0), (0, 0, 1), (0, 1, 0)) + try_add_tet(current, (1, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 1)) + + try_add_tet(current, (1, 0, 0), (0, 1, 0), (1, 0, 1), (1, 1, 0)) + try_add_tet(current, (0, 1, 1), (0, 1, 0), (1, 1, 0), (1, 0, 1)) + try_add_tet(current, (0, 1, 1), (1, 1, 1), (1, 0, 1), (1, 1, 0)) + return [node_tuples, result] diff --git a/meshmode/mesh/tools.py b/meshmode/mesh/tools.py index a67a0fb8a310c0b7626a66159b567ae85762759c..4b46bb83d6022710aed5505c3c21646929de8426 100644 --- a/meshmode/mesh/tools.py +++ b/meshmode/mesh/tools.py @@ -28,6 +28,8 @@ import numpy as np from pytools.spatial_btree import SpatialBinaryTreeBucket +# {{{ make_element_lookup_tree + def make_element_lookup_tree(mesh, eps=1e-12): from meshmode.mesh.processing import find_bounding_box bbox_min, bbox_max = find_bounding_box(mesh) @@ -46,3 +48,48 @@ def make_element_lookup_tree(mesh, eps=1e-12): tree.insert((igrp, iel_grp), (el_bbox_min, el_bbox_max)) return tree + +# }}} + + +# {{{ nd_quad_submesh + +def nd_quad_submesh(node_tuples): + """Return a list of tuples of indices into the node list that + generate a tesselation of the reference element. + + :arg node_tuples: A list of tuples *(i, j, ...)* of integers + indicating node positions inside the unit element. The + returned list references indices in this list. + + :func:`pytools.generate_nonnegative_integer_tuples_below` + may be used to generate *node_tuples*. + + See also :func:`modepy.tools.simplex_submesh`. + """ + + from pytools import single_valued, add_tuples + dims = single_valued(len(nt) for nt in node_tuples) + + node_dict = dict( + (ituple, idx) + for idx, ituple in enumerate(node_tuples)) + + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + result = [] + for current in node_tuples: + try: + result.append(tuple( + node_dict[add_tuples(current, offset)] + for offset in gnitb(2, dims))) + + except KeyError: + pass + + return result + +# }}} + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py index 8966861b168f5d740a099e5ec4ad80baccdba4b0..caa9c1beb3cdb7123a3c9798ce04a3c882b14f93 100644 --- a/meshmode/mesh/visualization.py +++ b/meshmode/mesh/visualization.py @@ -30,7 +30,8 @@ import numpy as np # {{{ draw_2d_mesh def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True, - draw_connectivity=False, **kwargs): + draw_nodal_adjacency=False, draw_face_numbers=False, + set_bounding_box=False, **kwargs): assert mesh.ambient_dim == 2 import matplotlib.pyplot as pt @@ -74,7 +75,7 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True, ha="center", va="center", color="blue", bbox=dict(facecolor='white', alpha=0.5, lw=0)) - if draw_connectivity: + if draw_nodal_adjacency: def global_iel_to_group_and_iel(global_iel): for igrp, grp in enumerate(mesh.groups): if global_iel < grp.nelements: @@ -83,7 +84,7 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True, raise ValueError("invalid element nr") - cnx = mesh.element_connectivity + cnx = mesh.nodal_adjacency nb_starts = cnx.neighbors_starts for iel_g in range(mesh.nelements): @@ -110,6 +111,42 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True, length_includes_head=True, color="green", head_width=1e-2, lw=1e-2) + if draw_face_numbers: + for igrp, grp in enumerate(mesh.groups): + for iel, el in enumerate(grp.vertex_indices): + elverts = mesh.vertices[:, el] + el_center = np.mean(elverts, axis=-1) + + for iface, fvi in enumerate(grp.face_vertex_indices()): + face_center = ( + 0.3*el_center + + + 0.7*np.mean(elverts[:, fvi], axis=-1)) + + pt.text(face_center[0], face_center[1], str(iface), fontsize=12, + ha="center", va="center", color="purple", + bbox=dict(facecolor='white', alpha=0.5, lw=0)) + + if set_bounding_box: + from meshmode.mesh.processing import find_bounding_box + lower, upper = find_bounding_box(mesh) + pt.xlim([lower[0], upper[0]]) + pt.ylim([lower[1], upper[1]]) + +# }}} + + +# {{{ draw_curve + +def draw_curve(mesh): + import matplotlib.pyplot as pt + pt.plot(mesh.vertices[0], mesh.vertices[1], "o") + + for i, group in enumerate(mesh.groups): + pt.plot( + group.nodes[0].ravel(), + group.nodes[1].ravel(), "-x", label="Group %d" % i) + # }}} # vim: foldmethod=marker diff --git a/meshmode/version.py b/meshmode/version.py index 92c29262fa4d3e5d7e00f6e6bef07af819a6d872..1da70a1e09bc16f4b7e677505af6cf4dd4502b86 100644 --- a/meshmode/version.py +++ b/meshmode/version.py @@ -1,2 +1,2 @@ -VERSION = (2014, 1) +VERSION = (2015, 1) VERSION_TEXT = ".".join(str(i) for i in VERSION) diff --git a/requirements.txt b/requirements.txt index a6902dd2d88736406f611ddf5d76b324987f5dcb..8e99cd7fdd25e578a5dc11c56f4c43506d18974f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,7 @@ numpy git+git://github.com/inducer/modepy git+git://github.com/pyopencl/pyopencl +git+git://github.com/inducer/islpy git+git://github.com/inducer/loopy # required by pytential, which is in turn needed for some tests diff --git a/setup.cfg b/setup.cfg index 6faef2e65abe138f9ab22c94452c43be0e52c075..512924240ea5f744029e9dd0395769ca66267ca1 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,3 @@ [flake8] -ignore = E126,E127,E128,E123,E226,E241,E242,E265 +ignore = E126,E127,E128,E123,E226,E241,E242,E265,W503,E402 max-line-length=85 diff --git a/test/circle.step b/test/circle.step new file mode 100644 index 0000000000000000000000000000000000000000..78f370557b0fd9f525f484a7d6cf91ba1ea78e29 --- /dev/null +++ b/test/circle.step @@ -0,0 +1,85 @@ +ISO-10303-21; +HEADER; +FILE_DESCRIPTION(('FreeCAD Model'),'2;1'); +FILE_NAME('/home/andreas/own/graphics/freecad/circle.step', + '2015-04-28T14:41:10',('FreeCAD'),('FreeCAD'), + 'Open CASCADE STEP processor 6.7','FreeCAD','Unknown'); +FILE_SCHEMA(('AUTOMOTIVE_DESIGN_CC2 { 1 2 10303 214 -1 1 5 4 }')); +ENDSEC; +DATA; +#1 = APPLICATION_PROTOCOL_DEFINITION('committee draft', + 'automotive_design',1997,#2); +#2 = APPLICATION_CONTEXT( + 'core data for automotive mechanical design processes'); +#3 = SHAPE_DEFINITION_REPRESENTATION(#4,#10); +#4 = PRODUCT_DEFINITION_SHAPE('','',#5); +#5 = PRODUCT_DEFINITION('design','',#6,#9); +#6 = PRODUCT_DEFINITION_FORMATION('','',#7); +#7 = PRODUCT('Circle001','Circle001','',(#8)); +#8 = MECHANICAL_CONTEXT('',#2,'mechanical'); +#9 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design'); +#10 = MANIFOLD_SURFACE_SHAPE_REPRESENTATION('',(#11,#15),#46); +#11 = AXIS2_PLACEMENT_3D('',#12,#13,#14); +#12 = CARTESIAN_POINT('',(0.,0.,0.)); +#13 = DIRECTION('',(0.,0.,1.)); +#14 = DIRECTION('',(1.,0.,-0.)); +#15 = SHELL_BASED_SURFACE_MODEL('',(#16)); +#16 = OPEN_SHELL('',(#17)); +#17 = ADVANCED_FACE('',(#18),#31,.T.); +#18 = FACE_BOUND('',#19,.F.); +#19 = EDGE_LOOP('',(#20)); +#20 = ORIENTED_EDGE('',*,*,#21,.T.); +#21 = EDGE_CURVE('',#22,#22,#24,.T.); +#22 = VERTEX_POINT('',#23); +#23 = CARTESIAN_POINT('',(2.E-04,0.,0.)); +#24 = SURFACE_CURVE('',#25,(#30),.PCURVE_S1.); +#25 = CIRCLE('',#26,2.E-04); +#26 = AXIS2_PLACEMENT_3D('',#27,#28,#29); +#27 = CARTESIAN_POINT('',(0.,0.,0.)); +#28 = DIRECTION('',(0.,0.,1.)); +#29 = DIRECTION('',(1.,0.,0.)); +#30 = PCURVE('',#31,#36); +#31 = PLANE('',#32); +#32 = AXIS2_PLACEMENT_3D('',#33,#34,#35); +#33 = CARTESIAN_POINT('',(2.E-04,0.,0.)); +#34 = DIRECTION('',(0.,0.,-1.)); +#35 = DIRECTION('',(-1.,0.,0.)); +#36 = DEFINITIONAL_REPRESENTATION('',(#37),#45); +#37 = ( BOUNDED_CURVE() B_SPLINE_CURVE(2,(#38,#39,#40,#41,#42,#43,#44), +.UNSPECIFIED.,.F.,.F.) B_SPLINE_CURVE_WITH_KNOTS((1,2,2,2,2,1),( + -2.094395102393,0.,2.094395102393,4.188790204786,6.28318530718, +8.377580409573),.UNSPECIFIED.) CURVE() GEOMETRIC_REPRESENTATION_ITEM() +RATIONAL_B_SPLINE_CURVE((1.,0.5,1.,0.5,1.,0.5,1.)) REPRESENTATION_ITEM( + '') ); +#38 = CARTESIAN_POINT('',(0.,0.)); +#39 = CARTESIAN_POINT('',(0.,3.464101615138E-04)); +#40 = CARTESIAN_POINT('',(3.E-04,1.732050807569E-04)); +#41 = CARTESIAN_POINT('',(6.E-04,4.898587196589E-20)); +#42 = CARTESIAN_POINT('',(3.E-04,-1.732050807569E-04)); +#43 = CARTESIAN_POINT('',(2.981555974335E-19,-3.464101615138E-04)); +#44 = CARTESIAN_POINT('',(0.,0.)); +#45 = ( GEOMETRIC_REPRESENTATION_CONTEXT(2) +PARAMETRIC_REPRESENTATION_CONTEXT() REPRESENTATION_CONTEXT('2D SPACE','' + ) ); +#46 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3) +GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#50)) GLOBAL_UNIT_ASSIGNED_CONTEXT( +(#47,#48,#49)) REPRESENTATION_CONTEXT('Context #1', + '3D Context with UNIT and UNCERTAINTY') ); +#47 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT($,.METRE.) ); +#48 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) ); +#49 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() ); +#50 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-10),#47, + 'distance_accuracy_value','confusion accuracy'); +#51 = PRODUCT_TYPE('part',$,(#7)); +#52 = MECHANICAL_DESIGN_GEOMETRIC_PRESENTATION_REPRESENTATION('',(#53), + #46); +#53 = STYLED_ITEM('color',(#54),#17); +#54 = PRESENTATION_STYLE_ASSIGNMENT((#55)); +#55 = SURFACE_STYLE_USAGE(.BOTH.,#56); +#56 = SURFACE_SIDE_STYLE('',(#57)); +#57 = SURFACE_STYLE_FILL_AREA(#58); +#58 = FILL_AREA_STYLE('',(#59)); +#59 = FILL_AREA_STYLE_COLOUR('',#60); +#60 = COLOUR_RGB('',0.800000011921,0.800000011921,0.800000011921); +ENDSEC; +END-ISO-10303-21; diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cacfa3b896691340ba8b29ac081c4036dbc14e4b..c7c2bc488932feec285147eaf2fc7fa7358deab4 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1,5 +1,4 @@ from __future__ import division, absolute_import, print_function -from six.moves import range __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" @@ -23,6 +22,7 @@ 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 import numpy.linalg as la import pyopencl as cl @@ -33,59 +33,343 @@ from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory, + PolynomialWarpAndBlendGroupFactory, + PolynomialEquidistantGroupFactory, + ) +from meshmode.mesh import BTAG_ALL +from meshmode.discretization.connection import \ + FRESTR_ALL_FACES, FRESTR_INTERIOR_FACES + import pytest import logging logger = logging.getLogger(__name__) -def test_boundary_interpolation(ctx_getter): +# {{{ circle mesh + +def test_circle_mesh(do_plot=False): + from meshmode.mesh.io import generate_gmsh, FileSource + print("BEGIN GEN") + mesh = generate_gmsh( + FileSource("circle.step"), 2, order=2, + force_ambient_dim=2, + other_options=[ + "-string", "Mesh.CharacteristicLengthMax = 0.05;"] + ) + print("END GEN") + print(mesh.nelements) + + from meshmode.mesh.processing import affine_map + mesh = affine_map(mesh, A=3*np.eye(2)) + + if do_plot: + from meshmode.mesh.visualization import draw_2d_mesh + draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True, + set_bounding_box=True) + import matplotlib.pyplot as pt + pt.show() + +# }}} + + +# {{{ convergence of boundary interpolation + +@pytest.mark.parametrize("group_factory", [ + InterpolatoryQuadratureSimplexGroupFactory, + PolynomialWarpAndBlendGroupFactory + ]) +@pytest.mark.parametrize("boundary_tag", [ + BTAG_ALL, + FRESTR_ALL_FACES, + FRESTR_INTERIOR_FACES, + ]) +@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_boundary_interpolation(ctx_getter, group_factory, boundary_tag, + mesh_name, dim, mesh_pars, per_face_groups): cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) - from meshmode.mesh.io import generate_gmsh, FileSource from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory - from meshmode.discretization.connection import make_boundary_restriction + from meshmode.discretization.connection import ( + make_face_restriction, check_connection) from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() order = 4 - for h in [1e-1, 3e-2, 1e-2]: - print("BEGIN GEN") - mesh = generate_gmsh( - FileSource("blob-2d.step"), 2, order=order, - force_ambient_dim=2, - other_options=[ - "-string", "Mesh.CharacteristicLengthMax = %s;" % h] - ) - print("END GEN") + def f(x): + return 0.1*cl.clmath.sin(30*x) + + for mesh_par in mesh_pars: + # {{{ get mesh + + if mesh_name == "blob": + assert dim == 2 + + h = mesh_par + + from meshmode.mesh.io import generate_gmsh, FileSource + print("BEGIN GEN") + mesh = generate_gmsh( + FileSource("blob-2d.step"), 2, order=order, + force_ambient_dim=2, + other_options=[ + "-string", "Mesh.CharacteristicLengthMax = %s;" % h] + ) + 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) + + h = 1/mesh_par + else: + raise ValueError("mesh_name not recognized") + + # }}} vol_discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) + group_factory(order)) print("h=%s -> %d elements" % ( h, sum(mgrp.nelements for mgrp in mesh.groups))) x = vol_discr.nodes()[0].with_queue(queue) - f = 0.1*cl.clmath.sin(30*x) + vol_f = f(x) + + bdry_connection = make_face_restriction( + vol_discr, group_factory(order), + boundary_tag, per_face_groups=per_face_groups) + check_connection(bdry_connection) + bdry_discr = bdry_connection.to_discr + + bdry_x = bdry_discr.nodes()[0].with_queue(queue) + bdry_f = f(bdry_x) + bdry_f_2 = bdry_connection(queue, vol_f) + + if mesh_name == "blob" and dim == 2: + mat = bdry_connection.full_resample_matrix(queue).get(queue) + bdry_f_2_by_mat = mat.dot(vol_f.get()) + + mat_error = la.norm(bdry_f_2.get(queue=queue) - bdry_f_2_by_mat) + assert mat_error < 1e-14, mat_error + + err = la.norm((bdry_f-bdry_f_2).get(), np.inf) + eoc_rec.add_data_point(h, err) + + print(eoc_rec) + assert ( + eoc_rec.order_estimate() >= order-0.5 + or eoc_rec.max_error() < 1e-14) + +# }}} + + +# {{{ boundary-to-all-faces connecttion + +@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(ctx_getter, mesh_name, dim, mesh_pars, + per_face_groups): + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + from meshmode.discretization import Discretization + from meshmode.discretization.connection import ( + make_face_restriction, make_face_to_all_faces_embedding, + check_connection) + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + order = 4 + + def f(x): + return 0.1*cl.clmath.sin(30*x) + + for mesh_par in mesh_pars: + # {{{ get mesh + + if mesh_name == "blob": + assert dim == 2 + + h = mesh_par + + from meshmode.mesh.io import generate_gmsh, FileSource + print("BEGIN GEN") + mesh = generate_gmsh( + FileSource("blob-2d.step"), 2, order=order, + force_ambient_dim=2, + other_options=[ + "-string", "Mesh.CharacteristicLengthMax = %s;" % h] + ) + 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) + + h = 1/mesh_par + else: + raise ValueError("mesh_name not recognized") + + # }}} + + vol_discr = Discretization(cl_ctx, mesh, + PolynomialWarpAndBlendGroupFactory(order)) + print("h=%s -> %d elements" % ( + h, sum(mgrp.nelements for mgrp in mesh.groups))) + + all_face_bdry_connection = make_face_restriction( + vol_discr, PolynomialWarpAndBlendGroupFactory(order), + FRESTR_ALL_FACES, per_face_groups=per_face_groups) + all_face_bdry_discr = all_face_bdry_connection.to_discr + + for ito_grp, ceg in enumerate(all_face_bdry_connection.groups): + for ibatch, batch in enumerate(ceg.batches): + assert np.array_equal( + batch.from_element_indices.get(queue), + np.arange(vol_discr.mesh.nelements)) + + if per_face_groups: + assert ito_grp == batch.to_element_face + else: + assert ibatch == batch.to_element_face + + all_face_x = all_face_bdry_discr.nodes()[0].with_queue(queue) + all_face_f = f(all_face_x) + + all_face_f_2 = all_face_bdry_discr.zeros(queue) + + for boundary_tag in [ + BTAG_ALL, + FRESTR_INTERIOR_FACES, + ]: + bdry_connection = make_face_restriction( + vol_discr, PolynomialWarpAndBlendGroupFactory(order), + boundary_tag, per_face_groups=per_face_groups) + bdry_discr = bdry_connection.to_discr + + bdry_x = bdry_discr.nodes()[0].with_queue(queue) + bdry_f = f(bdry_x) + + all_face_embedding = make_face_to_all_faces_embedding( + bdry_connection, all_face_bdry_discr) + + check_connection(all_face_embedding) + + all_face_f_2 += all_face_embedding(queue, bdry_f) + + err = la.norm((all_face_f-all_face_f_2).get(), np.inf) + eoc_rec.add_data_point(h, err) + + print(eoc_rec) + assert ( + eoc_rec.order_estimate() >= order-0.5 + or eoc_rec.max_error() < 1e-14) + +# }}} + + +# {{{ convergence of opposite-face interpolation + +@pytest.mark.parametrize("group_factory", [ + InterpolatoryQuadratureSimplexGroupFactory, + PolynomialWarpAndBlendGroupFactory + ]) +@pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ + ("blob", 2, [1e-1, 8e-2, 5e-2]), + ("warp", 2, [3, 5, 7]), + ("warp", 3, [3, 5]), + ]) +def test_opposite_face_interpolation(ctx_getter, group_factory, + mesh_name, dim, mesh_pars): + logging.basicConfig(level=logging.INFO) + + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) - bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction( - queue, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(order)) + from meshmode.discretization import Discretization + from meshmode.discretization.connection import ( + make_face_restriction, make_opposite_face_connection, + check_connection) + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + order = 5 + + def f(x): + return 0.1*cl.clmath.sin(30*x) + + for mesh_par in mesh_pars: + # {{{ get mesh + + if mesh_name == "blob": + assert dim == 2 + + h = mesh_par + + from meshmode.mesh.io import generate_gmsh, FileSource + print("BEGIN GEN") + mesh = generate_gmsh( + FileSource("blob-2d.step"), 2, order=order, + force_ambient_dim=2, + other_options=[ + "-string", "Mesh.CharacteristicLengthMax = %s;" % h] + ) + 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) + + h = 1/mesh_par + else: + raise ValueError("mesh_name not recognized") + + # }}} + + vol_discr = Discretization(cl_ctx, mesh, + group_factory(order)) + print("h=%s -> %d elements" % ( + h, sum(mgrp.nelements for mgrp in mesh.groups))) + + bdry_connection = make_face_restriction( + vol_discr, group_factory(order), + FRESTR_INTERIOR_FACES) + bdry_discr = bdry_connection.to_discr + + opp_face = make_opposite_face_connection(bdry_connection) + check_connection(opp_face) bdry_x = bdry_discr.nodes()[0].with_queue(queue) - bdry_f = 0.1*cl.clmath.sin(30*bdry_x) - bdry_f_2 = bdry_connection(queue, f) + bdry_f = f(bdry_x) + + bdry_f_2 = opp_face(queue, bdry_f) err = la.norm((bdry_f-bdry_f_2).get(), np.inf) eoc_rec.add_data_point(h, err) print(eoc_rec) - assert eoc_rec.order_estimate() >= order-0.5 + assert ( + eoc_rec.order_estimate() >= order-0.5 + or eoc_rec.max_error() < 1e-13) + +# }}} +# {{{ element orientation + def test_element_orientation(): from meshmode.mesh.io import generate_gmsh, FileSource @@ -114,6 +398,10 @@ def test_element_orientation(): assert ((mesh_orient < 0) == (flippy > 0)).all() +# }}} + + +# {{{ merge and map def test_merge_and_map(ctx_getter, visualize=False): from meshmode.mesh.io import generate_gmsh, FileSource @@ -126,10 +414,10 @@ def test_merge_and_map(ctx_getter, visualize=False): other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"] ) - from meshmode.mesh.processing import merge_dijsoint_meshes, affine_map + from meshmode.mesh.processing import merge_disjoint_meshes, affine_map mesh2 = affine_map(mesh, A=np.eye(2), b=np.array([5, 0])) - mesh3 = merge_dijsoint_meshes((mesh2, mesh)) + mesh3 = merge_disjoint_meshes((mesh2, mesh)) if visualize: from meshmode.discretization import Discretization @@ -145,6 +433,10 @@ def test_merge_and_map(ctx_getter, visualize=False): vis = make_visualizer(queue, discr, 1) vis.write_vtk_file("merged.vtu", []) +# }}} + + +# {{{ sanity checks: single element @pytest.mark.parametrize("dim", [2, 3]) @pytest.mark.parametrize("order", [1, 3]) @@ -154,21 +446,21 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False): cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) - from modepy.tools import UNIT_VERTICES - vertices = UNIT_VERTICES[dim].T.copy() + from modepy.tools import unit_vertices + vertices = unit_vertices(dim).T.copy() center = np.empty(dim, np.float64) center.fill(-0.5) import modepy as mp - from meshmode.mesh import SimplexElementGroup, Mesh + 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) - mesh = Mesh(vertices, [mg]) + mesh = Mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None) from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ @@ -196,9 +488,11 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False): # {{{ boundary discretization - from meshmode.discretization.connection import make_boundary_restriction - bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction( - queue, vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3)) + from meshmode.discretization.connection import make_face_restriction + bdry_connection = make_face_restriction( + vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3), + BTAG_ALL) + bdry_discr = bdry_connection.to_discr # }}} @@ -227,6 +521,60 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False): assert normal_outward_check.get().all(), normal_outward_check.get() +# }}} + + +# {{{ sanity check: volume interpolation on scipy/qhull delaunay meshes in nD + +@pytest.mark.parametrize("dim", [2, 3, 4]) +@pytest.mark.parametrize("order", [3]) +def test_sanity_qhull_nd(ctx_getter, dim, order): + pytest.importorskip("scipy") + + logging.basicConfig(level=logging.INFO) + + ctx = ctx_getter() + queue = cl.CommandQueue(ctx) + + from scipy.spatial import Delaunay + verts = np.random.rand(1000, dim) + dtri = Delaunay(verts) + + from meshmode.mesh.io import from_vertices_and_simplices + mesh = from_vertices_and_simplices(dtri.points.T, dtri.simplices, + fix_orientation=True) + + from meshmode.discretization import Discretization + low_discr = Discretization(ctx, mesh, + PolynomialEquidistantGroupFactory(order)) + high_discr = Discretization(ctx, mesh, + PolynomialEquidistantGroupFactory(order+1)) + + from meshmode.discretization.connection import make_same_mesh_connection + cnx = make_same_mesh_connection(high_discr, low_discr) + + def f(x): + return 0.1*cl.clmath.sin(x) + + x_low = low_discr.nodes()[0].with_queue(queue) + f_low = f(x_low) + + x_high = high_discr.nodes()[0].with_queue(queue) + f_high_ref = f(x_high) + + f_high_num = cnx(queue, f_low) + + err = (f_high_ref-f_high_num).get() + + err = la.norm(err, np.inf)/la.norm(f_high_ref.get(), np.inf) + + print(err) + assert err < 1e-2 + +# }}} + + +# {{{ sanity checks: ball meshes # python test_meshmode.py 'test_sanity_balls(cl._csc, "disk-radius-1.step", 2, 2, visualize=True)' # noqa @pytest.mark.parametrize(("src_file", "dim"), [ @@ -264,15 +612,15 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order, # {{{ discretizations and connections from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory vol_discr = Discretization(ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(quad_order)) - from meshmode.discretization.connection import make_boundary_restriction - bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction( - queue, vol_discr, - InterpolatoryQuadratureSimplexGroupFactory(quad_order)) + from meshmode.discretization.connection import make_face_restriction + bdry_connection = make_face_restriction( + vol_discr, + InterpolatoryQuadratureSimplexGroupFactory(quad_order), + BTAG_ALL) + bdry_discr = bdry_connection.to_discr # }}} @@ -342,6 +690,10 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order, print(surf_eoc_rec) assert surf_eoc_rec.order_estimate() >= mesh_order +# }}} + + +# {{{ rect/box mesh generation def test_rect_mesh(do_plot=False): from meshmode.mesh.generation import generate_regular_rect_mesh @@ -349,21 +701,41 @@ def test_rect_mesh(do_plot=False): if do_plot: from meshmode.mesh.visualization import draw_2d_mesh - draw_2d_mesh(mesh, fill=None, draw_connectivity=True) + draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True) import matplotlib.pyplot as pt pt.show() -def test_box_mesh(): +def test_box_mesh(ctx_getter, visualize=False): from meshmode.mesh.generation import generate_box_mesh - generate_box_mesh(3*(np.linspace(0, 1, 5),)) + mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),)) + 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, mesh, + PolynomialWarpAndBlendGroupFactory(1)) + + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(queue, discr, 1) + vis.write_vtk_file("box.vtu", []) + +# }}} + + +# {{{ as_python stringification def test_as_python(): - from meshmode.mesh.generation import make_curve_mesh, cloverleaf - mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, 100), order=3) + from meshmode.mesh.generation import generate_box_mesh + mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),)) - mesh.element_connectivity + # These implicitly compute these adjacency structures. + mesh.nodal_adjacency + mesh.facial_adjacency_groups from meshmode.mesh import as_python code = as_python(mesh) @@ -376,6 +748,10 @@ def test_as_python(): assert mesh == mesh_2 +# }}} + + +# {{{ test lookup tree for element finding def test_lookup_tree(do_plot=False): from meshmode.mesh.generation import make_curve_mesh, cloverleaf @@ -399,8 +775,12 @@ def test_lookup_tree(do_plot=False): with open("tree.dat", "w") as outf: tree.visualize(outf) -#construct vertex vertex_index -def get_vertex(mesh, vertex_index): +# }}} + + +# {{{ test refiner + +def _get_vertex(mesh, vertex_index): vertex = np.empty([len(mesh.vertices)]) for i_cur_dim, cur_dim_coords in enumerate(mesh.vertices): vertex[i_cur_dim] = cur_dim_coords[vertex_index] @@ -431,7 +811,7 @@ def get_some_mesh(): from meshmode.mesh.refinement import Refiner r = Refiner(mesh) times = random.randint(1, 2) - for time in xrange(times): + for time in range(times): flags = np.zeros(len(mesh.groups[0].vertex_indices)) ''' if time == 0: @@ -445,15 +825,18 @@ def get_some_mesh(): #flags[8] = 1 #flags[40] = 1 ''' - for i in xrange(0, len(flags)): + for i in range(0, len(flags)): flags[i] = random.randint(0, 1) mesh = r.refine(flags) + return mesh + @pytest.mark.parametrize("mesh_factory", [get_some_mesh]) -@pytest.mark.parametrize("num_rounds", [1,2,3]) +@pytest.mark.parametrize("num_rounds", [1, 2, 3]) def test_refiner_connectivity(mesh_factory, num_rounds): mesh = mesh_factory() + def group_and_iel_to_global_iel(igrp, iel): return mesh.groups[igrp].element_nr_base + iel @@ -463,17 +846,17 @@ def test_refiner_connectivity(mesh_factory, num_rounds): from meshmode.mesh.processing import find_bounding_box bbox_min, bbox_max = find_bounding_box(mesh) - extent = bbox_max-bbox_min cnx = mesh.element_connectivity nvertices_per_element = len(mesh.groups[0].vertex_indices[0]) #stores connectivity obtained from geometry using barycentric coordinates - connected_to_element_geometry = np.zeros([mesh.nelements, mesh.nelements], dtype=bool) + connected_to_element_geometry = np.zeros( + (mesh.nelements, mesh.nelements), dtype=bool) #store connectivity obtained from mesh's connectivity - connected_to_element_connectivity = np.zeros([mesh.nelements, mesh.nelements], dtype=bool) + connected_to_element_connectivity = np.zeros( + (mesh.nelements, mesh.nelements), dtype=bool) import six for igrp, grp in enumerate(mesh.groups): - iel_base = grp.element_nr_base for iel_grp in six.moves.range(grp.nelements): iel_g = group_and_iel_to_global_iel(igrp, iel_grp) @@ -481,23 +864,26 @@ def test_refiner_connectivity(mesh_factory, num_rounds): for nb_iel_g in cnx.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]: connected_to_element_connectivity[nb_iel_g, iel_g] = True connected_to_element_connectivity[iel_g, nb_iel_g] = True - + for vertex_index in grp.vertex_indices[iel_grp]: - vertex = get_vertex(mesh, vertex_index) - #check which elements touch this vertex + vertex = _get_vertex(mesh, vertex_index) + + # check which elements touch this vertex for bounding_igrp, bounding_iel in tree.generate_matches(vertex): if bounding_igrp == igrp and bounding_iel == iel_grp: continue bounding_grp = mesh.groups[bounding_igrp] - last_bounding_vertex = get_vertex(mesh, \ + last_bounding_vertex = _get_vertex(mesh, bounding_grp.vertex_indices[bounding_iel][nvertices_per_element-1]) + transformation = np.empty([len(mesh.vertices), nvertices_per_element-1]) vertex_transformed = vertex - last_bounding_vertex - for ibounding_vertex_index, bounding_vertex_index in \ - enumerate(bounding_grp.vertex_indices[bounding_iel][:nvertices_per_element-1]): - bounding_vertex = get_vertex(mesh, bounding_vertex_index) - transformation[:,ibounding_vertex_index] = bounding_vertex - last_bounding_vertex + for ibounding_vertex_index, bounding_vertex_index in enumerate( + bounding_grp.vertex_indices[bounding_iel][:nvertices_per_element-1]): + bounding_vertex = _get_vertex(mesh, bounding_vertex_index) + transformation[:, ibounding_vertex_index] = \ + bounding_vertex - last_bounding_vertex barycentric_coordinates, resid, _, _ = np.linalg.lstsq(transformation, vertex_transformed) bc = barycentric_coordinates @@ -507,9 +893,11 @@ def test_refiner_connectivity(mesh_factory, num_rounds): continue if ((bc > -tol).all() and np.sum(bc) < 1+tol): - connected_to_element_geometry[group_and_iel_to_global_iel(bounding_igrp, bounding_iel),\ + connected_to_element_geometry[ + group_and_iel_to_global_iel(bounding_igrp, bounding_iel), group_and_iel_to_global_iel(igrp, iel_grp)] = True - connected_to_element_geometry[group_and_iel_to_global_iel(igrp, iel_grp),\ + connected_to_element_geometry[ + group_and_iel_to_global_iel(igrp, iel_grp), group_and_iel_to_global_iel(bounding_igrp, bounding_iel)] = True ''' @@ -530,6 +918,33 @@ def test_refiner_connectivity(mesh_factory, num_rounds): ''' assert (connected_to_element_geometry == connected_to_element_connectivity).all() +# }}} + + +# {{{ test_nd_quad_submesh + +@pytest.mark.parametrize("dims", [2, 3, 4]) +def test_nd_quad_submesh(dims): + from meshmode.mesh.tools import nd_quad_submesh + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + node_tuples = list(gnitb(3, dims)) + + for i, nt in enumerate(node_tuples): + print(i, nt) + + assert len(node_tuples) == 3**dims + + elements = nd_quad_submesh(node_tuples) + + for e in elements: + print(e) + + assert len(elements) == 2**dims + +# }}} + + if __name__ == "__main__": import sys if len(sys.argv) > 1: