diff --git a/meshmode/interop/firedrake/__init__.py b/meshmode/interop/firedrake/__init__.py index 33daeb392b3276f855867ff79daaaea0ab3e2f64..890dbaf7db7c6025829120c983b199ce32d5783f 100644 --- a/meshmode/interop/firedrake/__init__.py +++ b/meshmode/interop/firedrake/__init__.py @@ -20,16 +20,12 @@ 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 +from meshmode.interop.firedrake.connection import FromFiredrakeConnection +from meshmode.interop.firedrake.mesh import import_firedrake_mesh - -# {{{ Helper functions to construct importers for firedrake objects +__all__ = ["FromFiredrakeConnection", "import_firedrake_mesh"] def _compute_cells_near_bdy(mesh, bdy_id): @@ -45,211 +41,3 @@ def _compute_cells_near_bdy(mesh, bdy_id): 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/connection.py b/meshmode/interop/firedrake/connection.py index 1e695e747e04cd0e6df80ac3996fe5c74aa7fd44..e3d950d380198582d2a3b231fcebdcc2d01e5a70 100644 --- a/meshmode/interop/firedrake/connection.py +++ b/meshmode/interop/firedrake/connection.py @@ -73,7 +73,6 @@ def _reorder_nodes(orient, nodes, flip_matrix, unflip=False): class FromFiredrakeConnection: - # FIXME : docs """ A connection created from a :mod:`firedrake` ``"CG"`` or ``"DG"`` function space which creates a corresponding @@ -83,10 +82,16 @@ class FromFiredrakeConnection: .. attribute:: to_discr The discretization corresponding to the firedrake function - space (DESCRIBE HERE) + space created with a + :class:`InterpolatoryQuadratureSimplexElementGroup`. """ def __init__(self, cl_ctx, fdrake_fspace): - # FIXME : docs + """ + :param cl_ctx: A :mod:`pyopencl` computing context + :param fdrake_fspace: A :mod:`firedrake` ``"CG"`` or ``"DG"`` + function space (of class :class:`WithGeometry`) built on + a mesh which is importable by :func:`import_firedrake_mesh`. + """ # Ensure fdrake_fspace is a function space with appropriate reference # element. from firedrake.functionspaceimpl import WithGeometry @@ -176,7 +181,19 @@ class FromFiredrakeConnection: return self._fspace_cache[dim] def from_firedrake(self, function, out=None): - """transport fiiredrake function""" + """ + transport firedrake function onto :attr:`to_discr` + + :param function: A :mod:`firedrake` function to transfer onto + :attr:`to_discr`. Its function space must have + the same family, degree, and mesh as ``self.from_fspace()``. + :param out: If *None* then ignored, otherwise a numpy array of the + shape *function.dat.data.shape.T* (i.e. + *(dim, nnodes)* or *(nnodes,)* in which :param:`function`'s + transported data is stored. + + :return: a numpy array holding the transported function + """ # make sure function is a firedrake function in an appropriate # function space from firedrake.function import Function @@ -212,7 +229,19 @@ class FromFiredrakeConnection: return out def from_meshmode(self, mm_field, out=None): - """transport meshmode field""" + """ + transport meshmode field from :attr:`to_discr` into an + appropriate firedrake function space. + + :param mm_field: A numpy array of shape *(nnodes,)* or *(dim, nnodes)* + representing a function on :attr:`to_distr`. + :param out: If *None* then ignored, otherwise a :mod:`firedrake` + function of the right function space for the transported data + to be stored in. + + :return: a :mod:`firedrake` :class:`Function` holding the transported + data. + """ if self._ufl_element.family() == 'Lagrange': raise ValueError("Cannot convert functions from discontinuous " " space (meshmode) to continuous firedrake " diff --git a/meshmode/interop/firedrake/function.py b/meshmode/interop/firedrake/function.py deleted file mode 100644 index 14ef3a4638956bde35798c32ead4d3a2a689dd3c..0000000000000000000000000000000000000000 --- a/meshmode/interop/firedrake/function.py +++ /dev/null @@ -1,118 +0,0 @@ -__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 deleted file mode 100644 index 85b3d6c22645f57ac2a14dff0e953307222b659f..0000000000000000000000000000000000000000 --- a/meshmode/interop/firedrake/function_space.py +++ /dev/null @@ -1,206 +0,0 @@ -__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 deleted file mode 100644 index e6fdec6cd2137c0831340bb13db206da3e21f9e5..0000000000000000000000000000000000000000 --- a/meshmode/interop/firedrake/function_space_coordless.py +++ /dev/null @@ -1,221 +0,0 @@ -__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 deleted file mode 100644 index c7b316a7b3da443c8647885ac0badee0eef945ad..0000000000000000000000000000000000000000 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ /dev/null @@ -1,265 +0,0 @@ -__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.py b/meshmode/interop/firedrake/mesh.py index 69e294478afa0bb0f39e20f06f927420fd9b9d83..26dbed764613587316d55886824c769edc330ea4 100644 --- a/meshmode/interop/firedrake/mesh.py +++ b/meshmode/interop/firedrake/mesh.py @@ -21,8 +21,7 @@ THE SOFTWARE. """ __doc__ = """ -.. autofunction:: get_firedrake_nodal_adjacency_group -.. autofunction:: get_firedrake_vertex_indices +.. autofunction:: import_firedrake_mesh """ from warnings import warn # noqa @@ -34,9 +33,25 @@ import six def _get_firedrake_nodal_info(fdrake_mesh_topology): - # FIXME: do docs """ - Get nodal adjacency, and vertex indices from a firedrake mesh topology + Get nodal adjacency and vertex indices corresponding + to a firedrake mesh topology. Note that as we do not use + geometric information, there is no guarantee that elements + have a positive orientation. + + The elements (in firedrake lingo, the cells) + are guaranteed to have the same numbering in :mod:`meshmode` + as :mod:`firedrdake` + + :param fdrake_mesh_topology: A :mod:`firedrake` instance of class + :class:`MeshTopology` or :class:`MeshGeometry`. + + :return: Returns *vertex_indices* as a numpy array of shape + *(nelements, ref_element.nvertices)* (as described by + the ``vertex_indices`` attribute of a :class:`MeshElementGroup`) + and a :class:`NodalAdjacency` constructed from + :param:`fdrake_mesh_topology` + as a tuple *(vertex_indices, nodal_adjacency)*. """ top = fdrake_mesh_topology.topology @@ -147,9 +162,19 @@ def _get_firedrake_boundary_tags(fdrake_mesh): return tuple(bdy_tags) -def _get_firedrake_facial_adjacency_groups(fdrake_mesh): - # FIXME: do docs - top = fdrake_mesh.topology +def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology): + """ + Return facial_adjacency_groups corresponding to + the given firedrake mesh topology. Note that as we do not + have geometric information, elements may need to be + flipped later. + + :param fdrake_mesh_topology: A :mod:`firedrake` instance of class + :class:`MeshTopology` or :class:`MeshGeometry`. + :return: A list of maps to :class:`FacialAdjacencyGroup`s as required + by a :mod:`meshmode` :class:`Mesh` + """ + top = fdrake_mesh_topology.topology # We only need one group # for interconnectivity and one for boundary connectivity. # The tricky part is moving from firedrake local facet numbering @@ -221,8 +246,8 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh): from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL ext_neighbors = np.zeros(ext_elements.shape, dtype=np.int32) for ifac, marker in enumerate(top.exterior_facets.markers): - ext_neighbors[ifac] = -(boundary_tag_bit(BTAG_ALL) \ - | boundary_tag_bit(BTAG_REALLY_ALL) \ + ext_neighbors[ifac] = -(boundary_tag_bit(BTAG_ALL) + | boundary_tag_bit(BTAG_REALLY_ALL) | boundary_tag_bit(marker)) exterior_grp = FacialAdjacencyGroup(igroup=0, ineighbor=None, @@ -240,12 +265,16 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh): def _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices, normals=None, no_normals_warn=True): - # FIXME : Fix docs """ 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. - Mesh must have co-dimension 0 or 1. + + :param fdrake_mesh: A :mod:`firedrake` instance of :class:`MeshGeometry` + + :param unflipped_group: A :class:`SimplexElementGroup` instance with + (potentially) some negatively oriented elements. + + :param vertices: The vertex coordinates as a numpy array of shape + *(ambient_dim, nvertices)* :param normals: _Only_ used if :param:`mesh` is a 1-surface embedded in 2-space. In this case, @@ -259,6 +288,10 @@ def _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices, :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*. + + :return: A numpy array, the *i*th element is > 0 if the *ith* element + is positively oriented, < 0 if negatively oriented. + Mesh must have co-dimension 0 or 1. """ # compute orientations tdim = fdrake_mesh.topological_dimension() @@ -309,11 +342,38 @@ def _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices, def import_firedrake_mesh(fdrake_mesh): # FIXME : docs """ - :return: A tuple (meshmode mesh, firedrake_orient). - firedrake_orient < 0 is True for any negatively - oriented firedrake cell (which was flipped by meshmode) - and False for any positively oriented firedrake cell - (whcih was not flipped by meshmode) + Create a :mod:`meshmode` :class:`Mesh` from a :mod:`firedrake` + :class:`MeshGeometry` with the same cells/elements, vertices, nodes, + mesh order, and facial adjacency. + + The vertex and node coordinates will be the same, as well + as the cell/element ordering. However, :mod:`firedrake` + does not require elements to be positively oriented, + so any negative elements are flipped + as in :func:`meshmode.processing.flip_simplex_element_group`. + + The flipped cells/elements are identified by the returned + *firedrake_orient* array + + :param fdrake_mesh: A :mod:`firedrake` :class:`MeshGeometry`. + This mesh **must** be in a space of ambient dimension + 1, 2, or 3 and have co-dimension of 0 or 1. + It must use a simplex as a reference element. + + Finally, its ``coordinates`` attribute must have a function + space whose *finat_element* associates a degree + of freedom with each vertex. In particular, + this means that the vertices of the mesh must have well-defined + coordinates. + For those unfamiliar with :mod:`firedrake`, you can + verify this by looking at the ``[0]`` entry of + ``fdrake_mesh.coordinates.function_space().finat_element.entity_dofs()``. + + :return: A tuple *(meshmode mesh, firedrake_orient)*. + ``firedrake_orient < 0`` is *True* for any negatively + oriented firedrake cell (which was flipped by meshmode) + and False for any positively oriented firedrake cell + (whcih was not flipped by meshmode). """ # Type validation from firedrake.mesh import MeshGeometry @@ -324,6 +384,7 @@ def import_firedrake_mesh(fdrake_mesh): assert fdrake_mesh.ufl_cell().is_simplex(), "Mesh must use simplex cells" gdim = fdrake_mesh.geometric_dimension() tdim = fdrake_mesh.topological_dimension() + assert gdim in [1, 2, 3], "Mesh must be in space of ambient dim 1, 2, or 3" assert gdim - tdim in [0, 1], "Mesh co-dimension must be 0 or 1" fdrake_mesh.init() diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py deleted file mode 100644 index da965cbe4148b7f264622b5a83b99d91d930bd11..0000000000000000000000000000000000000000 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ /dev/null @@ -1,461 +0,0 @@ -__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 deleted file mode 100644 index 3a5d7a2c1a68b3313c2f4d4c114d3f44dc5b5411..0000000000000000000000000000000000000000 --- a/meshmode/interop/firedrake/mesh_topology.py +++ /dev/null @@ -1,221 +0,0 @@ -__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 - -# }}}