Towards basic functionality

doc/discretization.rst

@@ -12,6 +12,11 @@ Composite polynomial discretization
 .. automodule:: meshmode.discretization.poly_element
 
 Resampling
------------------------------------
+----------
 
 .. automodule:: meshmode.discretization.resampling
+
+Visualization
+-------------
+
+.. automodule:: meshmode.discretization.visualization

doc/index.rst

@@ -5,6 +5,7 @@ Contents:
 
 .. toctree::
    :maxdepth: 2
+
    mesh
    discretization
 

doc/mesh.rst

@@ -16,4 +16,14 @@ Mesh generation
 
 .. automodule:: pytential.mesh.generation
 
+Mesh input/output
+-----------------
+
+.. automodule:: pytential.mesh.io
+
+Mesh connections (for interpolation)
+------------------------------------
+
+.. automodule:: pytential.mesh.connection
+
 .. vim: sw=4

meshmode/discretization/connection.py

+from __future__ import division
+
+__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
+
+
+__doc__ = """
+.. autoclass:: DiscretizationConnection
+
+.. autofunction:: make_same_mesh_connection
+
+Implementation details
+^^^^^^^^^^^^^^^^^^^^^^
+
+.. autoclass:: InterpolationGroup
+
+.. autoclass:: DiscretizationConnectionElementGroup
+"""
+
+
+class InterpolationGroup(object):
+    """
+    .. attribute:: source_element_indices
+
+        A :class:`numpy.ndarray` of length ``nelements``, containing the
+        element index from which this "*to*" element's data will be
+        interpolated.
+
+    .. attribute:: target_element_indices
+
+        A :class:`numpy.ndarray` of length ``nelements``, containing the
+        element index to which this "*to*" element's data will be
+        interpolated.
+
+    .. attribute:: source_interpolation_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.
+    """
+    def __init__(self, source_element_indices,
+            target_element_indices, source_interpolation_nodes):
+        self.source_element_indices = source_element_indices
+        self.target_element_indices = target_element_indices
+        self.source_interpolation_nodes = source_interpolation_nodes
+
+    @property
+    def nelements(self):
+        return len(self.source_element_indices)
+
+
+class DiscretizationConnectionElementGroup(object):
+    """
+    .. attribute:: interpolation_groups
+
+        A list of :class:`InterpolationGroup` instances.
+    """
+    def __init__(self, interpolation_groups):
+        self.interpolation_groups = interpolation_groups
+
+
+class DiscretizationConnection(object):
+    """
+    .. attribute:: from_discr
+
+    .. attribute:: to_discr
+
+    .. attribute:: groups
+
+        a list of :class:`MeshConnectionGroup` instances, with
+        a one-to-one correspondence to the groups in
+        :attr:`from_discr` and :attr:`to_discr`.
+    """
+
+    def __init__(self, from_discr, to_discr, groups):
+        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
+
+        self.from_discr = from_discr
+        self.to_discr = to_discr
+        self.groups = groups
+
+    def __call__(self, queue, field):
+        pass
+
+
+# {{{ constructor functions
+
+def make_same_mesh_connection(queue, from_discr, to_discr):
+    if from_discr.mesh is not to_discr.mesh:
+        raise ValueError("from_discr and to_discr must be based on "
+                "the same mesh")
+
+    assert queue.context == from_discr.cl_context
+    assert queue.context == to_discr.cl_context
+
+    groups = []
+    for fgrp, tgrp in zip(from_discr.groups, to_discr.groups):
+        all_elements = cl.array.arange(queue,
+                fgrp.nelements,
+                dtype=np.intp)
+        igroup = InterpolationGroup(
+                source_element_indices=all_elements,
+                target_element_indices=all_elements,
+                source_interpolation_nodes=tgrp.unit_nodes)
+
+        groups.append(
+                DiscretizationConnectionElementGroup(
+                    igroup))
+
+    return DiscretizationConnection(
+            from_discr, to_discr, groups)
+
+# }}}
+
+# vim: foldmethod=marker

meshmode/discretization/poly_element.py

@@ -107,16 +107,16 @@ class PolynomialElementGroupBase(object):
 # }}}
 
 
-# {{{ element group
+# {{{ concrete element groups
 
-class PolynomialElementGroup(PolynomialElementGroupBase):
+class PolynomialQuadratureElementGroup(PolynomialElementGroupBase):
     @memoize_method
     def _quadrature_rule(self):
         dims = self.mesh_el_group.dim
         if dims == 1:
             return mp.LegendreGaussQuadrature(self.order)
         else:
-            return mp.VioreanuRokhlinSimplexQuad(self.order, dims)
+            return mp.VioreanuRokhlinSimplexQuadrature(self.order, dims)
 
     @property
     @memoize_method
