diff --git a/doc/conf.py b/doc/conf.py index 866653625acb09a79319be6bec7c07872081f38f..8ff4e0bb27f8135e1dc077983c3de932a876cb5c 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -280,5 +280,6 @@ intersphinx_mapping = { 'https://documen.tician.de/pyopencl': None, 'https://documen.tician.de/meshpy': None, 'https://documen.tician.de/modepy': None, - 'https://documen.tician.de/loopy': None + 'https://documen.tician.de/loopy': None, + 'https://firedrakeproject.org/' : None } diff --git a/doc/images/firedrake_mesh_design.gv b/doc/images/firedrake_mesh_design.gv new file mode 100755 index 0000000000000000000000000000000000000000..b546539028893ec86138c327623c95c040a78bc2 --- /dev/null +++ b/doc/images/firedrake_mesh_design.gv @@ -0,0 +1,25 @@ +// created with graphviz2.38 dot + +digraph{ + // NODES + + top [label="Topological\nMesh"]; + ref [label="Reference\nElement"]; + fspace [label="Function Space"]; + coordless [label="Coordinateless\nFunction"]; + geo [label="Geometric\nMesh"]; + withgeo [label="With\nGeometry"]; + + // EDGES + + top -> fspace; + ref -> fspace; + + fspace -> coordless; + + top -> geo; + coordless -> geo [label="Mesh\nCoordinates"]; + + fspace -> withgeo; + geo -> withgeo; +} diff --git a/doc/images/firedrake_mesh_design.png b/doc/images/firedrake_mesh_design.png new file mode 100755 index 0000000000000000000000000000000000000000..f394f9e7f3d56138ff15850311e6bcd6dcdd7bbe Binary files /dev/null and b/doc/images/firedrake_mesh_design.png differ diff --git a/doc/index.rst b/doc/index.rst index 3a7a2ba71dab133f5f754ef2d213c01f23d97d32..0beee875eddae0e9b6d2b4d81ebb2d864bd8798f 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -9,6 +9,7 @@ Contents: mesh discretization connection + interop misc Indices and tables diff --git a/doc/interop.rst b/doc/interop.rst new file mode 100644 index 0000000000000000000000000000000000000000..41ad3b636ab1e6cccd4bf2289a901d1889f9c676 --- /dev/null +++ b/doc/interop.rst @@ -0,0 +1,94 @@ +interop +======= + +Interfacing with data outside of :mod:`meshmode`. + + +fiat +---- + +.. automodule:: meshmode.interop.fiat.simplex_cell + + +FInAT +----- + +.. automodule:: meshmode.interop.FInAT.lagrange_element + + +Firedrake +--------- + +TODO include this + + +Implementation Details +---------------------- + + +Firedrake Function Space Design +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + +In firedrake, meshes and function spaces have a close relationship. +In particular, due to some structure described in this +`firedrake pull request `_. +``fd2mm`` mimics this firedrake design style. + +In short, it is the idea +that every function space should have a mesh, and the coordinates of the mesh +should be representable as a function on that same mesh, which must live +on some function space on the mesh... etc. +Under the hood, we divide between topological and geometric objects, +roughly as so + +(1) A reference element defined using :mod:`FInAT` and :mod:`fiat` + is used to define what meshmode calls the unit nodes and unit + vertices. It is worth noting that :mod:`firedrake` does + not require a positive orientation of elements and that its + reference traingle is different than specified in :mod:`modepy`. + +(2) A ``MeshTopology`` which holds information about connectivity + and other topological properties, but nothing about geometry/coordinates + etc. + +(3) A ``FunctionSpace`` created from a ``FInAT`` element and a + ``MeshTopology`` which allows us to define functions + mapping the nodes (defined by the ``FInAT`` element) of + each element in the ``MeshTopology`` to some values + +(4) A ``CoordinatelessFunction`` (in the sense that its + *domain* has no coordinates) which is a function + in a ``FunctionSpace`` + +(5) A ``MeshGeometry`` created from a ``FunctionSpace`` + and a ``CoordinatelessFunction`` in that ``FunctionSpace`` + which maps each dof to its geometric coordinates. + +(6) A ``WithGeometry`` which is a ``FunctionSpace`` together + with a ``MeshGeometry``. This is the object returned + usually returned to the user by a call + to the :mod:`firedrake` function :func:`FunctionSpace`. + +(7) A ``Function`` is defined on a ``WithGeometry`` + +Thus, by the coordinates of a mesh geometry we mean + +(a) On the hidden back-end: a ``CoordinatelessFunction`` *f* on some function + space defined only on the mesh topology +(b) On the front-end: A ``Function`` with the values of *f* but defined + on a ``WithGeometry`` created from the ``FunctionSpace`` *f* lives in + and the ``MeshGeometry`` *f* defines. + +Basically, it's this picture (where a->b if b depends on a) + +.. warning:: + + In general, the ``FunctionSpace`` of the coordinates function + of a ``WithGeometry`` may not be the same ``FunctionSpace`` + as for functions which live in the ``WithGeometry``. + This picture + only shows how the class definitions depend on each other. + + +.. image:: images/firedrake_mesh_design.png diff --git a/meshmode/interop/FInAT/__init__.py b/meshmode/interop/FInAT/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..a79ad2ae3c12bc9b82afb3d17d15b26c81e906d7 --- /dev/null +++ b/meshmode/interop/FInAT/__init__.py @@ -0,0 +1,25 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 meshmode.interop.FInAT.lagrange_element import FinatLagrangeElementImporter + +__all__ = ['FinatLagrangeElementImporter'] diff --git a/meshmode/interop/FInAT/lagrange_element.py b/meshmode/interop/FInAT/lagrange_element.py new file mode 100644 index 0000000000000000000000000000000000000000..4b3d977c04ba8a0db2517f78290e7b7147d275b0 --- /dev/null +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -0,0 +1,218 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 numpy.linalg as la + +from meshmode.interop import ExternalImportHandler +from meshmode.interop.fiat import FIATSimplexCellImporter + + +__doc__ = """ +.. autoclass:: FinatLagrangeElementImporter + :members: +""" + + +class FinatLagrangeElementImporter(ExternalImportHandler): + """ + An importer for a FInAT element, usually instantiated from + ``some_instantiated_firedrake_function_space.finat_element`` + """ + def __init__(self, finat_element): + """ + :param finat_element: A FInAT element of type + :class:`finat.fiat_elements.Lagrange` or + :class:`finat.fiat_elements.DiscontinuousLagrange` + (or :class:`finat.spectral.GaussLegendre` or + :class:`finat.spectral.GaussLobattoLegendre` in 1D) + which uses affine mapping of the basis functions + (i.e. ``finat_element.mapping`` must be + ``"affine"``) + + :raises TypeError: If :param:`finat_element` is not of type + :class:`finat.fiat_elements.Lagrange` or + :class:`finat.fiat_elements.DiscontinuousLagrange` + (or the 1D classes listed above) + + :raises ValueError: If :param:`finat_element` does not + use affine mappings of the basis functions + """ + # {{{ Check and store input + + # Check types + from finat.fiat_elements import DiscontinuousLagrange, Lagrange + from finat.spectral import GaussLegendre, GaussLobattoLegendre + valid_types = (Lagrange, DiscontinuousLagrange, + GaussLegendre, GaussLobattoLegendre) + if not isinstance(finat_element, valid_types): + raise TypeError(":param:`finat_element` must be of type" + " `finat.fiat_elements.Lagrange`," + " `finat.fiat_elements.DiscontinuousLagrange`," + " or `finat.spectral.GaussLegendre` or" + " `finat.spectral.GaussLobattoLegendre` in 1D," + " but is instead of type `%s`" % type(finat_element)) + + if finat_element.mapping != 'affine': + raise ValueError("FInAT element must use affine mappings" + " of the bases") + # }}} + + super(FinatLagrangeElementImporter, self).__init__(finat_element) + + self.cell_importer = FIATSimplexCellImporter(finat_element.cell) + + # computed and stored once :meth:`unit_nodes`, :meth:`unit_vertices`, + # and :meth:`flip_matrix` are called + self._unit_nodes = None + self._unit_vertex_indices = None + self._flip_matrix = None + + def _compute_unit_vertex_indices_and_nodes(self): + """ + Compute the unit nodes, as well as the unit vertex indices, + if they have not already been computed. + """ + if self._unit_nodes is None or self._unit_vertex_indices is None: + # FIXME : This should work, but uses some private info + # {{{ Compute unit nodes + + # each point evaluator is a function p(f) evaluates f at a node, + # so we need to evaluate each point evaluator at the identity to + # recover the nodes + point_evaluators = self.data._element.dual.nodes + unit_nodes = [p(lambda x: x) for p in point_evaluators] + unit_nodes = np.array(unit_nodes).T + self._unit_nodes = \ + self.cell_importer.affinely_map_firedrake_to_meshmode(unit_nodes) + + # Is this safe?, I think so bc on a reference element + close = 1e-8 + # Get vertices as (dim, nunit_vertices) + unit_vertices = np.array(self.data.cell.vertices).T + unit_vertices = \ + self.cell_importer.affinely_map_firedrake_to_meshmode(unit_vertices) + self._unit_vertex_indices = [] + for n_ndx in range(self._unit_nodes.shape[1]): + for v_ndx in range(unit_vertices.shape[1]): + diff = self._unit_nodes[:, n_ndx] - unit_vertices[:, v_ndx] + if np.max(np.abs(diff)) < close: + self._unit_vertex_indices.append(n_ndx) + break + + self._unit_vertex_indices = np.array(self._unit_vertex_indices) + + # }}} + + def dim(self): + """ + :return: The dimension of the FInAT element's cell + """ + return self.cell_importer.data.get_dimension() + + def unit_vertex_indices(self): + """ + :return: A numpy integer array of indices + so that *self.unit_nodes()[self.unit_vertex_indices()]* + are the nodes of the reference element which coincide + with its vertices (this is possibly empty). + """ + self._compute_unit_vertex_indices_and_nodes() + return self._unit_vertex_indices + + def unit_nodes(self): + """ + :return: The unit nodes used by the FInAT element mapped + onto the appropriate :mod:`modepy` `reference + element `_ + as an array of shape *(dim, nunit_nodes)*. + """ + self._compute_unit_vertex_indices_and_nodes() + return self._unit_nodes + + def nunit_nodes(self): + """ + :return: The number of unit nodes. + """ + return self.unit_nodes().shape[1] + + def flip_matrix(self): + """ + :return: The matrix which should be applied to the + *(dim, nunitnodes)*-shaped array of nodes corresponding to + an element in order to change orientation - <-> +. + + The matrix will be *(dim, dim)* and orthogonal with + *np.float64* type entries. + """ + if self._flip_matrix is None: + # This is very similar to :mod:`meshmode` in processing.py + # the function :function:`from_simplex_element_group`, but + # we needed to use firedrake nodes + + from modepy.tools import barycentric_to_unit, unit_to_barycentric + + # Generate a resampling matrix that corresponds to the + # first two barycentric coordinates being swapped. + + bary_unit_nodes = unit_to_barycentric(self.unit_nodes()) + + flipped_bary_unit_nodes = bary_unit_nodes.copy() + flipped_bary_unit_nodes[0, :] = bary_unit_nodes[1, :] + flipped_bary_unit_nodes[1, :] = bary_unit_nodes[0, :] + flipped_unit_nodes = barycentric_to_unit(flipped_bary_unit_nodes) + + from modepy import resampling_matrix, simplex_best_available_basis + + flip_matrix = resampling_matrix( + simplex_best_available_basis(self.dim(), self.data.degree), + flipped_unit_nodes, self.unit_nodes()) + + flip_matrix[np.abs(flip_matrix) < 1e-15] = 0 + + # Flipping twice should be the identity + assert la.norm( + np.dot(flip_matrix, flip_matrix) + - np.eye(len(flip_matrix))) < 1e-13 + + self._flip_matrix = flip_matrix + + return self._flip_matrix + + def make_resampling_matrix(self, element_grp): + """ + :param element_grp: A + :class:`meshmode.discretization.InterpolatoryElementGroupBase` whose + basis functions span the same space as the FInAT element. + :return: A matrix which resamples a function sampled at + the firedrake unit nodes to a function sampled at + *element_grp.unit_nodes()* (by matrix multiplication) + """ + from meshmode.discretization import InterpolatoryElementGroupBase + assert isinstance(element_grp, InterpolatoryElementGroupBase), \ + "element group must be an interpolatory element group so that" \ + " can redistribute onto its nodes" + + from modepy import resampling_matrix + return resampling_matrix(element_grp.basis(), + new_nodes=element_grp.unit_nodes, + old_nodes=self.unit_nodes()) diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..2c1bf0727a0becb3a13d9d6dbad12b751d60db39 --- /dev/null +++ b/meshmode/interop/__init__.py @@ -0,0 +1,72 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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. +""" + +__doc__ = """ +Development Interface +--------------------- +.. autoclass:: ExternalDataHandler +.. autoclass:: ExternalExportHandler +.. autoclass:: ExternalImportHandler +""" + + +# {{{ Generic, most abstract class for transporting meshmode <-> external + +class ExternalDataHandler: + """ + A data handler takes data from meshmode and facilitates its use + in another package or the reverse: takes data from another package + and facilitates its use in meshmode. + + .. attribute:: data + + The object which needs to be interfaced either into meshmode or + out of meshmode. + Should not be modified after creation. + """ + def __init__(self, data): + self.data = data + +# }}} + + +# {{{ Define specific classes for meshmode -> external and meshmode <- external + +class ExternalExportHandler(ExternalDataHandler): + """ + Subclass of :class:`ExternalDataHandler` for meshmode -> external + data transfer + """ + pass + + +class ExternalImportHandler(ExternalDataHandler): + """ + Subclass of :class:`ExternalDataHandler` for external -> meshmode + data transfer + """ + pass + +# }}} + + +# vim: fdm=marker diff --git a/meshmode/interop/fiat/__init__.py b/meshmode/interop/fiat/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..57a42bfe47fb1ed78cb8be603326b222761add2f --- /dev/null +++ b/meshmode/interop/fiat/__init__.py @@ -0,0 +1,25 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 meshmode.interop.fiat.simplex_cell import FIATSimplexCellImporter + +__all__ = ['FIATSimplexCellImporter'] diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py new file mode 100644 index 0000000000000000000000000000000000000000..c6c0ca2ca4e4241a0baa99860a083214645b4de5 --- /dev/null +++ b/meshmode/interop/fiat/simplex_cell.py @@ -0,0 +1,122 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 numpy.linalg as la + +from meshmode.interop import ExternalImportHandler + + +__doc__ = """ +.. autoclass:: FIATSimplexCellImporter +""" + + +# {{{ Compute an affine mapping from given input/outputs + +def get_affine_mapping(reference_vects, vects): + r""" + Returns (mat, shift), + a matrix *mat* and vector *shift* which maps the + *i*th vector in *reference_vects* to the *i*th vector in *vects by + + ..math:: + + A ri + b -> vi, \qquad A = mat, b = shift + + :arg reference_vects: An np.array of *n* vectors of dimension *ref_dim* + :arg vects: An np.array of *n* vectors of dimension *dim*, with + *ref_dim* <= *dim*. + + NOTE : Should be shape *(ref_dim, nvectors)*, *(dim, nvectors)* respectively. + + *mat* will have shape *(dim, ref_dim)*, *shift* will have shape *(dim,)* + """ + # Make sure both have same number of vectors + ref_dim, num_vects = reference_vects.shape + assert num_vects == vects.shape[1] + + # Make sure d1 <= d2 (see docstring) + dim = vects.shape[0] + assert ref_dim <= dim + + # If there is only one vector, set M = I, b = vect - reference + if num_vects == 1: + mat = np.eye(dim, ref_dim) + shift = vects[:, 0] - np.matmul(mat, reference_vects[:, 0]) + else: + ref_span_vects = reference_vects[:, 1:] - reference_vects[:, 0, np.newaxis] + span_vects = vects[:, 1:] - vects[:, 0, np.newaxis] + mat = la.solve(ref_span_vects, span_vects) + shift = -np.matmul(mat, reference_vects[:, 0]) + vects[:, 0] + + return mat, shift + +# }}} + + +# {{{ Interoperator for FIAT's reference_element.Simplex + +class FIATSimplexCellImporter(ExternalImportHandler): + """ + Import handler for a :mod:`FIAT` simplex cell. + + .. attribute:: data + + An instance of :class:`fiat.FIAT.reference_element.Simplex`. + """ + def __init__(self, cell): + """ + :arg cell: a :class:`fiat.FIAT.reference_element.Simplex`. + """ + # Ensure this cell is actually a simplex + from FIAT.reference_element import Simplex + assert isinstance(cell, Simplex) + + super(FIATSimplexCellImporter, self).__init__(cell) + + # Stored as (dim, nunit_vertices) + from modepy.tools import unit_vertices + self._unit_vertices = unit_vertices(cell.get_dimension()).T + + # Maps FIAT reference vertices to :mod:`meshmode` + # unit vertices by x -> Ax + b, where A is :attr:`_mat` + # and b is :attr:`_shift` + reference_vertices = np.array(cell.vertices).T + self._mat, self._shift = get_affine_mapping(reference_vertices, + self._unit_vertices) + + def affinely_map_firedrake_to_meshmode(self, points): + """ + Map points on the firedrake reference simplex to + :mod:`modepy` + `unit coordinates `_. + + :arg points: *n* points on the reference simplex + as a numpy array of shape *(dim, n)* + :return: A numpy array of shape *(dim, n)* wherein the + firedrake refernece simplex has been affinely mapped + to the modepy reference simplex + """ + return np.matmul(self._mat, points) + self._shift[:, np.newaxis] + +# }}} diff --git a/meshmode/interop/firedrake/__init__.py b/meshmode/interop/firedrake/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..33daeb392b3276f855867ff79daaaea0ab3e2f64 --- /dev/null +++ b/meshmode/interop/firedrake/__init__.py @@ -0,0 +1,255 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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. +""" + +__doc__ = """ +.. autoclass:: FromFiredrakeConnection +.. autofunction:: import_firedrake_mesh +.. autofunction:: import_firedrake_function_space +""" + +import numpy as np + + +# {{{ Helper functions to construct importers for firedrake objects + + +def _compute_cells_near_bdy(mesh, bdy_id): + """ + Returns an array of the cell ids with >= 1 vertex on the + given bdy_id + """ + cfspace = mesh.coordinates.function_space() + cell_node_list = cfspace.cell_node_list + + boundary_nodes = cfspace.boundary_nodes(bdy_id, 'topological') + # Reduce along each cell: Is a vertex of the cell in boundary nodes? + cell_is_near_bdy = np.any(np.isin(cell_node_list, boundary_nodes), axis=1) + + return np.arange(cell_node_list.shape[0], dtype=np.int32)[cell_is_near_bdy] + + +def import_firedrake_mesh(fdrake_mesh, near_bdy=None): + """ + :param fdrake_mesh: an instance of :mod:`firedrake`'s mesh, i.e. + of type :class:`firedrake.mesh.MeshGeometry`. + + :param near_bdy: If *None* does nothing. Otherwise should be a + boundary id. In this case, only cell ids + with >= 1 vertex on the given bdy id are used. + + :returns: A :class:`FiredrakeMeshGeometryImporter` object which + is created appropriately from :param:`fdrake_mesh` + + :raises TypeError: if :param:`fdrake_mesh` is of the wrong type + """ + # {{{ Make sure have an initialized firedrake mesh + + from firedrake.mesh import MeshGeometry, MeshTopology + if not isinstance(fdrake_mesh, MeshGeometry): + if isinstance(fdrake_mesh, MeshTopology): + raise TypeError(":param:`fdrake_mesh` must have an associated" + " geometry, but is of type MeshTopology") + raise TypeError(":param:`fdrake_mesh` must be of type" + "`firedrake.mesh.MeshGeometry`, not %s " + % type(fdrake_mesh)) + + fdrake_mesh.init() + + # }}} + + # {{{ Use the coordinates function space to built an importer + + coords_fspace = fdrake_mesh.coordinates.function_space() + cells_to_use = None + if near_bdy is not None: + cells_to_use = _compute_cells_near_bdy(fdrake_mesh, near_bdy) + + from meshmode.interop.FInAT.lagrange_element import \ + FinatLagrangeElementImporter + from meshmode.interop.firedrake.mesh_topology import \ + FiredrakeMeshTopologyImporter + topology_importer = FiredrakeMeshTopologyImporter(fdrake_mesh, + cells_to_use=cells_to_use) + finat_elt_importer = FinatLagrangeElementImporter(coords_fspace.finat_element) + + from meshmode.interop.firedrake.function_space_coordless import \ + FiredrakeFunctionSpaceImporter, FiredrakeCoordinatelessFunctionImporter + + coords_fspace_importer = FiredrakeFunctionSpaceImporter(coords_fspace, + topology_importer, + finat_elt_importer) + coordinates_importer = \ + FiredrakeCoordinatelessFunctionImporter(fdrake_mesh.coordinates, + coords_fspace_importer) + + # }}} + + # FIXME Allow for normals and no_normals_warn? as in the constructor + # for FiredrakeMeshGeometryImporter s + from meshmode.interop.firedrake.mesh_geometry import \ + FiredrakeMeshGeometryImporter + return FiredrakeMeshGeometryImporter(fdrake_mesh, coordinates_importer) + + +def import_firedrake_function_space(cl_ctx, fdrake_fspace, mesh_importer): + """ + Builds a + :class:`FiredrakeWithGeometryImporter` built from the given + :mod:`firedrake` function space :param:`fdrake_fspace` + + :param cl_ctx: A pyopencl computing context. This input + is not checked. + + :param:`fdrake_fspace` An instance of class :mod:`firedrake`'s + :class:`WithGeometry` class, representing a function + space defined on a concrete mesh + + :param mesh_importer: An instance of class + :class:`FiredrakeMeshGeometryImporter` defined on the + mesh underlying :param:`fdrake_fspace`. + + :returns: An instance of class + :mod :`meshmode.interop.firedrake.function_space` + :class:`FiredrakeWithGeometryImporter` built from the given + :param:`fdrake_fspace` + + :raises TypeError: If any input is the wrong type + :raises ValueError: If :param:`mesh_importer` is built on a mesh + different from the one underlying :param:`fdrake_fspace`. + """ + # {{{ Input checking + + from firedrake.functionspaceimpl import WithGeometry + if not isinstance(fdrake_fspace, WithGeometry): + raise TypeError(":param:`fdrake_fspace` must be of type " + ":class:`firedrake.functionspaceimpl.WithGeometry`, " + "not %s " % type(fdrake_fspace)) + + from meshmode.interop.firedrake.mesh_geometry import \ + FiredrakeMeshGeometryImporter + if not isinstance(mesh_importer, FiredrakeMeshGeometryImporter): + raise TypeError(":param:`mesh_importer` must be of type " + ":class:`FiredrakeMeshGeometryImporter` not `%s`." + % type(mesh_importer)) + + if mesh_importer.data != fdrake_fspace.mesh(): + raise ValueError("``mesh_importer.data`` and ``fdrake_fspace.mesh()`` " + "must be identical") + + # }}} + + # {{{ Make an importer for the topological function space + + from meshmode.interop.FInAT import FinatLagrangeElementImporter + from meshmode.interop.firedrake.function_space_coordless import \ + FiredrakeFunctionSpaceImporter + mesh_importer.init(cl_ctx) + finat_elt_importer = \ + FinatLagrangeElementImporter(fdrake_fspace.finat_element) + topological_importer = FiredrakeFunctionSpaceImporter(fdrake_fspace, + mesh_importer, + finat_elt_importer) + + # }}} + + # now we can make the full importer + from meshmode.interop.firedrake.function_space import \ + FiredrakeWithGeometryImporter + return FiredrakeWithGeometryImporter(cl_ctx, + fdrake_fspace, + topological_importer, + mesh_importer) + +# }}} + + +class FromFiredrakeConnection: + """ + Creates a method of transporting functions on a Firedrake mesh + back and forth between firedrake and meshmode + + .. attribute:: with_geometry_importer + + An instance of class :class:`FiredrakeWithGeometryImporter` + that converts the firedrake function space into + meshmode and allows for field conversion. + """ + def __init__(self, cl_ctx, fdrake_function_space): + """ + :param cl_ctx: A computing context + :param fdrake_function_space: A firedrake function space (an instance of + class :class:`WithGeometry`) + """ + mesh_importer = import_firedrake_mesh(fdrake_function_space.mesh()) + self.with_geometry_importer = \ + import_firedrake_function_space(cl_ctx, + fdrake_function_space, + mesh_importer) + + def from_function_space(self): + """ + :returns: the firedrake function space this object was created from + """ + return self.with_geometry_importer.data + + def to_factory(self): + """ + :returns: An InterpolatoryQuadratureSimplexGroupFactory + of the appropriate degree to use + """ + return self.with_geometry_importer.factory() + + def to_discr(self): + """ + :returns: A meshmode discretization corresponding to the + firedrake function space + """ + return self.with_geometry_importer.discretization() + + def from_firedrake(self, fdrake_function): + """ + Convert a firedrake function on this function space to a numpy + array + + :param queue: A CommandQueue + :param fdrake_function: A firedrake function on the function space + of this connection + + :returns: A numpy array holding the data of this function as + a field for the corresponding meshmode mesh + """ + return self.with_geometry_importer.convert_function(fdrake_function) + + def from_meshmode(self, field, fdrake_function): + """ + Store a numpy array holding a field on the meshmode mesh into a + firedrake function + + :param field: A numpy array representing a field on the meshmode version + of this mesh + :param fdrake_function: The firedrake function to store the data in + """ + from meshmode.interop.firedrake.function import \ + FiredrakeFunctionImporter + importer = FiredrakeFunctionImporter(fdrake_function, + self.with_geometry_importer) + importer.set_from_field(field) diff --git a/meshmode/interop/firedrake/function.py b/meshmode/interop/firedrake/function.py new file mode 100644 index 0000000000000000000000000000000000000000..14ef3a4638956bde35798c32ead4d3a2a689dd3c --- /dev/null +++ b/meshmode/interop/firedrake/function.py @@ -0,0 +1,118 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 firedrake import Function + +from meshmode.interop import ExternalImportHandler +from meshmode.interop.firedrake.function_space import FiredrakeWithGeometryImporter +from meshmode.interop.firedrake.function_space_coordless import \ + FiredrakeCoordinatelessFunctionImporter + + +class FiredrakeFunctionImporter(ExternalImportHandler): + """ + An import handler for :mod:`firedrake` functions. + + .. attribute:: data + + A :mod:`firedrake` function + """ + def __init__(self, function, function_space_importer): + """ + :param function: A firedrake :class:`Function` or a + :class:`FiredrakeCoordinatelessFunctionImporter` + + :param function_space_importer: An instance of + :class:`FiredrakeWithGeometryImporter` + compatible with the function + space this function is built on + """ + # {{{ Check types + if not isinstance(function_space_importer, + FiredrakeWithGeometryImporter): + raise TypeError(":param:`function_space_importer` must be of " + "type FiredrakeWithGeometryImporter") + + if not isinstance(function, (Function, + FiredrakeCoordinatelessFunctionImporter)): + raise TypeError(":param:`function` must be one of " + "(:class:`firedrake.function.Function`," + " FiredrakeCoordinatelessFunctionImporter)") + # }}} + + super(FiredrakeFunctionImporter, self).__init__(function) + + self._function_space_importer = function_space_importer + if isinstance(function, Function): + self._topological_importer = \ + FiredrakeCoordinatelessFunctionImporter(function, + function_space_importer) + else: + self._topological_importer = function + + def function_space_importer(self): + """ + :returns: the :class:`FiredrakeWithGeometryImporter` instance + used for the function space this object's :attr:`data` lives on + """ + return self._function_space_importer + + @property + def topological_importer(self): + """ + :returns: The topological version of this object, i.e. the underlying + coordinateless importer (an instance of + :class:`FiredrakeCoordinatelessFunctionImporter`) + """ + return self._topological_importer + + def as_field(self): + """ + :returns: A numpy array holding the data of this function as a + field on the corresponding meshmode mesh + """ + return self.function_space_importer().convert_function(self) + + def set_from_field(self, field): + """ + Set function :attr:`data`'s data + from a field on the corresponding meshmode mesh. + + :param field: A numpy array representing a field on the corresponding + meshmode mesh + """ + # Handle 1-D case + if len(self.data.dat.data.shape) == 1 and len(field.shape) > 1: + field = field.reshape(field.shape[1]) + + # resample from nodes + group = self.function_space_importer().discretization().groups[0] + resampled = np.copy(field) + resampled_view = group.view(resampled) + resampling_mat = self.function_space_importer().resampling_mat(False) + np.matmul(resampled_view, resampling_mat.T, out=resampled_view) + + # reorder data + self.data.dat.data[:] = self.function_space_importer().reorder_nodes( + resampled, firedrake_to_meshmode=False)[:] diff --git a/meshmode/interop/firedrake/function_space.py b/meshmode/interop/firedrake/function_space.py new file mode 100644 index 0000000000000000000000000000000000000000..85b3d6c22645f57ac2a14dff0e953307222b659f --- /dev/null +++ b/meshmode/interop/firedrake/function_space.py @@ -0,0 +1,206 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 warnings import warn +import numpy as np + +from meshmode.interop import ExternalImportHandler +from meshmode.interop.firedrake.mesh_geometry import \ + FiredrakeMeshGeometryImporter +from meshmode.interop.firedrake.function_space_coordless import \ + FiredrakeFunctionSpaceImporter + +from meshmode.interop.firedrake.function_space_shared_data import \ + FiredrakeFunctionSpaceDataImporter + + +class FiredrakeWithGeometryImporter(ExternalImportHandler): + def __init__(self, + cl_ctx, + function_space, + function_space_importer, + mesh_importer): + """ + :param cl_ctx: A pyopencl context + + :param function_space: A firedrake :class:`WithGeometry` + + :param function_space_importer: An instance of class + :class:`FiredrakeFunctionSpaceImporter` for the + underlying topological function space of :param`function_space` + + :param mesh_importer: An instance of class + :class:`FiredrakeMeshGeometryImporter` for the mesh + that :param:`function_space` is built on + """ + # FIXME use on bdy + # {{{ Check input + from firedrake.functionspaceimpl import WithGeometry + if not isinstance(function_space, WithGeometry): + raise TypeError(":arg:`function_space` must be of type" + " :class:`firedrake.functionspaceimpl.WithGeometry") + + if not isinstance(function_space_importer, + FiredrakeFunctionSpaceImporter): + raise TypeError(":arg:`function_space_importer` must be of type" + " FiredrakeFunctionSpaceImporter") + + if not isinstance(mesh_importer, FiredrakeMeshGeometryImporter): + raise TypeError(":arg:`mesh_importer` must be of type" + " FiredrakeMeshGeometryImporter") + + # Make sure importers are good for given function space + assert function_space_importer.data == function_space.topological + assert mesh_importer.data == function_space.mesh() + + # }}} + + # Initialize as importer + super(FiredrakeWithGeometryImporter, self).__init__(function_space) + self._topology_importer = function_space_importer + self._mesh_importer = mesh_importer + self._cl_ctx = cl_ctx + + finat_element_importer = function_space_importer.finat_element_importer + self._shared_data = \ + FiredrakeFunctionSpaceDataImporter(cl_ctx, + mesh_importer, + finat_element_importer) + + mesh_order = mesh_importer.data.coordinates.\ + function_space().finat_element.degree + if mesh_order > self.data.finat_element.degree: + warn("Careful! When the mesh order is higher than the element" + " order. Conversion MIGHT work, but maybe not..." + " To be honest I really don't know.") + + # Used to convert between refernce node sets + self._resampling_mat_fd2mm = None + self._resampling_mat_mm2fd = None + + def __getattr__(self, attr): + return getattr(self._topology_importer, attr) + + def mesh_importer(self): + return self._mesh_importer + + def _reordering_array(self, firedrake_to_meshmode): + if firedrake_to_meshmode: + return self._shared_data.firedrake_to_meshmode() + return self._shared_data.meshmode_to_firedrake() + + def factory(self): + return self._shared_data.factory() + + def discretization(self): + return self._shared_data.discretization() + + def resampling_mat(self, firedrake_to_meshmode): + if self._resampling_mat_fd2mm is None: + element_grp = self.discretization().groups[0] + self._resampling_mat_fd2mm = \ + self.finat_element_importer.make_resampling_matrix(element_grp) + + self._resampling_mat_mm2fd = np.linalg.inv(self._resampling_mat_fd2mm) + + # return the correct resampling matrix + if firedrake_to_meshmode: + return self._resampling_mat_fd2mm + return self._resampling_mat_mm2fd + + def reorder_nodes(self, nodes, firedrake_to_meshmode=True): + """ + :arg nodes: An array representing function values at each of the + dofs, if :arg:`firedrake_to_meshmode` is *True*, should + be of shape (ndofs) or (ndofs, xtra_dims). + If *False*, should be of shape (ndofs) or (xtra_dims, ndofs) + :arg firedrake_to_meshmode: *True* iff firedrake->meshmode, *False* + if reordering meshmode->firedrake + """ + # {{{ Case where shape is (ndofs,), just apply reordering + + if len(nodes.shape) == 1: + return nodes[self._reordering_array(firedrake_to_meshmode)] + + # }}} + + # {{{ Else we have (xtra_dims, ndofs) or (ndofs, xtra_dims): + + # Make sure we have (xtra_dims, ndofs) ordering + if firedrake_to_meshmode: + nodes = nodes.T + + reordered_nodes = nodes[:, self._reordering_array(firedrake_to_meshmode)] + + # if converting mm->fd, change to (ndofs, xtra_dims) + if not firedrake_to_meshmode: + reordered_nodes = reordered_nodes.T + + return reordered_nodes + + # }}} + + def convert_function(self, function): + """ + Convert a firedrake function defined on this function space + to a meshmode form, reordering data as necessary + and resampling to the unit nodes in meshmode + + :param function: A firedrake function or a + :class:`FiredrakeFunctionImporter` instance + + :returns: A numpy array + """ + from meshmode.interop.firedrake.function import FiredrakeFunctionImporter + if isinstance(function, FiredrakeFunctionImporter): + function = function.data + + # FIXME: Check that function can be converted! Maybe check the + # shared data? We can just do the check below + # but it's a little too stiff + #assert function.function_space() == self.data + + nodes = function.dat.data + + # handle vector function spaces differently, hence the shape checks + + # {{{ Reorder the nodes to have positive orientation + # (and if a vector, now have meshmode [dims][nnodes] + # instead of firedrake [nnodes][dims] shape) + + if len(nodes.shape) > 1: + new_nodes = [self.reorder_nodes(nodes.T[i], True) for i in + range(nodes.shape[1])] + nodes = np.array(new_nodes) + else: + nodes = self.reorder_nodes(nodes, True) + + # }}} + + # {{{ Now convert to meshmode reference nodes + node_view = self.discretization().groups[0].view(nodes) + # Multiply each row (repping an element) by the resampler + np.matmul(node_view, self.resampling_mat(True).T, out=node_view) + + # }}} + + return nodes diff --git a/meshmode/interop/firedrake/function_space_coordless.py b/meshmode/interop/firedrake/function_space_coordless.py new file mode 100644 index 0000000000000000000000000000000000000000..e6fdec6cd2137c0831340bb13db206da3e21f9e5 --- /dev/null +++ b/meshmode/interop/firedrake/function_space_coordless.py @@ -0,0 +1,221 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 warnings import warn # noqa + +from meshmode.interop import ExternalImportHandler +from meshmode.interop.firedrake.mesh_topology import \ + FiredrakeMeshTopologyImporter +from meshmode.interop.FInAT import FinatLagrangeElementImporter + + +# {{{ Function space for coordinateless functions to live on + + +class FiredrakeFunctionSpaceImporter(ExternalImportHandler): + """ + This is a counterpart of :class:`firedrake.functionspaceimpl.FunctionSpace`, + + This is not what usually results from a call to + :func:`firedrake.functionspace.FunctionSpace`. + When someone calls :func:`firedrake.functionspace.FunctionSpace` + from the user side, they are usually talking about a + :class:`firedrake.functionspaceimpl.WithGeometry`. + + Generally, this class is here to match Firedrake's design + principles, i.e. so that we have something to put CoordinatelessFunction + import handlers on. + + Just like we have "topological" and "geometric" + meshes, one can think of think of this as a "topological" function space + whose counterpart is a "WithGeometry" + + .. attribute:: data + + An instance of :class:`firedrake.functionspaceimpl.FunctionSpace`, i.e. + a function space built to hold coordinateless function (functions + whose domain is just a set of points with some topology, but no + coordinates) + + .. attribute:: finat_element_importer + + An instance of + :class:`meshmode.interop.FInAT.FinatLagrangeElementImporter` which + is an importer for the *finat_element* of :attr:`data`. + """ + def __init__(self, function_space, mesh_importer, finat_element_importer): + """ + :param function_space: A :mod:`firedrake` + :class:`firedrake.functionspaceimpl.FunctionSpace` or + :class:`firedrake.functionspaceimpl.WithGeometry`. In the + latter case, the underlying ``FunctionSpace`` is extracted + from the ``WithGeometry``. + :param mesh_importer: An instance + of :class:`FiredrakeMeshTopology` created from the topological + mesh of :param:`function_space` + :param finat_element_importer: An instance + of :class:`FinatLagrangeElementIMporter` created from the + finat_element of :param:`function_space` + + :raises TypeError: If any of the arguments are the wrong type + :raises ValueError: If :param:`mesh_importer` or + :param:`finat_element_importer` are importing + a different mesh or finat element than the provided + :param:`function_space` is built on. + """ + # We want to ignore any geomery + function_space = function_space.topological + mesh_importer = mesh_importer.topological_importer + + # {{{ Some type-checking + from firedrake.functionspaceimpl import FunctionSpace, WithGeometry + + if not isinstance(function_space, (FunctionSpace, WithGeometry)): + raise TypeError(":param:`function_space` must be of type " + "``firedrake.functionspaceimpl.FunctionSpace`` " + "or ``firedrake.functionspaceimpl.WithGeometry`` ", + "not %s" % type(function_space)) + + if not isinstance(mesh_importer, FiredrakeMeshTopologyImporter): + raise TypeError(":param:`mesh_importer` must be either *None* " + "or of type :class:`meshmode.interop.firedrake." + "FiredrakeMeshTopologyImporter`") + if not function_space.mesh().topological == mesh_importer.data: + raise ValueError(":param:`mesh_importer`'s *data* attribute " + "must be the same mesh as returned by " + ":param:`function_space`'s *mesh()* method.") + + if not isinstance(finat_element_importer, FinatLagrangeElementImporter): + raise TypeError(":param:`finat_element_importer` must be either " + "*None* " + "or of type :class:`meshmode.interop.FInAT." + "FinatLagragneElementImporter`") + if not function_space.finat_element == finat_element_importer.data: + raise ValueError(":param:`finat_element_importer`'s *data* " + "attribute " + "must be the same finat element as " + ":param:`function_space`'s *finat_element*" + " attribute.") + + # }}} + + # finish initialization + super(FiredrakeFunctionSpaceImporter, self).__init__(function_space) + + self._mesh_importer = mesh_importer + self.finat_element_importer = finat_element_importer + + @property + def topological_importer(self): + """ + A reference to self for compatability with 'geometrical' function spaces + """ + return self + + def mesh_importer(self): + """ + Return this object's mesh importer + """ + return self._mesh_importer + +# }}} + + +# {{{ Container to hold the coordinateless functions themselves + + +class FiredrakeCoordinatelessFunctionImporter(ExternalImportHandler): + """ + A coordinateless function, i.e. a function defined on a set of + points which only have an associated topology, no associated + geometry. + + .. attribute:: data + + An instance of :mod:`firedrake` class + :class:`firedrake.function.CoordinatelessFunction`. + Note that a coordinateless function object in firedrake + references concrete, but mutable data. + """ + def __init__(self, function, function_space_importer): + """ + :param function: The instance of + :class:`firedrake.function.CoordinatelessFunction` + which this object is importing. Becomes the + :attr:`data` attribute. + + :param function_space_importer: An instance of + :class:`FiredrakeFunctionSpaceImporter` + which is importing the topological function space that + :param:`function` is built on. + + :raises TypeError: If either parameter is the wrong type + :raises ValueError: If :param:`function_space_importer` is an + importer for a firedrake function space which is not + identical to ``function.topological.function_space()`` + """ + # Throw out geometric information if there is any + function = function.topological + function_space_importer = function_space_importer.topological_importer + + # {{{ Some type-checking + + from firedrake.function import Function, CoordinatelessFunction + if not isinstance(function, (CoordinatelessFunction, Function)): + raise TypeError(":param:`function` must be of type " + "`firedrake.function.CoordinatelessFunction` " + " or `firedrdake.function.Function`") + + if not isinstance(function_space_importer, + FiredrakeFunctionSpaceImporter): + raise TypeError(":param:`function_space_importer` must be of type " + "`meshmode.interop.firedrake." + "FiredrakeFunctionSpaceImporter`.") + + if not function_space_importer.data == function.function_space(): + raise ValueError(":param:`function_space_importer`'s *data* " + "attribute and ``function.function_space()`` " + "must be identical.") + + super(FiredrakeCoordinatelessFunctionImporter, self).__init__(function) + + self._function_space_importer = function_space_importer + + def function_space_importer(self): + """ + Return the + :class:`meshmode.interop.firedrake.FiredrakeFunctionSpaceImporter` + instance being used to import the function space + of the underlying firedrake function this object represents. + """ + return self._function_space_importer + + @property + def topological_importer(self): + """ + A reference to self for compatability with functions that have + coordinates. + """ + return self + + +# }}} diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py new file mode 100644 index 0000000000000000000000000000000000000000..c7b316a7b3da443c8647885ac0badee0eef945ad --- /dev/null +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -0,0 +1,265 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 numpy.linalg as la +from decorator import decorator +from firedrake.functionspacedata import FunctionSpaceData +from meshmode.discretization import Discretization + +from meshmode.interop import ExternalImportHandler +from meshmode.interop.FInAT.lagrange_element import \ + FinatLagrangeElementImporter + + +@decorator +def cached(f, mesh_importer, key, *args, **kwargs): + """ + Exactly :func:`firedrake.functionspacedata.cached`, but + caching on mesh Geometry instead + + :param f: The function to cache. + :param mesh_importer: The mesh_importer to cache on (should have a + ``_shared_data_cache`` object). + :param key: The key to the cache. + :args args: Additional arguments to ``f``. + :kwargs kwargs: Additional keyword arguments to ``f``.""" + assert hasattr(mesh_importer, "_shared_data_cache") + cache = mesh_importer._shared_data_cache[f.__name__] + try: + return cache[key] + except KeyError: + result = f(mesh_importer, key, *args, **kwargs) + cache[key] = result + return result + + +def reorder_nodes(orient, nodes, flip_matrix, unflip=False): + """ + flips :param:`nodes` according to :param:`orient` + + :param orient: An array of shape *(nelements)* of orientations, + >0 for positive, <0 for negative + :param nodes: a *(nelements, nunit_nodes)* or + *(dim, nelements, nunit_nodes)* shaped array of nodes + :param flip_matrix: The matrix used to flip each negatively-oriented + element + :param unflip: If *True*, use transpose of :param:`flip_matrix` to + flip negatively-oriented elements + """ + # reorder nodes (Code adapted from + # meshmode.mesh.processing.flip_simplex_element_group) + + # ( round to int bc applying on integers) + flip_mat = np.rint(flip_matrix) + if unflip: + flip_mat = flip_mat.T + + # flipping twice should be identity + assert la.norm( + np.dot(flip_mat, flip_mat) + - np.eye(len(flip_mat))) < 1e-13 + + # }}} + + # {{{ flip nodes that need to be flipped, note that this point we act + # like we are in a DG space + + # if a vector function space, nodes array is shaped differently + if len(nodes.shape) > 2: + nodes[orient < 0] = np.einsum( + "ij,ejk->eik", + flip_mat, nodes[orient < 0]) + # Reshape to [nodes][vector dims] + nodes = nodes.reshape( + nodes.shape[0] * nodes.shape[1], nodes.shape[2]) + # pytential wants [vector dims][nodes] not [nodes][vector dims] + nodes = nodes.T.copy() + else: + nodes[orient < 0] = np.einsum( + "ij,ej->ei", + flip_mat, nodes[orient < 0]) + # convert from [element][unit_nodes] to + # global node number + nodes = nodes.flatten() + + +@cached +def reordering_array(mesh_importer, key, fspace_data): + """ + :param key: A tuple (finat_element_importer, firedrake_to_meshmode) + where *firedrake_to_meshmode* is a *bool*, *True* indicating + firedrake->meshmode reordering, *False* meshmode->firedrake + + :param fspace_data: A :mod:`firedrake` instance of shared + function space data, i.e. + :class:`firedrake.functionspacedata.FunctionSpaceData` + + :returns: an *np.array* that can reorder the data by composition, + see :meth:`reorder_nodes` in class :class:`FiredrakeWithGeometryImporter` + """ + finat_element_importer, firedrake_to_meshmode = key + assert isinstance(finat_element_importer, FinatLagrangeElementImporter) + + cell_node_list = fspace_data.entity_node_lists[mesh_importer.data.cell_set] + if mesh_importer.icell_to_fd is not None: + cell_node_list = cell_node_list[mesh_importer.icell_to_fd] + + num_fd_nodes = fspace_data.node_set.size + + nelements = cell_node_list.shape[0] + nunit_nodes = cell_node_list.shape[1] + num_mm_nodes = nelements * nunit_nodes + + if firedrake_to_meshmode: + nnodes = num_fd_nodes + else: + nnodes = num_mm_nodes + order = np.arange(nnodes) + + # Put into cell-node list if firedrake-to meshmode (so can apply + # flip-mat) + if firedrake_to_meshmode: + new_order = order[cell_node_list] + # else just need to reshape new_order so that can apply flip-mat + else: + new_order = order.reshape( + (order.shape[0]//nunit_nodes, nunit_nodes) + order.shape[1:]) + + flip_mat = finat_element_importer.flip_matrix() + reorder_nodes(mesh_importer.orientations(), new_order, flip_mat, + unflip=firedrake_to_meshmode) + new_order = new_order.flatten() + + # Resize new_order if going meshmode->firedrake and meshmode + # has duplicate nodes (e.g if used a CG fspace) + # + # TODO: this is done VERY LAZILY (NOT GOOD) + if not firedrake_to_meshmode and num_fd_nodes != num_mm_nodes: + newnew_order = np.zeros(num_fd_nodes, dtype=np.int32) + pyt_ndx = 0 + for nodes in cell_node_list: + for fd_index in nodes: + newnew_order[fd_index] = new_order[pyt_ndx] + pyt_ndx += 1 + + new_order = newnew_order + + # }}} + + return new_order + + +@cached +def get_factory(mesh_importer, degree): + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + return InterpolatoryQuadratureSimplexGroupFactory(degree) + + +@cached +def get_discretization(mesh_importer, key): + finat_element_importer, cl_ctx = key + assert isinstance(finat_element_importer, FinatLagrangeElementImporter) + + degree = finat_element_importer.data.degree + discretization = Discretization(cl_ctx, + mesh_importer.meshmode_mesh(), + get_factory(mesh_importer, degree)) + + return discretization + + +class FiredrakeFunctionSpaceDataImporter(ExternalImportHandler): + """ + This is not *quite* the usual thought of a + :class:`ExternalImportHandler`. + + It handles an "import" in the sense that a mesh & finat element + define a lot of the data of a function space that + can be shared between multiple function spaces on the + same mesh, so we use + this object to store data based on the mesh and finat element, + cached on the mesh geometry. + In contrast to firedrake's shared data, we need the mesh + geometry + + .. attribute:: data + + A tuple whose first entry is an instance of + class :class:`FiredrakeMeshGeometryImporter` and second + entry is an instance of class + :class:`FinatLagrangeElementImporter` + """ + # FIXME: Give two finat elts + def __init__(self, cl_ctx, mesh_importer, finat_element_importer): + """ + Note that :param:`mesh_importer` and :param:`finat_element_importer` + are used for checking + """ + if mesh_importer.topological_importer == mesh_importer: + raise TypeError(":param:`mesh_importer` is a " + "FiredrakeMeshTopologyImporter, but " + " must be a FiredrakeMeshGeometryImporter") + importer = (mesh_importer.data, finat_element_importer.data) + super(FiredrakeFunctionSpaceDataImporter, self).__init__(importer) + + self._fspace_data = FunctionSpaceData(mesh_importer.data, + finat_element_importer.data) + + self._cl_ctx = cl_ctx + self._mesh_importer = mesh_importer + self._finat_element_importer = finat_element_importer + self._discretization = None + + def firedrake_to_meshmode(self): + """ + See :func:`reordering_array` + """ + return reordering_array(self._mesh_importer, + (self._finat_element_importer, True), + self._fspace_data) + + def meshmode_to_firedrake(self): + """ + See :func:`reordering_array` + """ + return reordering_array(self._mesh_importer, + (self._finat_element_importer, False), + self._fspace_data) + + def discretization(self): + """ + Creates a discretization on the mesh returned by this object's + mesh importer's :func:`meshmode_mesh()` + """ + return get_discretization(self._mesh_importer, + (self._finat_element_importer, self._cl_ctx)) + + def factory(self): + """ + Returns a :mod:`meshmode` + :class:`InterpolatoryQuadratureSimplexGroupFactory` + of the same degree as this object's finat element. + """ + degree = self._finat_element_importer.data.degree + return get_factory(self._mesh_importer, degree) diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py new file mode 100644 index 0000000000000000000000000000000000000000..da965cbe4148b7f264622b5a83b99d91d930bd11 --- /dev/null +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -0,0 +1,461 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 warnings import warn # noqa +from collections import defaultdict +import numpy as np + +from meshmode.interop import ExternalImportHandler +from meshmode.interop.firedrake.function_space_coordless import \ + FiredrakeCoordinatelessFunctionImporter + + +class FiredrakeMeshGeometryImporter(ExternalImportHandler): + """ + This takes a :mod:`firedrake` :class:`MeshGeometry` + and converts its data so that :mod:`meshmode` can handle it. + + .. attribute:: data + + A :mod:`firedrake` :class:`MeshGeometry` instance + """ + + def __init__(self, + mesh, + coordinates_importer, + normals=None, + no_normals_warn=True): + """ + :param mesh: A :mod:`firedrake` :class:`MeshGeometry`. + We require that :aram:`mesh` have co-dimesnion + of 0 or 1. + Moreover, if :param:`mesh` is a 2-surface embedded in 3-space, + we _require_ that :function:`init_cell_orientations` + has been called already. + + :param coordinates_importer: An instance of + class :class:`FiredrakeCoordinatelessFunctionImporter` + to use, mapping the points of the topological mesh + to their coordinates (The function is coordinateless + in the sense that its *domain* has no coordinates) + + For other params see :meth:`orientations` + """ + super(FiredrakeMeshGeometryImporter, self).__init__(mesh) + + # {{{ Make sure input data is valid + + # Ensure is not a topological mesh + if mesh.topological == mesh: + raise TypeError(":param:`mesh` must be of type" + " :class:`firedrake.mesh.MeshGeometry`") + + # Ensure dimensions are in appropriate ranges + supported_dims = [1, 2, 3] + if mesh.geometric_dimension() not in supported_dims: + raise ValueError("Geometric dimension is %s. Geometric " + " dimension must be one of range %s" + % (mesh.geometric_dimension(), supported_dims)) + + # Raise warning if co-dimension is not 0 or 1 + co_dimension = mesh.geometric_dimension() - mesh.topological_dimension() + if co_dimension not in [0, 1]: + raise ValueError("Codimension is %s, but must be 0 or 1." % + (co_dimension)) + + # Ensure coordinates are coordinateless + if not isinstance(coordinates_importer, + FiredrakeCoordinatelessFunctionImporter): + raise ValueError(":param:`coordinates_importer` must be of type" + " FiredrakeCoordinatelessFunctionImporter") + + fspace_importer = coordinates_importer.function_space_importer() + topology_importer = fspace_importer.mesh_importer() + + if topology_importer.data != mesh.topology: + raise ValueError("Topology :param:`coordinates` lives on must be " + "the same " + "topology that :param:`mesh` lives on") + + # }}} + + # For sharing data like in firedrake + self._shared_data_cache = defaultdict(dict) + + # Store input information + self._coordinates_importer = coordinates_importer + self._topology_importer = topology_importer + + self._normals = normals + self._no_normals_warn = no_normals_warn + + # To be computed later + self._vertex_indices = None + self._vertices = None + self._nodes = None + self._group = None + self._orient = None + self._facial_adjacency_groups = None + self._meshmode_mesh = None + + def callback(cl_ctx): + """ + Finish initialization by creating a coordinates function importer + for public access on this mesh which is not "coordinateless" (i.e. + its domain has coordinates) + """ + from meshmode.interop.firedrake.function import \ + FiredrakeFunctionImporter + from meshmode.interop.firedrake.function_space import \ + FiredrakeWithGeometryImporter + from firedrake import Function + + coordinates_fs = self.data.coordinates.function_space() + coordinates_fs_importer = \ + self._coordinates_importer.function_space_importer() + + fspace_importer = \ + FiredrakeWithGeometryImporter(cl_ctx, + coordinates_fs, + coordinates_fs_importer, + self) + f = Function(fspace_importer.data, val=self._coordinates_importer.data) + self._coordinates_function_importer = \ + FiredrakeFunctionImporter(f, fspace_importer) + + del self._callback + + self._callback = callback + + def initialized(self): + return not hasattr(self, '_callback') + + def init(self, cl_ctx): + if not self.initialized(): + self._callback(cl_ctx) + + def __getattr__(self, attr): + """ + If can't find an attribute, look in the underlying + topological mesh + (Done like :class:`firedrake.function.MeshGeometry`) + """ + return getattr(self._topology_importer, attr) + + @property + def coordinates_importer(self): + """ + Return coordinates as a function + + PRECONDITION: Have called :meth:`init` + """ + try: + return self._coordinates_function_importer + except AttributeError: + raise AttributeError("No coordinates function, have you finished" + " initializing this object" + " (i.e. have you called :meth:`init`)?") + + def _compute_vertex_indices_and_vertices(self): + if self._vertex_indices is None: + fspace_importer = self.coordinates_importer.function_space_importer() + finat_element_importer = fspace_importer.finat_element_importer + + from finat.fiat_elements import Lagrange + from finat.spectral import GaussLobattoLegendre + if not isinstance(finat_element_importer.data, + (Lagrange, GaussLobattoLegendre)): + raise TypeError("Coordinates must live in a continuous space " + " (Lagrange or GaussLobattoLegendre), not %s" + % type(finat_element_importer.data)) + + # Convert cell node list of mesh to vertex list + unit_vertex_indices = finat_element_importer.unit_vertex_indices() + cfspace = self.data.coordinates.function_space() + if self.icell_to_fd is not None: + cell_node_list = cfspace.cell_node_list[self.icell_to_fd] + else: + cell_node_list = cfspace.cell_node_list + + vertex_indices = cell_node_list[:, unit_vertex_indices] + + # Get maps newnumbering->old and old->new (new numbering comes + # from removing the non-vertex + # nodes) + vert_ndx_to_fd_ndx = np.unique(vertex_indices.flatten()) + fd_ndx_to_vert_ndx = dict(zip(vert_ndx_to_fd_ndx, + np.arange(vert_ndx_to_fd_ndx.shape[0], + dtype=np.int32) + )) + # Get vertices array + vertices = np.real( + self.data.coordinates.dat.data[vert_ndx_to_fd_ndx]) + + #:mod:`meshmode` wants shape to be [ambient_dim][nvertices] + if len(vertices.shape) == 1: + # 1 dim case, (note we're about to transpose) + vertices = vertices.reshape(vertices.shape[0], 1) + vertices = vertices.T.copy() + + # Use new numbering on vertex indices + vertex_indices = np.vectorize(fd_ndx_to_vert_ndx.get)(vertex_indices) + + # store vertex indices and vertices + self._vertex_indices = vertex_indices + self._vertices = vertices + + def vertex_indices(self): + """ + A numpy array of shape *(nelements, nunitvertices)* + holding the vertex indices associated to each element + """ + self._compute_vertex_indices_and_vertices() + return self._vertex_indices + + def vertices(self): + """ + Return the mesh vertices as a numpy array of shape + *(dim, nvertices)* + """ + self._compute_vertex_indices_and_vertices() + return self._vertices + + def nodes(self): + """ + Return the mesh nodes as a numpy array of shape + *(dim, nelements, nunitnodes)* + """ + if self._nodes is None: + coords = self.data.coordinates.dat.data + cfspace = self.data.coordinates.function_space() + + if self.icell_to_fd is not None: + cell_node_list = cfspace.cell_node_list[self.icell_to_fd] + else: + cell_node_list = cfspace.cell_node_list + self._nodes = np.real(coords[cell_node_list]) + + # reshape for 1D so that [nelements][nunit_nodes][dim] + if len(self._nodes.shape) != 3: + self._nodes = np.reshape(self._nodes, self._nodes.shape + (1,)) + + # Change shape to [dim][nelements][nunit_nodes] + self._nodes = np.transpose(self._nodes, (2, 0, 1)) + + return self._nodes + + def _get_unflipped_group(self): + """ + Return an instance of :class:meshmode.mesh.SimplexElementGroup` + but with elements which are not guaranteed to have a positive + orientation. + + This is important because :meth:`group` requires + the result of :meth:`orientations` to guarantee each + element has positive orientations, and + :meth:`orientations` (usually) needs an unflipped group to compute + the orientations, so we need this function to prevent a recursive + loop of (orientations calling group calling orientations calling...etc.) + """ + # Cache the group we create since :meth:`group` and :meth:`orientations` + # may both use it. One should note that once :meth:`group` terminates, + # this attribute is deleted (once :attr:`_group` is computed then + # :attr:`_orient` must also have been computed, so we no longer have any + # use for :attr:`_unflipped_group`) + if not hasattr(self, "_unflipped_group"): + from meshmode.mesh import SimplexElementGroup + + fspace_importer = self.coordinates_importer.function_space_importer() + finat_element_importer = fspace_importer.finat_element_importer + + self._unflipped_group = SimplexElementGroup( + finat_element_importer.data.degree, + self.vertex_indices(), + self.nodes(), + dim=self.cell_dimension(), + unit_nodes=finat_element_importer.unit_nodes()) + + return self._unflipped_group + + def group(self): + """ + Return an instance of :class:meshmode.mesh.SimplexElementGroup` + corresponding to the mesh :attr:`data` + """ + if self._group is None: + from meshmode.mesh.processing import flip_simplex_element_group + self._group = flip_simplex_element_group(self.vertices(), + self._get_unflipped_group(), + self.orientations() < 0) + # We don't need this anymore + del self._unflipped_group + + return self._group + + def orientations(self): + """ + Return the orientations of the mesh elements: + an array, the *i*th element is > 0 if the *ith* element + is positively oriented, < 0 if negatively oriented + + :param normals: _Only_ used if :param:`mesh` is a 1-surface + embedded in 2-space. In this case, + - If *None* then + all elements are assumed to be positively oriented. + - Else, should be a list/array whose *i*th entry + is the normal for the *i*th element (*i*th + in :param:`mesh`*.coordinate.function_space()*'s + :attribute:`cell_node_list`) + + :param no_normals_warn: If *True*, raises a warning + if :param:`mesh` is a 1-surface embedded in 2-space + and :param:`normals` is *None*. + """ + if self._orient is None: + # compute orientations + tdim = self.data.topological_dimension() + gdim = self.data.geometric_dimension() + + orient = None + if gdim == tdim: + # We use :mod:`meshmode` to check our orientations + from meshmode.mesh.processing import \ + find_volume_mesh_element_group_orientation + + orient = find_volume_mesh_element_group_orientation( + self.vertices(), + self._get_unflipped_group()) + + if tdim == 1 and gdim == 2: + # In this case we have a 1-surface embedded in 2-space + orient = np.ones(self.nelements()) + if self._normals: + for i, (normal, vertices) in enumerate(zip( + np.array(self._normals), self.vertices())): + if np.cross(normal, vertices) < 0: + orient[i] = -1.0 + elif self._no_normals_warn: + warn("Assuming all elements are positively-oriented.") + + elif tdim == 2 and gdim == 3: + # In this case we have a 2-surface embedded in 3-space + orient = self.data.cell_orientations().dat.data + r""" + Convert (0 \implies negative, 1 \implies positive) to + (-1 \implies negative, 1 \implies positive) + """ + orient *= 2 + orient -= np.ones(orient.shape, dtype=orient.dtype) + + self._orient = orient + #Make sure the mesh fell into one of the cases + """ + NOTE : This should be guaranteed by previous checks, + but is here anyway in case of future development. + """ + assert self._orient is not None + + return self._orient + + def face_vertex_indices_to_tags(self): + """ + Returns a dict mapping frozensets + of vertex indices (which identify faces) to a + list of boundary tags + """ + finat_element = self.data.coordinates.function_space().finat_element + exterior_facets = self.data.exterior_facets + + # fvi_to_tags maps frozenset(vertex indices) to tags + fvi_to_tags = {} + # maps faces to local vertex indices + connectivity = finat_element.cell.connectivity[(self.cell_dimension()-1, 0)] + + # Compatability for older versions of firedrake + try: + local_fac_number = exterior_facets.local_facet_number + except AttributeError: + local_fac_number = exterior_facets.local_facet_dat.data + + for i, (icells, ifacs) in enumerate(zip(exterior_facets.facet_cell, + local_fac_number)): + # Compatability for older versions of firedrake + try: + iter(ifacs) + except TypeError: + ifacs = [ifacs] + + for icell, ifac in zip(icells, ifacs): + # If necessary, convert to new cell numbering + if self.fd_to_icell is not None: + if icell not in self.fd_to_icell: + continue + else: + icell = self.fd_to_icell[icell] + + # record face vertex indices to tag map + cell_vertices = self.vertex_indices()[icell] + facet_indices = connectivity[ifac] + fvi = frozenset(cell_vertices[list(facet_indices)]) + fvi_to_tags.setdefault(fvi, []) + fvi_to_tags[fvi].append(exterior_facets.markers[i]) + + # }}} + + return fvi_to_tags + + def facial_adjacency_groups(self): + """ + Return a list of :mod:`meshmode` :class:`FacialAdjacencyGroups` + as used in the construction of a :mod:`meshmode` :class:`Mesh` + """ + # {{{ Compute facial adjacency groups if not already done + + if self._facial_adjacency_groups is None: + from meshmode.mesh import _compute_facial_adjacency_from_vertices + + self._facial_adjacency_groups = _compute_facial_adjacency_from_vertices( + [self.group()], + self.bdy_tags(), + np.int32, np.int8, + face_vertex_indices_to_tags=self.face_vertex_indices_to_tags()) + + # }}} + + return self._facial_adjacency_groups + + def meshmode_mesh(self): + """ + PRECONDITION: Have called :meth:`init` + """ + if self._meshmode_mesh is None: + assert self.initialized(), \ + "Must call :meth:`init` before :meth:`meshmode_mesh`" + + from meshmode.mesh import Mesh + self._meshmode_mesh = \ + Mesh(self.vertices(), [self.group()], + boundary_tags=self.bdy_tags(), + nodal_adjacency=self.nodal_adjacency(), + facial_adjacency_groups=self.facial_adjacency_groups()) + + return self._meshmode_mesh diff --git a/meshmode/interop/firedrake/mesh_topology.py b/meshmode/interop/firedrake/mesh_topology.py new file mode 100644 index 0000000000000000000000000000000000000000..3a5d7a2c1a68b3313c2f4d4c114d3f44dc5b5411 --- /dev/null +++ b/meshmode/interop/firedrake/mesh_topology.py @@ -0,0 +1,221 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 warnings import warn # noqa +import numpy as np + +from meshmode.interop import ExternalImportHandler + + +# {{{ ImportHandler for firedrake's MeshTopology class + +class FiredrakeMeshTopologyImporter(ExternalImportHandler): + """ + An Importer for :class:`firedrake.mesh.MeshTopology`. + Holds the topological (as opposed to geometric) information + about a mesh. + """ + + def __init__(self, mesh, cells_to_use=None): + """ + :param mesh: An instance :mod:`firedrake` :class:`MeshTopology` or + :class:`MeshGeometry`. If an instance of + :class:`MeshGeometry`, uses its underlying topology. + :param cells_to_use: Either + + * *None*, in which case this argument is ignored + * An array of cell ids, in which case those are the + only cells for which information is gathered/converted + + We require that :param:`mesh` have co-dimesnion + of 0 or 1. + Moreover, if :param:`mesh` is a 2-surface embedded in 3-space, + we _require_ that :function:`init_cell_orientations` + has been called already. + + :raises TypeError: If :param:`mesh` is not of :mod:`firedrake` + :class:`MeshTopology` or :class:`MeshGeometry` + """ + top = mesh.topological # convert geometric to topological + + # {{{ Check input types + from firedrake.mesh import MeshTopology + if not isinstance(top, MeshTopology): + raise TypeError(":param:`mesh` must be of type " + ":class:`firedrake.mesh.MeshTopology` or " + ":class:`firedrake.mesh.MeshGeometry`") + # }}} + + super(FiredrakeMeshTopologyImporter, self).__init__(top) + + # Ensure has simplex-type elements + if not top.ufl_cell().is_simplex(): + raise ValueError(":param:`mesh` must have simplex type elements, " + "%s is not a simplex" % (mesh.ufl_cell())) + + # Ensure dimensions are in appropriate ranges + supported_dims = [1, 2, 3] + if self.cell_dimension() not in supported_dims: + raise ValueError("Cell dimension is %s. Cell dimension must be one of" + " %s" % (self.cell_dimension(), supported_dims)) + + self._nodal_adjacency = None + self.icell_to_fd = cells_to_use # Map cell index -> fd cell index + self.fd_to_icell = None # Map fd cell index -> cell index + + # perform checks on :param:`cells_to_use` if not *None* + if self.icell_to_fd is not None: + assert np.unique(self.icell_to_fd).shape == self.icell_to_fd.shape + self.fd_to_icell = dict(zip(self.icell_to_fd, + np.arange(self.icell_to_fd.shape[0], + dtype=np.int32) + )) + + @property + def topology_importer(self): + """ + A reference to *self*, for compatability with mesh geometry + importers + """ + return self + + @property + def topological_importer(self): + """ + A reference to *self*, for compatability with mesh geometry + importers + """ + return self + + def cell_dimension(self): + """ + Return the dimension of the cells used by this topology + """ + return self.data.cell_dimension() + + def nelements(self): + """ + Return the number of cells in this topology + """ + if self.icell_to_fd is None: + num_cells = self.data.num_cells() + else: + num_cells = self.icell_to_fd.shape[0] + + return num_cells + + def nunit_vertices(self): + """ + Return the number of unit vertices on the reference element + """ + return self.data.ufl_cell().num_vertices() + + def bdy_tags(self): + """ + Return a tuple of bdy tags as requested in + the construction of a :mod:`meshmode` :class:`Mesh` + + The tags used are :class:`meshmode.mesh.BTAG_ALL`, + :class:`meshmode.mesh.BTAG_REALLY_ALL`, and + any markers in the mesh topology's exterior facets + (see :attr:`firedrake.mesh.MeshTopology.exterior_facets.unique_markers`) + """ + from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL + bdy_tags = [BTAG_ALL, BTAG_REALLY_ALL] + + unique_markers = self.data.exterior_facets.unique_markers + if unique_markers is not None: + bdy_tags += list(unique_markers) + + return tuple(bdy_tags) + + def nodal_adjacency(self): + """ + Returns a :class:`meshmode.mesh.NodalAdjacency` object + representing the nodal adjacency of this mesh + """ + if self._nodal_adjacency is None: + # TODO... not sure how to get around the private access + plex = self.data._plex + + # dmplex cell Start/end and vertex Start/end. + c_start, c_end = plex.getHeightStratum(0) + v_start, v_end = plex.getDepthStratum(0) + + # TODO... not sure how to get around the private access + to_fd_id = np.vectorize(self.data._cell_numbering.getOffset)( + np.arange(c_start, c_end, dtype=np.int32)) + + element_to_neighbors = {} + verts_checked = set() # dmplex ids of vertex checked + + # If using all cells, loop over them all + if self.icell_to_fd is None: + range_ = range(c_start, c_end) + # Otherwise, just the ones you're using + else: + isin = np.isin(to_fd_id, self.icell_to_fd) + range_ = np.arange(c_start, c_end, dtype=np.int32)[isin] + + # For each cell + for cell_id in range_: + # For each vertex touching the cell (that haven't already seen) + for vert_id in plex.getTransitiveClosure(cell_id)[0]: + if v_start <= vert_id < v_end and vert_id not in verts_checked: + verts_checked.add(vert_id) + cells = [] + # Record all cells touching that vertex + support = plex.getTransitiveClosure(vert_id, + useCone=False)[0] + for other_cell_id in support: + if c_start <= other_cell_id < c_end: + cells.append(to_fd_id[other_cell_id - c_start]) + + # If only using some cells, clean out extraneous ones + # and relabel them to new id + cells = set(cells) + if self.fd_to_icell is not None: + cells = set([self.fd_to_icell[fd_ndx] + for fd_ndx in cells + if fd_ndx in self.fd_to_icell]) + + # mark cells as neighbors + for cell_one in cells: + element_to_neighbors.setdefault(cell_one, set()) + element_to_neighbors[cell_one] |= cells + + # Create neighbors_starts and neighbors + neighbors = [] + neighbors_starts = np.zeros(self.nelements() + 1, dtype=np.int32) + for iel in range(len(element_to_neighbors)): + elt_neighbors = element_to_neighbors[iel] + neighbors += list(elt_neighbors) + neighbors_starts[iel+1] = len(neighbors) + + neighbors = np.array(neighbors, dtype=np.int32) + + from meshmode.mesh import NodalAdjacency + self._nodal_adjacency = NodalAdjacency(neighbors_starts=neighbors_starts, + neighbors=neighbors) + return self._nodal_adjacency + +# }}} diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py new file mode 100644 index 0000000000000000000000000000000000000000..f3337c4a831c4bacbb2478358648b152d8c889d1 --- /dev/null +++ b/test/test_firedrake_interop.py @@ -0,0 +1,233 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__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 + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl + as pytest_generate_tests) + +from meshmode.interop.firedrake import FromFiredrakeConnection + +import pytest + +import logging +logger = logging.getLogger(__name__) + +# skip testing this module if cannot import firedrake +firedrake = pytest.importorskip("firedrake") + +from firedrake import ( + UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh, + FunctionSpace, VectorFunctionSpace, Function, + SpatialCoordinate, Constant) + + +CLOSE_ATOL = 10**-12 + + +@pytest.fixture(params=[1, 2, 3], ids=["1D", "2D", "3D"]) +def fdrake_mesh(request): + dim = request.param + if dim == 1: + return UnitIntervalMesh(100) + if dim == 2: + return UnitSquareMesh(10, 10) + if dim == 3: + return UnitCubeMesh(5, 5, 5) + return None + + +@pytest.fixture(params=["CG", "DG"]) +def fdrake_family(request): + return request.param + + +@pytest.fixture(params=[1, 2, 3], ids=["P^1", "P^2", "P^3"]) +def fdrake_degree(request): + return request.param + + +# {{{ Basic conversion checks for the function space + +def test_discretization_consistency(ctx_factory, fdrake_mesh, fdrake_degree): + """ + While nodes may change, vertex conversion should be *identical* up to + reordering, ensure this is the case for DG spaces. Also ensure the + meshes have the same basic properties and the function space/discretization + agree across firedrake vs meshmode + """ + # get fdrake_verts (shaped like (nverts, dim)) + fdrake_verts = fdrake_mesh.coordinates.dat.data + if fdrake_mesh.geometric_dimension() == 1: + fdrake_verts = fdrake_verts[:, np.newaxis] + + # Get meshmode vertices (shaped like (dim, nverts)) + fdrake_fspace = FunctionSpace(fdrake_mesh, 'DG', fdrake_degree) + cl_ctx = ctx_factory() + fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_fspace) + to_discr = fdrake_connection.to_discr() + meshmode_verts = to_discr.mesh.vertices + + # Ensure the meshmode mesh has one group and make sure both + # meshes agree on some basic properties + assert len(to_discr.mesh.groups) == 1 + fdrake_mesh_fspace = fdrake_mesh.coordinates.function_space() + fdrake_mesh_order = fdrake_mesh_fspace.finat_element.degree + assert to_discr.mesh.groups[0].order == fdrake_mesh_order + assert to_discr.mesh.groups[0].nelements == fdrake_mesh.num_cells() + assert to_discr.mesh.nvertices == fdrake_mesh.num_vertices() + + # Ensure that the vertex sets are identical up to reordering + # Nb: I got help on this from stack overflow: + # https://stackoverflow.com/questions/38277143/sort-2d-numpy-array-lexicographically # noqa: E501 + lex_sorted_mm_verts = meshmode_verts[:, np.lexsort(meshmode_verts)] + lex_sorted_fdrake_verts = fdrake_verts[np.lexsort(fdrake_verts.T)] + np.testing.assert_array_equal(lex_sorted_mm_verts, lex_sorted_fdrake_verts.T) + + # Ensure the discretization and the firedrake function space agree on + # some basic properties + finat_elt = fdrake_fspace.finat_element + assert len(to_discr.groups) == 1 + assert to_discr.groups[0].order == finat_elt.degree + assert to_discr.groups[0].nunit_nodes == finat_elt.space_dimension() + assert to_discr.nnodes == fdrake_fspace.node_count + +# }}} + + +# {{{ Double check functions are being transported correctly + +def alternating_sum_fd(spat_coord): + """ + Return an expression x1 - x2 + x3 -+... + """ + return sum( + [(-1)**i * spat_coord for i, spat_coord in enumerate(spat_coord)] + ) + + +def alternating_sum_mm(nodes): + """ + Take the *(dim, nnodes)* array nodes and return an array + holding the alternating sum of the coordinates of each node + """ + alternator = np.ones(nodes.shape[0]) + alternator[1::2] *= -1 + return np.matmul(alternator, nodes) + + +# In 1D/2D/3D check constant 1, +# projection to x1, x1/x1+x2/x1+x2+x3, and x1/x1-x2/x1-x2+x3. +# This should show that projection to any coordinate in 1D/2D/3D +# transfers correctly. +test_functions = [ + (lambda spat_coord: Constant(1.0), lambda nodes: np.ones(nodes.shape[1])), + (lambda spat_coord: spat_coord[0], lambda nodes: nodes[0, :]), + (sum, lambda nodes: np.sum(nodes, axis=0)), + (alternating_sum_fd, alternating_sum_mm) +] + + +@pytest.mark.parametrize("fdrake_f_expr,meshmode_f_eval", test_functions) +def test_function_transfer(ctx_factory, + fdrake_mesh, fdrake_family, fdrake_degree, + fdrake_f_expr, meshmode_f_eval): + """ + Make sure creating a function then transporting it is the same + (up to resampling error) as creating a function on the transported + mesh + """ + fdrake_fspace = FunctionSpace(fdrake_mesh, fdrake_family, fdrake_degree) + spat_coord = SpatialCoordinate(fdrake_mesh) + + fdrake_f = Function(fdrake_fspace).interpolate(fdrake_f_expr(spat_coord)) + + cl_ctx = ctx_factory() + fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_fspace) + transported_f = fdrake_connection.from_firedrake(fdrake_f) + + to_discr = fdrake_connection.to_discr() + with cl.CommandQueue(cl_ctx) as queue: + nodes = to_discr.nodes().get(queue=queue) + meshmode_f = meshmode_f_eval(nodes) + + np.testing.assert_allclose(transported_f, meshmode_f, atol=CLOSE_ATOL) + +# }}} + + +# {{{ Idempotency tests fd->mm->fd and (fd->)mm->fd->mm for connection + +def check_idempotency(fdrake_connection, fdrake_function): + """ + Make sure fd->mm->fd and mm->fd->mm are identity for DG spaces + """ + fdrake_fspace = fdrake_connection.from_function_space() + + # Test for idempotency fd->mm->fd + mm_field = fdrake_connection.from_firedrake(fdrake_function) + fdrake_function_copy = Function(fdrake_fspace) + fdrake_connection.from_meshmode(mm_field, fdrake_function_copy) + + np.testing.assert_allclose(fdrake_function_copy.dat.data, + fdrake_function.dat.data, + atol=CLOSE_ATOL) + + # Test for idempotency (fd->)mm->fd->mm + mm_field_copy = fdrake_connection.from_firedrake(fdrake_function_copy) + np.testing.assert_allclose(mm_field_copy, mm_field, atol=CLOSE_ATOL) + + +def test_scalar_idempotency(ctx_factory, fdrake_mesh, fdrake_degree): + """ + Make sure fd->mm->fd and mm->fd->mm are identity for scalar DG spaces + """ + fdrake_fspace = FunctionSpace(fdrake_mesh, 'DG', fdrake_degree) + + # Make a function with unique values at each node + fdrake_unique = Function(fdrake_fspace) + fdrake_unique.dat.data[:] = np.arange(fdrake_unique.dat.data.shape[0]) + + # test idempotency + cl_ctx = ctx_factory() + fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_fspace) + check_idempotency(fdrake_connection, fdrake_unique) + + +def test_vector_idempotency(ctx_factory, fdrake_mesh, fdrake_degree): + """ + Make sure fd->mm->fd and mm->fd->mm are identity for vector DG spaces + """ + fdrake_vfspace = VectorFunctionSpace(fdrake_mesh, 'DG', fdrake_degree) + + # Make a function with unique values at each node + xx = SpatialCoordinate(fdrake_vfspace.mesh()) + fdrake_unique = Function(fdrake_vfspace).interpolate(xx) + + # test idempotency + cl_ctx = ctx_factory() + fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_vfspace) + check_idempotency(fdrake_connection, fdrake_unique) + +# }}}