@@ -128,12 +128,25 @@ class PolynomialElementGroup(PolynomialElementGroupBase):
     def weights(self):
         return self._quadrature_rule().weights
 
+
+class PolynomialWarpAndBlendElementGroup(PolynomialElementGroupBase):
+    @property
+    @memoize_method
+    def unit_nodes(self):
+        dims = self.mesh_el_group.dim
+        return mp.warp_and_blend_nodes(dims, self.order)
+
+    @property
+    @memoize_method
+    def weights(self):
+        raise NotImplementedError()
+
 # }}}
 
 
 # {{{ discretization
 
-class PolynomialElementDiscretization(Discretization):
+class PolynomialElementDiscretizationBase(Discretization):
     """An (unstructured) composite polynomial discretization without
     any specific opinion on how to evaluate layer potentials.
 
@@ -280,9 +293,16 @@ class PolynomialElementDiscretization(Discretization):
 
         return result
 
-    group_class = PolynomialElementGroup
-
 # }}}
 
 
+class PolynomialQuadratureElementDiscretization(
+        PolynomialElementDiscretizationBase):
+    group_class = PolynomialQuadratureElementGroup
+
+
+class PolynomialWarpAndBlendElementDiscretization(
+        PolynomialElementDiscretizationBase):
+    group_class = PolynomialWarpAndBlendElementGroup
+
 # vim: fdm=marker

meshmode/discretization/visualization.py

+from __future__ import division
+
+__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
+from pytools import memoize_method
+import pyopencl as cl
+
+__doc__ = """
+
+.. autofunction:: make_visualizer
+
+.. autoclass:: Visualizer
+"""
+
+
+def separate_by_real_and_imag(data, real_only):
+    for name, field in data:
+        from pytools.obj_array import log_shape
+        ls = log_shape(field)
+
+        if ls != () and ls[0] > 1:
+            assert len(ls) == 1
+            from pytools.obj_array import (
+                    oarray_real_copy, oarray_imag_copy,
+                    with_object_array_or_scalar)
+
+            if field[0].dtype.kind == "c":
+                if real_only:
+                    yield (name,
+                            with_object_array_or_scalar(oarray_real_copy, field))
+                else:
+                    yield (name+"_r",
+                            with_object_array_or_scalar(oarray_real_copy, field))
+                    yield (name+"_i",
+                            with_object_array_or_scalar(oarray_imag_copy, field))
+            else:
+                yield (name, field)
+        else:
+            # ls == ()
+            if field.dtype.kind == "c":
+                yield (name+"_r", field.real.copy())
+                yield (name+"_i", field.imag.copy())
+            else:
+                yield (name, field)
+
+
+class Visualizer(object):
+    def __init__(self, discr, vis_discr, connection):
+        self.discr = discr
+        self.vis_discr = vis_discr
+        self.connection = connection
+
+    @memoize_method
+    def _vis_connectivity(self):
+        """
+        :return: an array of shape
+            ``(vis_discr.nelements,nsubelements,primitive_element_size)``
+        """
+        # Assume that we're using modepy's default node ordering.
+
+        from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \
+                as gnitstam, single_valued
+        vis_order = single_valued(
+                group.order for group in self.vis_discr.groups)
+        node_tuples = list(gnitstam(vis_order, self.vis_discr.dim))
+
+        from modepy.tools import submesh
+        el_connectivity = np.array(
+                submesh(node_tuples),
+                dtype=np.intp)
+
+        nelements = sum(group.nelements for group in self.vis_discr.groups)
+        vis_connectivity = np.empty(
+                (nelements,) + el_connectivity.shape, dtype=np.intp)
+
+        el_nr_base = 0
+        for group in self.vis_discr.groups:
+            assert len(node_tuples) == group.nunit_nodes
+            vis_connectivity[el_nr_base:el_nr_base+group.nelements] = (
+                    np.arange(
+                        el_nr_base*group.nunit_nodes,
+                        (el_nr_base+group.nelements)*group.nunit_nodes,
+                        group.nunit_nodes
+                        )[:, np.newaxis, np.newaxis]
+                    + el_connectivity)
+            el_nr_base += group.nelements
+
+        return vis_connectivity
+
+    def show_scalar_in_mayavi(self, field, **kwargs):
+        import mayavi.mlab as mlab
+
+        do_show = kwargs.pop("do_show", True)
+
+        with cl.CommandQueue(self.vis_discr.cl_context) as queue:
+            nodes = self.vis_discr.nodes().with_queue(queue).get()
+
+        assert nodes.shape[0] == self.vis_discr.ambient_dim
+        #mlab.points3d(nodes[0], nodes[1], 0*nodes[0])
+
+        vis_connectivity = self._vis_connectivity()
+
+        if self.vis_discr.dim == 2:
+            nodes = list(nodes)
+            # pad to 3D with zeros
+            while len(nodes) < 3:
+                nodes.append(0*nodes[0])
+
+            args = tuple(nodes) + (vis_connectivity.reshape(-1, 3),)
+
+            mlab.triangular_mesh(*args, **kwargs)
+
+        else:
+            raise RuntimeError("meshes of bulk dimension %d are currently "
+                    "unsupported" % self.vis_discr.dim)
+
+        if do_show:
+            mlab.show()
+
+    def write_vtk_file(self, file_name, fields_and_names, compressor=None,
+            real_only=False):
+
+        from pyvisfile.vtk import (
+                UnstructuredGrid, DataArray,
+                AppendedDataXMLGenerator,
+                VTK_LINE, VTK_TRIANGLE, VTK_TETRA,
+                VF_LIST_OF_COMPONENTS)
+        el_types = {
+                1: VTK_LINE,
+                2: VTK_TRIANGLE,
+                3: VTK_TETRA,
+                }
+
+        el_type = el_types[self.vis_discr.dim]
+
+        with cl.CommandQueue(self.vis_discr.cl_context) as queue:
+            nodes = self.vis_discr.nodes().with_queue(queue).get()
+
+        connectivity = self._vis_connectivity()
+
+        nprimitive_elements = (
+                connectivity.shape[0]
+                * connectivity.shape[1])
+
+        grid = UnstructuredGrid(
+                (self.vis_discr.nnodes,
+                    DataArray("points",
+                        nodes.reshape(self.vis_discr.ambient_dim, -1),
+                        vector_format=VF_LIST_OF_COMPONENTS)),
+                cells=connectivity.reshape(-1),
+                cell_types=np.asarray([el_type] * nprimitive_elements,
+                    dtype=np.uint8))
+
+        # for name, field in separate_by_real_and_imag(cell_data, real_only):
+        #     grid.add_celldata(DataArray(name, field,
+        #         vector_format=VF_LIST_OF_COMPONENTS))
+
+        for name, field in separate_by_real_and_imag(fields_and_names, real_only):
+            grid.add_pointdata(DataArray(name, field,
+                vector_format=VF_LIST_OF_COMPONENTS))
+
+        from os.path import exists
+        if exists(file_name):
+            raise RuntimeError("output file '%s' already exists"
+                    % file_name)
+
+        with open(file_name, "w") as outf:
+            AppendedDataXMLGenerator(compressor)(grid).write(outf)
+
+
+def make_visualizer(queue, discr, vis_order):
+    from meshmode.discretization.poly_element import \
+            PolynomialWarpAndBlendElementDiscretization
+    vis_discr = PolynomialWarpAndBlendElementDiscretization(
+            discr.cl_context, discr.mesh, order=vis_order,
+            real_dtype=discr.real_dtype)
+    from meshmode.discretization.connection import \
+            make_same_mesh_connection
+    cnx = make_same_mesh_connection(queue, discr, vis_discr)
+
+    return Visualizer(discr, vis_discr, cnx)
+
+# vim: foldmethod=marker

meshmode/mesh/__init__.py

@@ -82,6 +82,10 @@ class MeshElementGroup(Record):
         """
 
         if unit_nodes is None:
+            if dim is None:
+                raise TypeError("'dim' must be passed "
+                        "if 'unit_nodes' is not passed")
+
             unit_nodes = mp.warp_and_blend_nodes(dim, order)
 
         dims = unit_nodes.shape[0]

meshmode/mesh/io.py

+from __future__ import division
+
+__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
+
+from meshpy.gmsh_reader import (  # noqa
+        GmshMeshReceiverBase, FileSource, LiteralSource)
+
+
+__doc__ = """
+
+.. autoclass:: FileSource
+.. autoclass:: LiteralSource
+
+.. autofunction:: read_gmsh
+.. autofunction:: generate_gmsh
+
+"""
+
+
+# {{{ gmsh receiver
+
+class GmshMeshReceiver(GmshMeshReceiverBase):
+    def __init__(self):
+        # Use data fields similar to meshpy.triangle.MeshInfo and
+        # meshpy.tet.MeshInfo
+        self.points = None
+        self.elements = None
+        self.element_types = None
+        self.element_markers = None
+        self.tags = None
+
+    def set_up_nodes(self, count):
+        # Preallocate array of nodes within list; treat None as sentinel value.
+        # Preallocation not done for performance, but to assign values at indices
+        # in random order.
+        self.points = [None] * count
+
+    def add_node(self, node_nr, point):
+        self.points[node_nr] = point
+
+    def finalize_nodes(self):
+        self.points = np.array(self.points, dtype=np.float64)
+
+    def set_up_elements(self, count):
+        # Preallocation of arrays for assignment elements in random order.
+        self.element_vertices = [None] * count
+        self.element_nodes = [None] * count
+        self.element_types = [None] * count
+        self.element_markers = [None] * count
+        self.tags = []
+
+    def add_element(self, element_nr, element_type, vertex_nrs,
+            lexicographic_nodes, tag_numbers):
+        self.element_vertices[element_nr] = vertex_nrs
+        self.element_nodes[element_nr] = lexicographic_nodes
+        self.element_types[element_nr] = element_type
+        self.element_markers[element_nr] = tag_numbers
+
+    def finalize_elements(self):
+        pass
+
+    def add_tag(self, name, index, dimension):
+        pass
+
+    def finalize_tags(self):
+        pass
+
+    def get_mesh(self):
+        el_type_hist = {}
+        for el_type in self.element_types:
+            el_type_hist[el_type] = el_type_hist.get(el_type, 0) + 1
+
+        groups = self.groups = []
+
+        ambient_dim = self.points.shape[-1]
+
+        mesh_bulk_dim = max(
+                el_type.dimensions for el_type in el_type_hist.iterkeys())
+
+        # {{{ build vertex numbering
+
+        vertex_index_gmsh_to_mine = {}
+        for el_vertices, el_type in zip(
+                self.element_vertices, self.element_types):
+            for gmsh_vertex_nr in el_vertices:
+                if gmsh_vertex_nr not in vertex_index_gmsh_to_mine:
+                    vertex_index_gmsh_to_mine[gmsh_vertex_nr] = \
+                            len(vertex_index_gmsh_to_mine)
+
+        # }}}
+
+        # {{{ build vertex array
+
+        gmsh_vertex_indices, my_vertex_indices = \
+                zip(*vertex_index_gmsh_to_mine.iteritems())
+        vertices = np.empty(
+                (ambient_dim, len(vertex_index_gmsh_to_mine)), dtype=np.float64)
+        vertices[:, np.array(my_vertex_indices, np.intp)] = \
+                self.points[np.array(gmsh_vertex_indices, np.intp)].T
+
+        # }}}
+
+        from meshmode.mesh import MeshElementGroup, Mesh
+
+        for group_el_type, ngroup_elements in el_type_hist.iteritems():
+            if group_el_type.dimensions != mesh_bulk_dim:
+                continue
+
+            nodes = np.empty((ambient_dim, ngroup_elements, el_type.node_count()),
+                    np.float64)
+            el_vertex_count = group_el_type.vertex_count()
+            vertex_indices = np.empty(
+                    (ngroup_elements, el_vertex_count),
+                    np.float64)
+            i = 0
+
+            for el_vertices, el_nodes, el_type in zip(
+                    self.element_vertices, self.element_nodes, self.element_types):
+                if el_type is not group_el_type:
+                    continue
+
+                nodes[:, i] = self.points[el_nodes].T
+                vertex_indices[i] = [vertex_index_gmsh_to_mine[v_nr]
+                        for v_nr in el_vertices]
+
+                i += 1
+
+            unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(),
+                    dtype=np.float64).T/group_el_type.order)*2 - 1
+
+            groups.append(MeshElementGroup(
+                group_el_type.order,
+                vertex_indices,
+                nodes,
+                unit_nodes=unit_nodes
+                ))
+
+        return Mesh(vertices, groups)
+
+# }}}
+
+
+def read_gmsh(filename, force_ambient_dimension=None):
+    """Read a gmsh mesh file from *filename* and return a
+    :class:`meshmode.mesh.Mesh`.
+
+    :arg force_dimension: if not None, truncate point coordinates to
+        this many dimensions.
+    """
+    from meshpy.gmsh_reader import read_gmsh
+    recv = GmshMeshReceiver()
+    read_gmsh(recv, filename, force_dimension=force_ambient_dimension)
+    return recv.get_mesh()
+
+
+def generate_gmsh(source, dimensions, order=None, other_options=[],
+        extension="geo", gmsh_executable="gmsh", force_ambient_dimension=None):
+    """Run :cmd:`gmsh` on the input given by *source*, and return a
+    :class:`meshmode.mesh.Mesh` based on the result.
+
+    :arg source: an instance of either :class:`FileSource` or
+        :class:`LiteralSource`
+    :arg force_ambient_dimension: if not None, truncate point coordinates to
+        this many dimensions.
+    """
+    recv = GmshMeshReceiver()
+
+    from meshpy.gmsh import GmshRunner
+    from meshpy.gmsh_reader import parse_gmsh
+    with GmshRunner(source, dimensions, order=order,
+            other_options=other_options, extension=extension,
+            gmsh_executable=gmsh_executable) as runner:
+        parse_gmsh(recv, runner.output_file,
+                force_dimension=force_ambient_dimension)
+
+    return recv.get_mesh()