From 11f3a2871af755f92481529ba5909de2f898ddd1 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 25 May 2020 11:12:46 -0500 Subject: [PATCH 01/33] Set up basic interface for object conversion --- doc/conf.py | 3 +- meshmode/interop/__init__.py | 57 ++++++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 1 deletion(-) create mode 100644 meshmode/interop/__init__.py diff --git a/doc/conf.py b/doc/conf.py index 86665362..8ff4e0bb 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/meshmode/interop/__init__.py b/meshmode/interop/__init__.py new file mode 100644 index 00000000..e3ca48fc --- /dev/null +++ b/meshmode/interop/__init__.py @@ -0,0 +1,57 @@ +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from abc import ABC + +__doc__ = """ +Development Interface +--------------------- +.. autoclass:: Interoperator +""" + + +class Interoperator(ABC): + """ + An object which handles the interoperation of an external + object with its closest meshmode representation + """ + + def __init__(self): + pass + + def as_meshmode(self): + pass + + def as_external(self): + pass + + def __hash__(self): + pass + + def __eq__(self, other): + pass + + def __neq__(self, other): + return not self.__eq__(other) + + +# vim: fdm=marker -- GitLab From bd5f322c94a687972b1eef02f81a6d69a477fc33 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 11:27:53 -0500 Subject: [PATCH 02/33] Update interface to be unidirectional --- meshmode/interop/__init__.py | 82 ++++++++++++++++++++++++++++++------ 1 file changed, 69 insertions(+), 13 deletions(-) diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py index e3ca48fc..c08a36e1 100644 --- a/meshmode/interop/__init__.py +++ b/meshmode/interop/__init__.py @@ -25,33 +25,89 @@ from abc import ABC __doc__ = """ Development Interface --------------------- -.. autoclass:: Interoperator +.. autoclass:: ExternalTransporter +.. autoclass:: ExternalImporter +.. autoclass:: ExternalExporter """ +# {{{ Generic, most abstract class for transporting meshmode <-> external -class Interoperator(ABC): - """ - An object which handles the interoperation of an external - object with its closest meshmode representation +class ExternalTransporter(ABC): """ + .. attribute:: from_data + + The object which needs to be transported either to meshmode or + from meshmode - def __init__(self): - pass + .. attribute:: to_data - def as_meshmode(self): - pass + The "transported" object, i.e. the - def as_external(self): - pass + This attribute does not exist at instantiation time. + If exporting (resp. importing) from meshmode + then we are using an :class:`ExternalExporter` + (resp. :class:`ExternalImporter`) instance. :attr:`to_data` is + computed with a call to :fun:`ExternalExporter.export` + (resp. :fun:`ExternalImporter.import`). + + :raises ValueError: if :attr:`to_data` is accessed before creation. + :raises NotImplementedError: if :meth:`validate_to_data` is called + without an implementation. + """ + def __init__(self, from_data): + self.from_data = from_data + + def validate_to_data(self): + """ + Validate :attr:`to_data` + + :return: *True* if :attr:`to_data` has been computed and is valid + and *False* otherwise + """ + raise NotImplementedError("*validate_to_data* method not implemented " \ + "for object of type %s" % type(self)) def __hash__(self): - pass + return hash((type(self), self.from_data)) def __eq__(self, other): - pass + return isinstance(other, type(self)) and \ + isinstance(self, type(other)) and \ + self.from_data == other.from_data def __neq__(self, other): return not self.__eq__(other) + def __getattr__(self, attr): + if attr != 'to_data': + return super(ExternalTransporter, self).__getattr__(attr) + raise ValueError("Attribute *to_data* has not yet been computed. " \ + "An object of class *ExternalExporter* (resp. " \ + "*ExternalImporter*) must call *export()* " \ + "(resp. *import()*) to compute attribute *to_data*") + +# }}} + + +# {{{ Define specific classes for meshmode -> external and meshmode <- external + +class ExternalExporter(ExternalTransporter): + def export(self): + """ + Compute :attr:`to_data` from :attr:`from_data` + """ + raise NotImplementedError("*export* method not implemented " \ + "for type %s" % type(self)) + +class ExternalImporter(ExternalTransporter): + def import(self): + """ + Compute :attr:`to_data` from :attr:`from_data` + """ + raise NotImplementedError("*import* method not implemented " \ + "for type %s" % type(self)) + +# }}} + # vim: fdm=marker -- GitLab From 70657ca9e56c3fea77d3efd3b336797e8cd91ad2 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 11:41:52 -0500 Subject: [PATCH 03/33] Added an importer for FIAT's simplex cell --- meshmode/interop/fiat/simplex_cell.py | 123 ++++++++++++++++++++++++++ 1 file changed, 123 insertions(+) create mode 100644 meshmode/interop/fiat/simplex_cell.py diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py new file mode 100644 index 00000000..51c61a04 --- /dev/null +++ b/meshmode/interop/fiat/simplex_cell.py @@ -0,0 +1,123 @@ +from warnings import warn # noqa +import numpy as np +import numpy.linalg as la + +from meshmode.mesh.interop import ExternalImporter + + +__doc__ = """ +Interoperators +-------------- +.. autoclass:: SimplexCellAnalog + :members: +""" + +# {{{ 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(ExternalImporter): + """ + Importer for a :mod:`FIAT` simplex cell. + There is no data created from the :attr:`from_data`. Instead, + the data is simply used to obtain FIAT's reference + nodes according to :mod:`modepy`'s reference coordinates + using :meth:`make_modepy_points`. + + .. attribute:: from_data + + An instance of :class:`fiat.FIAT.reference_element.Simplex`. + + .. attribute:: to_data + + :raises ValueError: If accessed. + """ + 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(SimplexCellInteroperator, 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 make_modepy_points(self, dim, entity_id, order): + """ + Constructs a lattice of points on the *entity_id*th facet + of dimension *dim*. + + Args are exactly as in + :meth:`fiat.FIAT.reference_element.Cell.make_points` + (see `FIAT docs `_), + but the unit nodes are (affinely) mapped to :mod:`modepy` + `unit coordinates `_. + + :arg dim: Dimension of the facet we are constructing points on. + :arg entity_id: identifier to determine which facet of + dimension *dim* to construct the points on. + :arg order: Number of points to include in each direction. + + :return: an *np.array* of shape *(dim, npoints)* holding the + coordinates of each of the ver + """ + points = self.from_data.make_points(dim, entity_id, order) + if not points: + return points + points = np.array(points) + # Points is (npoints, dim) so have to transpose + return (np.matmul(self._mat, points.T) + self._shift[:, np.newaxis]).T + +# }}} -- GitLab From 4ff6048f106e80514ecc9ffabc55258579e0a637 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 11:42:44 -0500 Subject: [PATCH 04/33] Added an __init__.py file for FIAT directory --- meshmode/interop/fiat/__init__.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 meshmode/interop/fiat/__init__.py diff --git a/meshmode/interop/fiat/__init__.py b/meshmode/interop/fiat/__init__.py new file mode 100644 index 00000000..f8797739 --- /dev/null +++ b/meshmode/interop/fiat/__init__.py @@ -0,0 +1,24 @@ +from future import division + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + -- GitLab From 07d0f416ecf5dae25938fbd87affdd6147ecc143 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 11:44:37 -0500 Subject: [PATCH 05/33] Fixed docstring --- meshmode/interop/fiat/simplex_cell.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py index 51c61a04..b582094a 100644 --- a/meshmode/interop/fiat/simplex_cell.py +++ b/meshmode/interop/fiat/simplex_cell.py @@ -6,10 +6,7 @@ from meshmode.mesh.interop import ExternalImporter __doc__ = """ -Interoperators --------------- -.. autoclass:: SimplexCellAnalog - :members: +.. autoclass:: FIATSimplexCellImporter """ # {{{ Compute an affine mapping from given input/outputs -- GitLab From 7491da5e2b930f06d7f005571ae9a1b3ca50f053 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 11:58:19 -0500 Subject: [PATCH 06/33] Some flake8 fixes --- meshmode/interop/__init__.py | 36 ++++++++++++++++++++----------- meshmode/interop/fiat/__init__.py | 3 --- 2 files changed, 23 insertions(+), 16 deletions(-) diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py index c08a36e1..91ddfd76 100644 --- a/meshmode/interop/__init__.py +++ b/meshmode/interop/__init__.py @@ -30,6 +30,7 @@ Development Interface .. autoclass:: ExternalExporter """ + # {{{ Generic, most abstract class for transporting meshmode <-> external class ExternalTransporter(ABC): @@ -45,10 +46,10 @@ class ExternalTransporter(ABC): This attribute does not exist at instantiation time. If exporting (resp. importing) from meshmode - then we are using an :class:`ExternalExporter` - (resp. :class:`ExternalImporter`) instance. :attr:`to_data` is - computed with a call to :fun:`ExternalExporter.export` - (resp. :fun:`ExternalImporter.import`). + then we are using an :class:`ExternalExporter` + (resp. :class:`ExternalImporter`) instance. :attr:`to_data` is + computed with a call to :fun:`ExternalExporter.export_data` + (resp. :fun:`ExternalImporter.import_data`). :raises ValueError: if :attr:`to_data` is accessed before creation. :raises NotImplementedError: if :meth:`validate_to_data` is called @@ -64,7 +65,7 @@ class ExternalTransporter(ABC): :return: *True* if :attr:`to_data` has been computed and is valid and *False* otherwise """ - raise NotImplementedError("*validate_to_data* method not implemented " \ + raise NotImplementedError("*validate_to_data* method not implemented " "for object of type %s" % type(self)) def __hash__(self): @@ -81,10 +82,10 @@ class ExternalTransporter(ABC): def __getattr__(self, attr): if attr != 'to_data': return super(ExternalTransporter, self).__getattr__(attr) - raise ValueError("Attribute *to_data* has not yet been computed. " \ - "An object of class *ExternalExporter* (resp. " \ - "*ExternalImporter*) must call *export()* " \ - "(resp. *import()*) to compute attribute *to_data*") + raise ValueError("Attribute *to_data* has not yet been computed. " + "An object of class *ExternalExporter* (resp. " + "*ExternalImporter*) must call *export_data()* " + "(resp. *import_data()*) to compute attribute *to_data*") # }}} @@ -92,19 +93,28 @@ class ExternalTransporter(ABC): # {{{ Define specific classes for meshmode -> external and meshmode <- external class ExternalExporter(ExternalTransporter): - def export(self): + """ + Subclass of :class:`ExternalTransporter` for meshmode -> external + data transfer + """ + def export_data(self): """ Compute :attr:`to_data` from :attr:`from_data` """ - raise NotImplementedError("*export* method not implemented " \ + raise NotImplementedError("*export_data* method not implemented " "for type %s" % type(self)) + class ExternalImporter(ExternalTransporter): - def import(self): + """ + Subclass of :class:`ExternalTransporter` for external -> meshmode + data transfer + """ + def import_data(self): """ Compute :attr:`to_data` from :attr:`from_data` """ - raise NotImplementedError("*import* method not implemented " \ + raise NotImplementedError("*import_data* method not implemented " "for type %s" % type(self)) # }}} diff --git a/meshmode/interop/fiat/__init__.py b/meshmode/interop/fiat/__init__.py index f8797739..908c7cd2 100644 --- a/meshmode/interop/fiat/__init__.py +++ b/meshmode/interop/fiat/__init__.py @@ -1,5 +1,3 @@ -from future import division - __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" __license__ = """ @@ -21,4 +19,3 @@ 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. """ - -- GitLab From fc9baf080fe6bc838c292e9427a6b4f56df72153 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 12:26:22 -0500 Subject: [PATCH 07/33] Some more minor flake8 changes --- meshmode/interop/fiat/__init__.py | 4 ++++ meshmode/interop/fiat/simplex_cell.py | 22 ++++++++++++++-------- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/meshmode/interop/fiat/__init__.py b/meshmode/interop/fiat/__init__.py index 908c7cd2..3536a1ea 100644 --- a/meshmode/interop/fiat/__init__.py +++ b/meshmode/interop/fiat/__init__.py @@ -19,3 +19,7 @@ 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 index b582094a..f2b8366b 100644 --- a/meshmode/interop/fiat/simplex_cell.py +++ b/meshmode/interop/fiat/simplex_cell.py @@ -1,14 +1,14 @@ -from warnings import warn # noqa import numpy as np import numpy.linalg as la -from meshmode.mesh.interop import ExternalImporter +from meshmode.interop import ExternalImporter __doc__ = """ .. autoclass:: FIATSimplexCellImporter """ + # {{{ Compute an affine mapping from given input/outputs def get_affine_mapping(reference_vects, vects): @@ -57,17 +57,23 @@ def get_affine_mapping(reference_vects, vects): class FIATSimplexCellImporter(ExternalImporter): """ Importer for a :mod:`FIAT` simplex cell. - There is no data created from the :attr:`from_data`. Instead, - the data is simply used to obtain FIAT's reference + + There is no obvious meshmode counterpart for :attr:`from_data`. + The data is simply used to obtain FIAT's reference nodes according to :mod:`modepy`'s reference coordinates - using :meth:`make_modepy_points`. + using :meth:`make_points`. + + In particular, methods + :meth:`import_data` and :meth:`validate_to_data` + are *NOT* implemented + and there is no :attr:`to_data` to be computed. .. attribute:: from_data An instance of :class:`fiat.FIAT.reference_element.Simplex`. .. attribute:: to_data - + :raises ValueError: If accessed. """ def __init__(self, cell): @@ -78,7 +84,7 @@ class FIATSimplexCellImporter(ExternalImporter): from FIAT.reference_element import Simplex assert isinstance(cell, Simplex) - super(SimplexCellInteroperator, self).__init__(cell) + super(FIATSimplexCellImporter, self).__init__(cell) # Stored as (dim, nunit_vertices) from modepy.tools import unit_vertices @@ -91,7 +97,7 @@ class FIATSimplexCellImporter(ExternalImporter): self._mat, self._shift = get_affine_mapping(reference_vertices, self._unit_vertices) - def make_modepy_points(self, dim, entity_id, order): + def make_points(self, dim, entity_id, order): """ Constructs a lattice of points on the *entity_id*th facet of dimension *dim*. -- GitLab From d67546125459dc7a10561610d58bbfe2651f4052 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 12:26:43 -0500 Subject: [PATCH 08/33] Added an importer for FInAT Lagrange elements --- meshmode/interop/FInAT/__init__.py | 23 +++ meshmode/interop/FInAT/lagrange_element.py | 196 +++++++++++++++++++++ 2 files changed, 219 insertions(+) create mode 100644 meshmode/interop/FInAT/__init__.py create mode 100644 meshmode/interop/FInAT/lagrange_element.py diff --git a/meshmode/interop/FInAT/__init__.py b/meshmode/interop/FInAT/__init__.py new file mode 100644 index 00000000..f53a2135 --- /dev/null +++ b/meshmode/interop/FInAT/__init__.py @@ -0,0 +1,23 @@ +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +__all__ = ['FinatLagrangeElementImporter'] diff --git a/meshmode/interop/FInAT/lagrange_element.py b/meshmode/interop/FInAT/lagrange_element.py new file mode 100644 index 00000000..c74b2e53 --- /dev/null +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -0,0 +1,196 @@ +import numpy as np +import numpy.linalg as la +import six + +from meshmode.interop import ExternalImporter +from meshmode.interop.fiat import FIATSimplexCellImporter + + +__doc__ = """ +.. autoclass:: FinatLagrangeElementImporter + :members: +""" + + +class FinatLagrangeElementImporter(ExternalImporter): + """ + An importer for a FInAT element, usually instantiated from + ``some_instantiated_firedrake_function_space.finat_element`` + + This importer does not have an obvious meshmode counterpart, + so is instead used for its method. + In particular, methods :meth:`import_data` + and :meth:`validate_to_data` are *NOT* implemented + and there is no :attr:`to_data` to be computed. + """ + def __init__(self, finat_element): + """ + :param finat_element: A FInAT element of type + :class:`finat.fiat_elements.Lagrange` or + :class:`finat.fiat_elements.DiscontinuousLagrange` + 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` + + :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 + if not isinstance(finat_element, (Lagrange, DiscontinuousLagrange)): + raise TypeError(":param:`finat_element` must be of type" + " `finat.fiat_elements.Lagrange` or" + " `finat.fiat_elements.DiscontinuousLagrange`", + " not 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: + # {{{ Compute unit nodes + node_nr_to_coords = {} + unit_vertex_indices = [] + + # Get unit nodes + for dim, element_nrs in six.iteritems( + self.analog().entity_support_dofs()): + for element_nr, node_list in six.iteritems(element_nrs): + # Get the nodes on the element (in meshmode reference coords) + pts_on_element = self.cell_importer.make_points( + dim, element_nr, self.analog().degree) + # Record any new nodes + i = 0 + for node_nr in node_list: + if node_nr not in node_nr_to_coords: + node_nr_to_coords[node_nr] = pts_on_element[i] + i += 1 + # If is a vertex, store the index + if dim == 0: + unit_vertex_indices.append(node_nr) + + # store vertex indices + self._unit_vertex_indices = np.array(sorted(unit_vertex_indices)) + + # Convert unit_nodes to array, then change to (dim, nunit_nodes) + # from (nunit_nodes, dim) + unit_nodes = np.array([node_nr_to_coords[i] for i in + range(len(node_nr_to_coords))]) + self._unit_nodes = unit_nodes.T.copy() + + # }}} + + def dim(self): + """ + :return: The dimension of the FInAT element's cell + """ + return self.cell_importer.from_data.get_dimension() + + def unit_vertex_indices(self): + """ + :return: An array of shape *(dim+1,)* of indices + so that *self.unit_nodes()[self.unit_vertex_indices()]* + are the vertices of the reference element. + """ + 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.analog().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()) -- GitLab From ec0715190f029d07e9c1f88d55377c94720a0ed1 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Thu, 28 May 2020 15:28:21 -0500 Subject: [PATCH 09/33] Added a mesh topology import handler --- meshmode/interop/firedrake/mesh_topology.py | 200 ++++++++++++++++++++ 1 file changed, 200 insertions(+) create mode 100644 meshmode/interop/firedrake/mesh_topology.py diff --git a/meshmode/interop/firedrake/mesh_topology.py b/meshmode/interop/firedrake/mesh_topology.py new file mode 100644 index 00000000..d4ec055e --- /dev/null +++ b/meshmode/interop/firedrake/mesh_topology.py @@ -0,0 +1,200 @@ +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 + + # variable names follow a dmplex naming convention for + # cell Start/end and vertex Start/end. + cStart, cEnd = plex.getHeightStratum(0) + vStart, vEnd = 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(cStart, cEnd, 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(cStart, cEnd) + # Otherwise, just the ones you're using + else: + isin = np.isin(to_fd_id, self.icell_to_fd) + range_ = np.arange(cStart, cEnd, 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 vStart <= vert_id < vEnd 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 cStart <= other_cell_id < cEnd: + cells.append(to_fd_id[other_cell_id - cStart]) + + # 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 + +# }}} -- GitLab From 4dd0d6adf1a139a28a3ec4dde156a69eef1968c7 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 11:15:54 -0500 Subject: [PATCH 10/33] Added coordinateless functions implementation --- .../coordinateless_functions_impl.py | 195 ++++++++++++++++++ 1 file changed, 195 insertions(+) create mode 100644 meshmode/interop/firedrake/coordinateless_functions_impl.py diff --git a/meshmode/interop/firedrake/coordinateless_functions_impl.py b/meshmode/interop/firedrake/coordinateless_functions_impl.py new file mode 100644 index 00000000..0a4e7bc4 --- /dev/null +++ b/meshmode/interop/firedrake/coordinateless_functions_impl.py @@ -0,0 +1,195 @@ +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. + """ + # {{{ 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() == 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.") + + # }}} + + # We want to ignore any geometry and then finish initialization + function_space = function_space.topological + 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()`` + """ + # {{{ 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.") + + function = function.topological + function_space_importer = function_space_importer.topological_importer + + 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 + + +# }}} -- GitLab From aa9bd6d93a5e74839f09eec554b9035c0b43f45c Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 11:16:04 -0500 Subject: [PATCH 11/33] Forgot to add __init__.py --- meshmode/interop/firedrake/__init__.py | 33 ++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 meshmode/interop/firedrake/__init__.py diff --git a/meshmode/interop/firedrake/__init__.py b/meshmode/interop/firedrake/__init__.py new file mode 100644 index 00000000..67712b23 --- /dev/null +++ b/meshmode/interop/firedrake/__init__.py @@ -0,0 +1,33 @@ +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +__doc__ = """ +.. autoclass:: FiredrakeConnection +""" + +from meshmode.discretization.connection import DiscretizationConnection + + +# TODO +class FiredrakeConnection(DiscretizationConnection): + def __init__(self, cl_ctx, fdrake_function_space): + pass -- GitLab From 1c3b93bf70be68a579e8243b0bba68d345380860 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 11:16:46 -0500 Subject: [PATCH 12/33] Changed name of underlying structure from transporter to handler --- meshmode/interop/FInAT/__init__.py | 2 + meshmode/interop/FInAT/lagrange_element.py | 18 ++--- meshmode/interop/__init__.py | 80 ++++++---------------- meshmode/interop/fiat/simplex_cell.py | 24 ++----- 4 files changed, 34 insertions(+), 90 deletions(-) diff --git a/meshmode/interop/FInAT/__init__.py b/meshmode/interop/FInAT/__init__.py index f53a2135..3708afa4 100644 --- a/meshmode/interop/FInAT/__init__.py +++ b/meshmode/interop/FInAT/__init__.py @@ -20,4 +20,6 @@ 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 index c74b2e53..ce71d6b7 100644 --- a/meshmode/interop/FInAT/lagrange_element.py +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -2,7 +2,7 @@ import numpy as np import numpy.linalg as la import six -from meshmode.interop import ExternalImporter +from meshmode.interop import ExternalImportHandler from meshmode.interop.fiat import FIATSimplexCellImporter @@ -12,16 +12,10 @@ __doc__ = """ """ -class FinatLagrangeElementImporter(ExternalImporter): +class FinatLagrangeElementImporter(ExternalImportHandler): """ An importer for a FInAT element, usually instantiated from ``some_instantiated_firedrake_function_space.finat_element`` - - This importer does not have an obvious meshmode counterpart, - so is instead used for its method. - In particular, methods :meth:`import_data` - and :meth:`validate_to_data` are *NOT* implemented - and there is no :attr:`to_data` to be computed. """ def __init__(self, finat_element): """ @@ -76,11 +70,11 @@ class FinatLagrangeElementImporter(ExternalImporter): # Get unit nodes for dim, element_nrs in six.iteritems( - self.analog().entity_support_dofs()): + self.data.entity_support_dofs()): for element_nr, node_list in six.iteritems(element_nrs): # Get the nodes on the element (in meshmode reference coords) pts_on_element = self.cell_importer.make_points( - dim, element_nr, self.analog().degree) + dim, element_nr, self.data.degree) # Record any new nodes i = 0 for node_nr in node_list: @@ -106,7 +100,7 @@ class FinatLagrangeElementImporter(ExternalImporter): """ :return: The dimension of the FInAT element's cell """ - return self.cell_importer.from_data.get_dimension() + return self.cell_importer.data.get_dimension() def unit_vertex_indices(self): """ @@ -162,7 +156,7 @@ class FinatLagrangeElementImporter(ExternalImporter): from modepy import resampling_matrix, simplex_best_available_basis flip_matrix = resampling_matrix( - simplex_best_available_basis(self.dim(), self.analog().degree), + simplex_best_available_basis(self.dim(), self.data.degree), flipped_unit_nodes, self.unit_nodes()) flip_matrix[np.abs(flip_matrix) < 1e-15] = 0 diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py index 91ddfd76..bd6d9257 100644 --- a/meshmode/interop/__init__.py +++ b/meshmode/interop/__init__.py @@ -25,97 +25,59 @@ from abc import ABC __doc__ = """ Development Interface --------------------- -.. autoclass:: ExternalTransporter -.. autoclass:: ExternalImporter -.. autoclass:: ExternalExporter +.. autoclass:: ExternalDataHandler +.. autoclass:: ExternalExportHandler +.. autoclass:: ExternalImportHandler """ # {{{ Generic, most abstract class for transporting meshmode <-> external -class ExternalTransporter(ABC): +class ExternalDataHandler(ABC): """ - .. attribute:: from_data + 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. - The object which needs to be transported either to meshmode or - from meshmode + .. attribute:: data - .. attribute:: to_data - - The "transported" object, i.e. the - - This attribute does not exist at instantiation time. - If exporting (resp. importing) from meshmode - then we are using an :class:`ExternalExporter` - (resp. :class:`ExternalImporter`) instance. :attr:`to_data` is - computed with a call to :fun:`ExternalExporter.export_data` - (resp. :fun:`ExternalImporter.import_data`). - - :raises ValueError: if :attr:`to_data` is accessed before creation. - :raises NotImplementedError: if :meth:`validate_to_data` is called - without an implementation. + The object which needs to be interfaced either into meshmode or + out of meshmode. + Should not be modified after creation. """ - def __init__(self, from_data): - self.from_data = from_data - - def validate_to_data(self): - """ - Validate :attr:`to_data` - - :return: *True* if :attr:`to_data` has been computed and is valid - and *False* otherwise - """ - raise NotImplementedError("*validate_to_data* method not implemented " - "for object of type %s" % type(self)) + def __init__(self, data): + self.data = data def __hash__(self): - return hash((type(self), self.from_data)) + return hash((type(self), self.data)) def __eq__(self, other): return isinstance(other, type(self)) and \ isinstance(self, type(other)) and \ - self.from_data == other.from_data + self.data == other.data def __neq__(self, other): return not self.__eq__(other) - def __getattr__(self, attr): - if attr != 'to_data': - return super(ExternalTransporter, self).__getattr__(attr) - raise ValueError("Attribute *to_data* has not yet been computed. " - "An object of class *ExternalExporter* (resp. " - "*ExternalImporter*) must call *export_data()* " - "(resp. *import_data()*) to compute attribute *to_data*") - # }}} # {{{ Define specific classes for meshmode -> external and meshmode <- external -class ExternalExporter(ExternalTransporter): +class ExternalExportHandler(ExternalDataHandler): """ - Subclass of :class:`ExternalTransporter` for meshmode -> external + Subclass of :class:`ExternalDataHandler` for meshmode -> external data transfer """ - def export_data(self): - """ - Compute :attr:`to_data` from :attr:`from_data` - """ - raise NotImplementedError("*export_data* method not implemented " - "for type %s" % type(self)) + pass -class ExternalImporter(ExternalTransporter): +class ExternalImportHandler(ExternalDataHandler): """ - Subclass of :class:`ExternalTransporter` for external -> meshmode + Subclass of :class:`ExternalDataHandler` for external -> meshmode data transfer """ - def import_data(self): - """ - Compute :attr:`to_data` from :attr:`from_data` - """ - raise NotImplementedError("*import_data* method not implemented " - "for type %s" % type(self)) + pass # }}} diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py index f2b8366b..bdcabd83 100644 --- a/meshmode/interop/fiat/simplex_cell.py +++ b/meshmode/interop/fiat/simplex_cell.py @@ -1,7 +1,7 @@ import numpy as np import numpy.linalg as la -from meshmode.interop import ExternalImporter +from meshmode.interop import ExternalImportHandler __doc__ = """ @@ -54,27 +54,13 @@ def get_affine_mapping(reference_vects, vects): # {{{ Interoperator for FIAT's reference_element.Simplex -class FIATSimplexCellImporter(ExternalImporter): +class FIATSimplexCellImporter(ExternalImportHandler): """ - Importer for a :mod:`FIAT` simplex cell. + Import handler for a :mod:`FIAT` simplex cell. - There is no obvious meshmode counterpart for :attr:`from_data`. - The data is simply used to obtain FIAT's reference - nodes according to :mod:`modepy`'s reference coordinates - using :meth:`make_points`. - - In particular, methods - :meth:`import_data` and :meth:`validate_to_data` - are *NOT* implemented - and there is no :attr:`to_data` to be computed. - - .. attribute:: from_data + .. attribute:: data An instance of :class:`fiat.FIAT.reference_element.Simplex`. - - .. attribute:: to_data - - :raises ValueError: If accessed. """ def __init__(self, cell): """ @@ -116,7 +102,7 @@ class FIATSimplexCellImporter(ExternalImporter): :return: an *np.array* of shape *(dim, npoints)* holding the coordinates of each of the ver """ - points = self.from_data.make_points(dim, entity_id, order) + points = self.data.make_points(dim, entity_id, order) if not points: return points points = np.array(points) -- GitLab From 30f04b46261b34d7575cdfe02addfa815b182943 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 11:51:08 -0500 Subject: [PATCH 13/33] Added an importer for mesh geometries --- meshmode/interop/firedrake/mesh_geometry.py | 408 ++++++++++++++++++++ 1 file changed, 408 insertions(+) create mode 100644 meshmode/interop/firedrake/mesh_geometry.py diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py new file mode 100644 index 00000000..313f7c3e --- /dev/null +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -0,0 +1,408 @@ +from warnings import warn # noqa +from collections import defaultdict +import numpy as np + +from meshmode.interop import ExternalImportHandler +from meshmode.interop.firedrake.coordinateless_functions_impl 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 != 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.functionspace import \ + FiredrakeWithGeometryImporter + from firedrake import Function + + coordinates_fs = self.data.coordinates.function_space() + coordinates_fs_importer = \ + self._coordinates_importer.function_space_importer() + + V_importer = FiredrakeWithGeometryImporter(cl_ctx, + coordinates_fs, + coordinates_fs_importer, + self) + f = Function(V_importer.data, val=self._coordinates_a.data) + self._coordinates_function_importer = \ + FiredrakeFunctionImporter(f, V_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 + + # 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 group(self): + """ + Return an instance of :class:meshmode.mesh.SimplexElementGroup` + corresponding to the mesh :attr:`data` + """ + if self._group is None: + from meshmode.mesh import SimplexElementGroup + from meshmode.mesh.processing import flip_simplex_element_group + + fspace_importer = self.coordinates_importer.function_space_importer() + finat_element_importer = fspace_importer.finat_element_importer + + # IMPORTANT that set :attr:`_group` because + # :meth:`orientations` may call :meth:`group` + self._group = SimplexElementGroup( + finat_element_importer.data.degree, + self.vertex_indices(), + self.nodes(), + dim=self.cell_dimension(), + unit_nodes=finat_element_importer.unit_nodes()) + + self._group = flip_simplex_element_group(self.vertices(), self._group, + self.orientations() < 0) + + 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.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 -- GitLab From a2ab04f23ff534bd85e72e0ac8b0c08f26e2da7e Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 12:12:14 -0500 Subject: [PATCH 14/33] Added function space shared data --- ...ns_impl.py => function_space_coordless.py} | 0 .../firedrake/function_space_shared_data.py | 243 ++++++++++++++++++ meshmode/interop/firedrake/mesh_geometry.py | 2 +- 3 files changed, 244 insertions(+), 1 deletion(-) rename meshmode/interop/firedrake/{coordinateless_functions_impl.py => function_space_coordless.py} (100%) create mode 100644 meshmode/interop/firedrake/function_space_shared_data.py diff --git a/meshmode/interop/firedrake/coordinateless_functions_impl.py b/meshmode/interop/firedrake/function_space_coordless.py similarity index 100% rename from meshmode/interop/firedrake/coordinateless_functions_impl.py rename to meshmode/interop/firedrake/function_space_coordless.py 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 00000000..e5032b9b --- /dev/null +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -0,0 +1,243 @@ +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_analog, key, *args, **kwargs): + """ + Exactly :func:`firedrake.functionspacedata.cached`, but + caching on mesh Geometry instead + + :param f: The function to cache. + :param mesh_analog: The mesh_analog 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_analog, "_shared_data_cache") + cache = mesh_analog._shared_data_cache[f.__name__] + try: + return cache[key] + except KeyError: + result = f(mesh_analog, key, *args, **kwargs) + cache[key] = result + return result + + +def reorder_nodes(orient, nodes, flip_matrix, unflip=False): + """ + :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 + + flips :param:`nodes` + """ + # 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_anlog, 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:`fd2mm.function_space.FunctionSpaceAnalog.reorder_nodes` + """ + 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_analog, degree): + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + return InterpolatoryQuadratureSimplexGroupFactory(degree) + + +@cached +def get_discretization(mesh_analog, key): + finat_element_importer, cl_ctx = key + assert isinstance(finat_element_importer, FinatLagrangeElementImporter) + + degree = finat_element_importer.data.degree + discretization = Discretization(cl_ctx, + mesh_analog.meshmode_mesh(), + get_factory(mesh_analog, degree)) + + return discretization + + +class FunctionSpaceDataAnalog(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_analog` and :param:`finat_element_analog` + are used for analog-checking + """ + if mesh_importer.topological_a == mesh_importer: + raise TypeError(":param:`mesh_importer` is a " + "FiredrakeMeshTopologyImporter, but " + " must be a FiredrakeMeshGeometryImporter") + importer = (mesh_importer.importer(), finat_element_importer.importer()) + super(FunctionSpaceDataAnalog, self).__init__(importer) + + self._fspace_data = FunctionSpaceData(mesh_importer.importer(), + finat_element_importer.importer()) + + 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 index 313f7c3e..23d803a7 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -3,7 +3,7 @@ from collections import defaultdict import numpy as np from meshmode.interop import ExternalImportHandler -from meshmode.interop.firedrake.coordinateless_functions_impl import \ +from meshmode.interop.firedrake.function_space_coordless import \ FiredrakeCoordinatelessFunctionImporter -- GitLab From eedefb9a48261a127f0f33d69cf503dd24aa58fd Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 13:02:59 -0500 Subject: [PATCH 15/33] Added a function space importer --- meshmode/interop/firedrake/function_space.py | 184 ++++++++++++++++++ .../firedrake/function_space_shared_data.py | 2 +- meshmode/interop/firedrake/mesh_geometry.py | 2 +- 3 files changed, 186 insertions(+), 2 deletions(-) create mode 100644 meshmode/interop/firedrake/function_space.py diff --git a/meshmode/interop/firedrake/function_space.py b/meshmode/interop/firedrake/function_space.py new file mode 100644 index 00000000..ca4bfd63 --- /dev/null +++ b/meshmode/interop/firedrake/function_space.py @@ -0,0 +1,184 @@ +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 Analog + 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.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_a, 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_a.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_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index e5032b9b..9644f0f1 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -167,7 +167,7 @@ def get_discretization(mesh_analog, key): return discretization -class FunctionSpaceDataAnalog(ExternalImportHandler): +class FiredrakeFunctionSpaceDataImporter(ExternalImportHandler): """ This is not *quite* the usual thought of a :class:`ExternalImportHandler`. diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py index 23d803a7..cd5c1d43 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -103,7 +103,7 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): """ from meshmode.interop.firedrake.function import \ FiredrakeFunctionImporter - from meshmode.interop.firedrake.functionspace import \ + from meshmode.interop.firedrake.function_space import \ FiredrakeWithGeometryImporter from firedrake import Function -- GitLab From e8f56cb89269dfca74127d975b7fda61dac69af4 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 13:17:19 -0500 Subject: [PATCH 16/33] Added a function importer --- meshmode/interop/firedrake/function.py | 96 ++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 meshmode/interop/firedrake/function.py diff --git a/meshmode/interop/firedrake/function.py b/meshmode/interop/firedrake/function.py new file mode 100644 index 00000000..ce601ec7 --- /dev/null +++ b/meshmode/interop/firedrake/function.py @@ -0,0 +1,96 @@ +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)[:] -- GitLab From 75e6337a83a3562d0b2540c98eb041d123ea4ffe Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 13:21:05 -0500 Subject: [PATCH 17/33] Some leftover changes from fd2mm --- meshmode/interop/firedrake/function_space.py | 2 +- .../firedrake/function_space_shared_data.py | 26 +++++++++---------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/meshmode/interop/firedrake/function_space.py b/meshmode/interop/firedrake/function_space.py index ca4bfd63..95f4e297 100644 --- a/meshmode/interop/firedrake/function_space.py +++ b/meshmode/interop/firedrake/function_space.py @@ -52,7 +52,7 @@ class FiredrakeWithGeometryImporter(ExternalImportHandler): # }}} - # Initialize as Analog + # Initialize as importer super(FiredrakeWithGeometryImporter, self).__init__(function_space) self._topology_importer = function_space_importer self._mesh_importer = mesh_importer diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index 9644f0f1..c58bbcd1 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -10,23 +10,23 @@ from meshmode.interop.FInAT.lagrange_element import \ @decorator -def cached(f, mesh_analog, key, *args, **kwargs): +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_analog: The mesh_analog to cache on (should have a + :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_analog, "_shared_data_cache") - cache = mesh_analog._shared_data_cache[f.__name__] + assert hasattr(mesh_importer, "_shared_data_cache") + cache = mesh_importer._shared_data_cache[f.__name__] try: return cache[key] except KeyError: - result = f(mesh_analog, key, *args, **kwargs) + result = f(mesh_importer, key, *args, **kwargs) cache[key] = result return result @@ -93,7 +93,7 @@ def reordering_array(mesh_importer, key, fspace_data): :class:`firedrake.functionspacedata.FunctionSpaceData` :returns: an *np.array* that can reorder the data by composition, - see :meth:`fd2mm.function_space.FunctionSpaceAnalog.reorder_nodes` + see :meth:`reorder_nodes` in class :class:`FiredrakeWithGeometryImporter` """ finat_element_importer, firedrake_to_meshmode = key assert isinstance(finat_element_importer, FinatLagrangeElementImporter) @@ -148,21 +148,21 @@ def reordering_array(mesh_importer, key, fspace_data): @cached -def get_factory(mesh_analog, degree): +def get_factory(mesh_importer, degree): from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory return InterpolatoryQuadratureSimplexGroupFactory(degree) @cached -def get_discretization(mesh_analog, key): +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_analog.meshmode_mesh(), - get_factory(mesh_analog, degree)) + mesh_importer.meshmode_mesh(), + get_factory(mesh_importer, degree)) return discretization @@ -191,15 +191,15 @@ class FiredrakeFunctionSpaceDataImporter(ExternalImportHandler): # FIXME: Give two finat elts def __init__(self, cl_ctx, mesh_importer, finat_element_importer): """ - Note that :param:`mesh_analog` and :param:`finat_element_analog` - are used for analog-checking + Note that :param:`mesh_importer` and :param:`finat_element_importer` + are used for checking """ if mesh_importer.topological_a == mesh_importer: raise TypeError(":param:`mesh_importer` is a " "FiredrakeMeshTopologyImporter, but " " must be a FiredrakeMeshGeometryImporter") importer = (mesh_importer.importer(), finat_element_importer.importer()) - super(FunctionSpaceDataAnalog, self).__init__(importer) + super(FiredrakeFunctionSpaceDataImporter, self).__init__(importer) self._fspace_data = FunctionSpaceData(mesh_importer.importer(), finat_element_importer.importer()) -- GitLab From 5a213f3340a0e2eafd1e174600474b48b99f16ed Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 15:31:56 -0500 Subject: [PATCH 18/33] Finished FromFiredrakeConnection --- meshmode/interop/firedrake/__init__.py | 231 ++++++++++++++++++++++++- 1 file changed, 226 insertions(+), 5 deletions(-) diff --git a/meshmode/interop/firedrake/__init__.py b/meshmode/interop/firedrake/__init__.py index 67712b23..11902b75 100644 --- a/meshmode/interop/firedrake/__init__.py +++ b/meshmode/interop/firedrake/__init__.py @@ -21,13 +21,234 @@ THE SOFTWARE. """ __doc__ = """ -.. autoclass:: FiredrakeConnection +.. autoclass:: FromFiredrakeConnection +.. autofunction:: import_firedrake_mesh +.. autofunction:: import_firedrake_function_space """ -from meshmode.discretization.connection import DiscretizationConnection +import numpy as np -# TODO -class FiredrakeConnection(DiscretizationConnection): +# {{{ 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_elt) + 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): - pass + """ + :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) -- GitLab From 947668b8fe18ba5d61a4071aaaf366f4ac6dee81 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 16:57:00 -0500 Subject: [PATCH 19/33] Fixed copyrights --- meshmode/interop/FInAT/__init__.py | 2 +- meshmode/interop/FInAT/lagrange_element.py | 22 +++++++++++++++++++ meshmode/interop/__init__.py | 2 +- meshmode/interop/fiat/__init__.py | 2 +- meshmode/interop/fiat/simplex_cell.py | 22 +++++++++++++++++++ meshmode/interop/firedrake/__init__.py | 2 +- meshmode/interop/firedrake/function.py | 22 +++++++++++++++++++ meshmode/interop/firedrake/function_space.py | 22 +++++++++++++++++++ .../firedrake/function_space_coordless.py | 22 +++++++++++++++++++ .../firedrake/function_space_shared_data.py | 22 +++++++++++++++++++ meshmode/interop/firedrake/mesh_geometry.py | 22 +++++++++++++++++++ meshmode/interop/firedrake/mesh_topology.py | 22 +++++++++++++++++++ 12 files changed, 180 insertions(+), 4 deletions(-) diff --git a/meshmode/interop/FInAT/__init__.py b/meshmode/interop/FInAT/__init__.py index 3708afa4..a79ad2ae 100644 --- a/meshmode/interop/FInAT/__init__.py +++ b/meshmode/interop/FInAT/__init__.py @@ -1,4 +1,4 @@ -__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy diff --git a/meshmode/interop/FInAT/lagrange_element.py b/meshmode/interop/FInAT/lagrange_element.py index ce71d6b7..6a80ed9f 100644 --- a/meshmode/interop/FInAT/lagrange_element.py +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -1,3 +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. +""" + import numpy as np import numpy.linalg as la import six diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py index bd6d9257..089e4471 100644 --- a/meshmode/interop/__init__.py +++ b/meshmode/interop/__init__.py @@ -1,4 +1,4 @@ -__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy diff --git a/meshmode/interop/fiat/__init__.py b/meshmode/interop/fiat/__init__.py index 3536a1ea..57a42bfe 100644 --- a/meshmode/interop/fiat/__init__.py +++ b/meshmode/interop/fiat/__init__.py @@ -1,4 +1,4 @@ -__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py index bdcabd83..cf0c1622 100644 --- a/meshmode/interop/fiat/simplex_cell.py +++ b/meshmode/interop/fiat/simplex_cell.py @@ -1,3 +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. +""" + import numpy as np import numpy.linalg as la diff --git a/meshmode/interop/firedrake/__init__.py b/meshmode/interop/firedrake/__init__.py index 11902b75..47c56e6a 100644 --- a/meshmode/interop/firedrake/__init__.py +++ b/meshmode/interop/firedrake/__init__.py @@ -1,4 +1,4 @@ -__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy diff --git a/meshmode/interop/firedrake/function.py b/meshmode/interop/firedrake/function.py index ce601ec7..14ef3a46 100644 --- a/meshmode/interop/firedrake/function.py +++ b/meshmode/interop/firedrake/function.py @@ -1,3 +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. +""" + import numpy as np from firedrake import Function diff --git a/meshmode/interop/firedrake/function_space.py b/meshmode/interop/firedrake/function_space.py index 95f4e297..998bed3b 100644 --- a/meshmode/interop/firedrake/function_space.py +++ b/meshmode/interop/firedrake/function_space.py @@ -1,3 +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 warnings import warn import numpy as np diff --git a/meshmode/interop/firedrake/function_space_coordless.py b/meshmode/interop/firedrake/function_space_coordless.py index 0a4e7bc4..3d785bc9 100644 --- a/meshmode/interop/firedrake/function_space_coordless.py +++ b/meshmode/interop/firedrake/function_space_coordless.py @@ -1,3 +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 warnings import warn # noqa from meshmode.interop import ExternalImportHandler diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index c58bbcd1..1682ffdd 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -1,3 +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. +""" + import numpy as np import numpy.linalg as la from decorator import decorator diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py index cd5c1d43..715e2336 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -1,3 +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 warnings import warn # noqa from collections import defaultdict import numpy as np diff --git a/meshmode/interop/firedrake/mesh_topology.py b/meshmode/interop/firedrake/mesh_topology.py index d4ec055e..1b16df9f 100644 --- a/meshmode/interop/firedrake/mesh_topology.py +++ b/meshmode/interop/firedrake/mesh_topology.py @@ -1,3 +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 warnings import warn # noqa import numpy as np -- GitLab From 5ed91ce0f3694f4db35c918cd17d6eeece1156e1 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 18:03:16 -0500 Subject: [PATCH 20/33] renamed to fix flake8 errors --- meshmode/interop/firedrake/mesh_geometry.py | 13 +++++++------ meshmode/interop/firedrake/mesh_topology.py | 19 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py index 715e2336..735bbc6a 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -133,13 +133,14 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): coordinates_fs_importer = \ self._coordinates_importer.function_space_importer() - V_importer = FiredrakeWithGeometryImporter(cl_ctx, - coordinates_fs, - coordinates_fs_importer, - self) - f = Function(V_importer.data, val=self._coordinates_a.data) + fspace_importer = \ + FiredrakeWithGeometryImporter(cl_ctx, + coordinates_fs, + coordinates_fs_importer, + self) + f = Function(fspace_importer.data, val=self._coordinates_a.data) self._coordinates_function_importer = \ - FiredrakeFunctionImporter(f, V_importer) + FiredrakeFunctionImporter(f, fspace_importer) del self._callback diff --git a/meshmode/interop/firedrake/mesh_topology.py b/meshmode/interop/firedrake/mesh_topology.py index 1b16df9f..3a5d7a2c 100644 --- a/meshmode/interop/firedrake/mesh_topology.py +++ b/meshmode/interop/firedrake/mesh_topology.py @@ -157,39 +157,38 @@ class FiredrakeMeshTopologyImporter(ExternalImportHandler): # TODO... not sure how to get around the private access plex = self.data._plex - # variable names follow a dmplex naming convention for - # cell Start/end and vertex Start/end. - cStart, cEnd = plex.getHeightStratum(0) - vStart, vEnd = plex.getDepthStratum(0) + # 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(cStart, cEnd, dtype=np.int32)) + 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(cStart, cEnd) + 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(cStart, cEnd, dtype=np.int32)[isin] + 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 vStart <= vert_id < vEnd and vert_id not in verts_checked: + 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 cStart <= other_cell_id < cEnd: - cells.append(to_fd_id[other_cell_id - cStart]) + 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 -- GitLab From 4e97cc81220e1ad1ed30baa80e0696261a96ee8f Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 1 Jun 2020 18:53:13 -0500 Subject: [PATCH 21/33] Added some docs about interop --- doc/images/firedrake_mesh_design.gv | 25 ++++++++ doc/images/firedrake_mesh_design.png | Bin 0 -> 43833 bytes doc/index.rst | 1 + doc/interop.rst | 91 +++++++++++++++++++++++++++ 4 files changed, 117 insertions(+) create mode 100755 doc/images/firedrake_mesh_design.gv create mode 100755 doc/images/firedrake_mesh_design.png create mode 100644 doc/interop.rst diff --git a/doc/images/firedrake_mesh_design.gv b/doc/images/firedrake_mesh_design.gv new file mode 100755 index 00000000..b5465390 --- /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 GIT binary patch literal 43833 zcmZ6z2V9Qt8$SF*k&vV!?M0djY0%z9B~7JGSq)mWho+`VMo9~mlp;x*Dj5wbQA$e^ z?Tz=i{J#JH`+nZQ?yi-j(Xkv(c@`hIR8gs!ekY75W{2L z_-sj4eRfAGSC+MrXU`U#eNpO>QCiwD*YU~rlJmTx^X#e8?hbR$)5ns&R;{}e5P7?i zfthME?WK1i?|yCaZ*mApoz!X^S@fEqBqW|cu64L~@#ee5o9CrN7*<`j*><93ppCw2 z_lDZLYlIgz)b2J9V)nrwi)QG^1>~G^a#Hg3^-X#4Y0uuhA(fSWDKA(r(?n1cKWqaQc@2eKK%BW`$lT2l8MRI^78WCQc`OI=>4CVn3#lhbac>@2RQda6OV3i z{CLdQueK$nrBsxe+1c#+FL;ax4<1ZRO$FuV?)?1u^X=T++&d=m8QEoJ+wR@Fr+Mg* zn!5VslB%uD%v!lcMMd!}wps?sOit_9um3hW8u)IcwyKIkv9`5cB{uJWusnp}sF4ws z*Nj?QR~K!(=KFq~%!&$jz0(eYv*t;zX`iGadFaUWk_8)xUjJ0=mMQA#Eg zQBhIx*RF+^lfp6)>2jVw`+A>WVu(=foAlZyZtzEkoce?J+pFh|3Ix()StaR;oa$RWGE~PVPM96*rS@E{Fw^zh= zr=_J$Oiw31zUu4a!_7eJ*f9R6yj(UbCx=%+Ky9hseXK*Yq@<*(wsy5)?qM|-m)$`@ zL8ec#>70&MZ;epJ&DPvA%yn6sapK>((^JqTai`wPx`-{MkE`}cN&UOTQ0JfZE}zR> zukAa{a^hXum0gB8iY6)Jl^KByL3t*HUkhLSQ9gQpWa-&TKQIVe3b zG5FfGZFO~ZeV?D6zx&|9<&>1IC*E8Sy?vWo_oj5byUomP1g4gjR(Nb|1MWO7Ie871 z@AWl_)q769Rn^p_MGQMOu)jAyc5E{nThxm`OY1daIZ23+Uwrs*_v^HSY*^Q}O82vE zpFc+nn>^s&v11KBqh@Zt{mhv&53Q?LiI|trw&&`ZM2U4&tBS7wGydtx_MySOa7df`_lE;64*PvwJ-mOdXKl`zkHSv5JadB}JLqkIm z%W@|17b}uDC{7J*?BUnT&X3S9u1vp+h>hJurFeNGn|j9sOC3f=#&A=S91dswQ>Rb! zbw04<`@7Jf5T2O0g-WraqM~tUx{tEto;^+t?V3hLMuTBoe%yF|Y6>e3#^U$9LQN#S zamVf1ckhOu`CJx{Ti-&dd-pE(%*WzT9Ol-ToYuQ6E7-%SvBDP#*nS2|U2kt%{_*G6 zEy^!MAeXd8=eT!wcWWZ-&7OIPvLYu^Delm{rs%gcy0fjNMNuxT`_e=omD!^c8zm$q zd6C`hKNdUiOGzacG|;&$PB~0XHH+(5eKN{3j$fYt)gF@_f0tZ@@#AFjNYXMgn{g_qYSxR#A=z{n?Tyr4?}u}>ZrwT_3$c_9T=F-4 z{w$7UVLy{LU-J_Z6&3j6GriIE{CPfh`Y=RVFv7DXMy=Ff%AH$C{XPoF(2aO}xBf#l@my0*52Gd<-=C66Ds#ni`XT!~NjTZ-uF z?q)<7GP1I2nVXwCHXPB>-nB6Q>&wpdEF#>L#vdOaY5Dk6NV3%GCttm)I{fB_1|<$> z;wDaR7_LvEHT8hEB(}kRyt_23xH!RsU;Q!?=2U+e7o&W-9`c&${@4w4q6=8J!FD|b z0U!BvFE6h+9A!bbx#;VWQBi!oy}hbdR_XaxFZ^1fzdIqsKRt2Ysd9M(k6lBI1C6SR z{mGMwNIRQx!OG>Ia#Ba`bA&RJWnkRDf4{B0y=vvscj3pYb(-;X34d2c6jZqx7*3u% z>Grj5U3dK85p#3x;kUPRF1Wilt>yBmDyQA6d-}BS;!L|4IJ86? zEhmm0JJy5)q*-PmW>(_3xv9C?P*qG_g_(sV;bgjRvWql(@rgIvuZUbtO--%)@F5Ag zKcP8AygjD2vA?_YbX`x+EooWVgd6*v+hSC1Si+Wbbxe`tdq-Pu=U+(Ym)}lK?DmH9oWkrU10FBlFUJXn0@qvVwz3A&n^+{l*Hly#dn?F0xY;k$PQA}{a7Y1TgH8jF}7W%_pHMUbx-1@2l#h$o}9bHS$e%pLvUHmPb zTlZvt;V2%KW~cvH;r-}U@<)^0!?$=3aK@5Tt!jLmM06c6_-`ZZ0jQsm|3m3^IJ9N6{D zv#@(66HtY5=gyrZ&)vO{@8MbKtM@BHUywgOKE6%tHjamfm#$~ejh=%KexJO3|9)gh zSeWNGsUM@R!gSYR903BeNhq~D^01; zz@1y>q~yyb{5+)P@4Nhy%WV4J?dFAqP$yi!E@Y{&G+j3LQO^0+%;T83A@5y@n*6)~ zG&_*5^qZ#MD;1_M6c z+_*96DUIN1t|a@sf@j`yL#_!87vyi9x8mw-KmWa1wmW4c`^@~$9=?}I%HvOO-|N&{ zY0K7$NKa4i-e5RZpx@N(EilS;_C@im`=g7ur^m(y><@Jwb4@NOUg8xO*DG(zKTKoA z&ln#NKq+?VS56XtcIQrzsi|q4=&{FleLjXX{=7VTPdL*zul(7A8$Tbfx^&UPbF4~{ zbAOzBiTFDmW{>f1HR@G?#y9UjZ!npAf6uV4s!H)ezinj4C%V>7#WJ(#2{-F0Nst{qG}T?)tXfB~Gh*Z60Wq`_boA3o_HZK}BReeCKPF&&Srdv?`3L*krY{7#mZlgm4^ zzWCw8IFx)tj$5sZr71HGo+dK_?ep&^zTczuoG(amS?}~ zetmotMi%S(_3cK(|91Y!H-5RRw2M|I+RJZFm`%r*gnZ&+mTg&_opr^NkIg1{-5Tlf zI~qS&Pqvc>@ z`|;1bg9ZlU6{~+#1+KO<|B$e#Y2J>6!_UjBJUr3RJG@ia#Q*$-KO`qBXXp=C2-Igx z&x{|_PFS9vnv%URdSSob^~!~>>m}U=HM2d7!s~5qPoC6}I@_J~erP5uxAS=0t(~X0 z73I6D7SHbRwixMf(~C`Xrpq+B?p-};Y`q#DnhotJv?Q*u^7V^2Szcb6oU`l<4mG4H;wO`;KOZZQ>&3z^%WACR=O-XaO1i2pjAB6!$XU z$#Y-Y-ZVESVh!hu#HW8W57^%d$|!6e`zfV-(MsczIsHbo>3lFN zd5hjs!fmiote%sSUr^xMS4Bm4@7pz*^E=+&l*xSW05H6cii%3xa!uyPRjXDpaB=C$ zU2eOuYvJtKvt{mMEH9&3Rrq4HuU@^{xN_{n+^;XO0Fb;23NLEMQIQwMpL}^}-^8fh zbo58(EK^~!vg=Wh2D|2#taWdkd#P3ZBrdmklj z=hcO~9F){P%VrOS(W=RQC z?d6SzuU1m4t3d=HgMkFe;ZWOkLGFvbkr54L_r85b-rI}}Qk7|XnQox8d>EgytlF)enVY+o zf`(i>)gRz0;O}_;$!%r;0hASIU)D3_Og?sYb`(W0MXzdick8qjZ>6QyPJLnjoHJi^ zzq$m!yIyhS1}RIr%7Ow>d`n+nUy(1iW7n=-)!39GbZp%QcQ2yDpa>Sk_SDr;tE;OA z=Qt-9?0;ooYD$Y>yT3el^`1^O9y;vC4YqSVsSDrS`1Y-F_sY-*@1E~J*S1v5=9HJ4^Od*% z{QNW+z?{$#F?X0G8~XKb$uff}x&Ky~hKCxdqwTNfuil>VlJQ?#=tAi}F6NEPofZ(Y zC=2N8v(6T?`&t)K<~9_0B!y>_a`=Uokioxm3fTVS-dv}#Llze7XLL=~^?qJn(fN_X-4YiC7$X?AZVhv4 z&%u-Im0a(Kc3M|g_o4F_d83hxQ^T)QDFHyUo=<;G4K=qdJCOB97h2WtLYn5Aa-I)r ze%6WC*y{55x)2q*a4$J{Qw4Lv2Hg{ci z+o*4r`|PYjL4HHt@>K4w`#6Pf^KBH`sr z^(3+Rmsi+G3a+bDO%i)jSbmJ_;*AGS+)v!5vD!RsG%{;eEDyd1F?8X1Id#sH8Dx4tFJdO z`*Y0sxf=MK_ikgQg};0~)r`_kw+fxE{9Y~(vMr%CQzJ3O$}iP}Dd}6g+(i1{l`!D& z`)8vxH_H}ACvfEXI$fC`+J1ZdnH#^!zVA+H-@vzTGc#-4RUSW*_aBnlzd!Ngqmx@e#JYjTbU1Njmzi*JNWP6j&!JE4r81ui zzxB9}l%!jEHn&{Z@7&A(`0-f$sJtqeItk+-SO^78S0@t-nN*n_4XZe1Gl6jYRMY zcMlb%fOGmoBD|?SYO9lm%~5U%3Glhs;+TbvTVn3G{XQojsmemR$oNtJi=xMmMygz6y3I38g2 z75rt%z;11qA6*jUcUeTuvywM%GK$L!UO?6aGY~g9(2zG99~J!_jP+(DUGg!4oaGLu zxw#!62TBHJKE1sqM;He0{?M(OmX;jj($VbOw{LGeUHDv`w8uDP zP7N7{c~lgyU%%cc_p3q$tiw85TJHVlKI_+=VcnP*AFq_7mwuR#5l`K+U&Y=w$^vu3 zprfPn_4{{iG0SpQIr=>mlfvTyoSdA5Kwy{ii1uJF<({5t{GAQ9b|afq-}euBsLSAx zYp?G+oz=O^Va*c|$Q6J8EQ9G*b8ryAbwUb`J*lax3h*eQx&b}}Qbga+pD|#Wd?u>b zf~3#jh#FlTuNk5zzetUS>+MgoL5}^zfVQ zJ#BX7o!D|zufv4OL>_2L7CpA6r!nmDW2vdyO$rPkb|pcyD)~#P9ZHFd`$s7v9?{nC z8MFPV>ZCDyH7O|``@aR;4}s-5hi+)9s00_YG~Crce}2!$$7k7GFI?bXT3WIndXx6IQXMvhvlsH|Q@95p?(sSU5Jl=mLEvcEBS?IsXp_qOaj%%G>IS2Hl# z2)aj0O^s$;)4U)j=Q3d|Lf(yR0VOpBcA_mNtK+E0XFpe03Gb0xp5Yl8JQ5NTJbA_i zyid;c9Qp8-RmR@IA%5TK4ulR0AXXnns$|a$uw{LSe zPU|-}HSrMA5eI4mhwRmVY$;*)M6M`}zqTL;7#z;vz1%e3ix)2vZn5Ie(ir)(Q@*Oh z);n0Dv_Q8aWCd7dF9?#;V`F18h=Lq^?qojizJ2@RQ5q9A^Qkk4rRWf9j-E?imoShT zY)tqUV2eg0w1kKU55!zuU3tJ^IxJJKRTMp3kN!?m>fC3|9fz)~g8wOc{CEes&#ht> z)}iW62tp7SJDs1;vJ{n;3Q!ssCI_1KfAPuu&?B{)hb5{B1*NINdpaDF6S*@&SAPnW zR)i2W_ckQhSzZF$%kH(wpjVLn7fm7>(>{7 zC0hacSmybQ6C8Q~Bnzkn)$Q$x9tTR3Zr!@_Ptu5CR|4C9c@tO6P;<%v__vVEOulii z?Rw*MewQFZN`oIhoGyQ5{Z_HtMMazOOV_{_LciIFl!hwV7ISY7Dop&JN9AY%RI)nW3R>84_;(KmV5Z*Ndw65iT8$jL>01lpd7Qw;r1}ii|5J7EFka) zW=7kkL0ckOP>H2PN!Y!t`tpSWX9XUBUG0;-e z($Ogs9V|Gwe{4+d?h6W}3zvbH8!gIQr~uo;qoP!h^&lB$9aG%0L$_t=-NRx>-RfWY z=$2B%ER~^O;ZBHjh*DY$Nf?Jd1d0q18AnIao^Jws+ZvPa0!Fm))hp!_Cw8FU2!+;Q z@%$n)i27}k_Ov*SYq1l$X$M$I!n|_j%Iif2<%g#`#quqmv4E2ALb}8)8ZFzG=$|_$ zhFY;x)NgTdkvzbL4I%02T#$wE0PDeIcN?q~QBdIePp)-ylN38%5&r6xHhLEua4@v= z^e;=C+D-T0iF7)I z(j_~Hd)FWWR6Kv4k(c+W&;eArZG9wjcx0rqmzQ+qpQY1iScOMZC2kE34)zg=5bKXe zC@wCBUT%Hrl#u8jR4*4-S0d2jmX?2W3J05tyhO%DBO45c9GQ}5-#%t^BX{rK4Zz`5 zk6~X2W*u}nGZ7-gH}g3*bZ7_X^ffEfDxEv`phWQmz9K5R5p0}_np)@u*+50^->VR< zVG!eAzkY3f?p*KCllb_}!1!x!-pk2ZQ(9VjH$PuhO-(K33z4z~1>;y2?%x+0wR`jW zb>Hytj-FHCw?d#+-n@B}sr%hi|G$0}zvn*nO)nzzb1^V7-g)S~$G7wS`w-+S9%!@O z+HV^h;~s2*UcXJ=N7}-|qHl1JnvIPONvsO)0K~w<;^H7^ji`)k-`~AgP@o3-89oC2 zBS(U!rrc2bGXAHnLzB^fLttuZss-9HWw(Sx;In6Phz|*{QwvK=HJ?7+aNlx2LdkOC z66Cmm&z~)bhhQo3#*L_;pjD6*bfR3Ki9MNaA!&wrbZpGi^d4@d_5FM5zSfWY0$j9a zmp8ET{#KUQzdsC_hJ?nMuH8hDFBiD;>fF7A1V#vg+t+U5tSfeCeh4<5a7bqsb>~rg zvI?A<-a$y`VWKlIHolClJ=ggJ&WnJC0VNHM+}adH)EZ=GXtvKRM`mYcvaAH*TX1o6 zOE`BHnX#t6o(5(KZSS|ywvK13f8}9q=nGu5)^m_tg$_U5%5RZI@8B#Gxk=xf;CxWD zDOD{k3@H5+3UCB1brAhA`h?~8t*9Ok-G;Pr*v@{ELS}7^shd&4?}OadH#Y9b=y#j{ zX$hdTLGQpMHBv)XI~yyjlAc~+?bbjN5r+=hI5|aRsD!;Lb!fi3)_7v*>h@^w}IytN}0n{$b0=$hccrc(O;tJba@|QesxA z!Zpux*HcsbfeKzWCGL!e_A4$f{^QGwbPtge6x@}yqhn*_^JH5ZtHY>`lr^5!DPAj+ z4X#l35B#2K_kyZX1UmSD&qADsNDk3`U`TOluokrW-#7Z)o^mYxTDu1f3p76Fk$34D z>gsFMV>z!t8q22)M4R>W1^W<1@$2ituQUZL?4qX4vn#NVQ%q;yTeU;emPe<*z?$B97`3yaqt zU8iKQ)^#uYg=0cPY39a1g@M%t@4*SCQ`%=?H?F4cjwyEs(dZm&h3+O`opi8Uj>TWo4FR1UUoekh_Mx%%U3j(U{IKrl3hyoN4$iZq8M%zbx^w$>0FW); zM`c@4pYfM;m840d94Kd?%{%~_u?x)zw6XZa#MR{3b`=F%8W(XPOH%?$Cra6c7u>8* z`X4QTW`p*F`@jKq=ZeWS#63c!-J-q91_qlchHZf?dJb`%!)uXpy=Hw4Me{CQ5-3WN>bD=v>ThIBEquTYi~$g zM#HC1Y-o}eI$!(=)l>M*1WZ6529AVVJGhw2&bd?AA2PRmTps-A$8{!iXpuxQZ zFjZ-A`G(*)z=@%Zdmd}EF|9VC?g%TC2Qw8~RwjL>>me0DTpe6gE<>-FyNmAHiTbRu zrm+%(4S>=m$l+c)nZzMt3^LUTBr+{M#S9+N zNy9J|iXels>H2gVn)Fusxw$zYlNy}#NZ<=wm_2a$3=W-LZs*T~i4TOehZ6fXVckP4 z7bWar_26n24vs4b!%RfG<*$h=e-%$Gd_2q@`X3*%-j6 zwz`@MrQ%^xQ8n7o-LSVsYOxyq3*4Z0$CcmApdNg^vcH#tjDij=Hgx?~G>r`Ki~t}d zK0WUC)Z*Ls?{~7ZVJI0@^yPh^0)pFPiZ?HKPxZzdPr5#Q$`O^ShiZ^H(^Ggo8Si6lc< zrrDuQwupT%oDOIw6tosVNS_o@exrl93)R zk5;3NjZa7jf*~d2{{8#SNu|IwH8nNw&-upy90FU)e$L6ui-K7ojQNu2OFH@58D4pL z&Xe!bl4PdT_!)_V(Q7n2MUa1t6Z#E^S=pHZBr+R-W+p?4lA@q@yAJ|EG=Q?C7JITB zs|`Jcv|K*(J)T0k9nVJFb9I)#14NNk*l?dOGWCxgO}fp+iH>8=|Ll=@H{8P*f?%xL zlp=0rY*5^L=!Vpxi%&I*OrxKC@w?&1dszxenfr7X+n%zc#>U!ze-`|@!4R}F%}e=g zOozLtC-J`pXlO7$D=X7KT-Wka>;wN=r1J^GZpv(-DlUw}r3Px+y^b$I8Go_;z{uc2E#^0c}C55G-eC*zS;r!&KB7%uEtSa8ztyl_Uxez^b99hx+q}HmPoLHxGbjSp#Vu^q$a>W5=S%Er^NH zpL}~$>C~xHDdGX7n!$OoIzVHj5L^kd`=VSRKK0GaBG4_3ei?Vp%ga+abm+?KRS$$< ztp$_P0>KEV@4ud)-*fpz0ld1WPp^j)XkyVvC^fa0R3@$J!u$7wA|lp}OK+gT6QS>` zLEUwo{Hg{krnX*F-`H3b5DpPf3tSH#e9pQ4TO5xHvCSg9a3SOcJQ(@^MW1ckw*3cP zIf^ySzZC=Zq;B=Sm!D6~Lk&UL@!UCetQOK&D;%WXzn|K>Z(jyHyq)5A=LkB*%>lZ@ zAY0qry`5N;|8I9uCog{Wp%};IQ8ha|KA>>oCkD$u0BW-F%^NkenRIk?S$5Nhbangs z`Y7n&(3{ov^jMlN{P^*M9AiASkYyl}faS}f)p6d({7XwpZUbh8hK33X2`OVe!C+*+ zc=YHIQ62Bx;T=_qgku^=J_vq9(&G#duOXeGw|6y`{HJd8nfy@kFt|}{p)TczDrmx_ z&4;uP!dewJK1oT*xYX1LG%T zsj)N<`z$QZ^b|Ov1K?*SEjj$QVC2vR73=I3B83CaQ4MQr?pMvtHs{VAmZNXqbm-6_ z)S@+nDfqqQ4kIidvS&d-!EQ;(1eWV_Y7HMhre1k{+*BZh{Ps4)Yea$pD!LYiCo{d+SLQwSJdcsy!bTH-v8s;{HvX?*=!)xzTD%(Vz` zK{cpg$k+7{@Bvdc6exJ0*vb&}Brh+Y-)G3n66HGlnw3!QVARN2gOs!T82muzk3d^VSCIiyE}mAhKYJ&?KmVBv-JXQjYrzSZDVE(X7C$jFKE&TDlL+ z{AECxJ7s0pJoY{<(o10HAKCEB=dogrKTMvR85x7%RPX2dyRa5{>7mb}oAEc+`%vaw zvDUELa$_C<_e+9g+|RZqcz>(vjluU|UP<9J3VuEUZ31{JI*XXZ+2Z(L@BY zf67jutEaFRtzcgLQ)`HZMNqP_JY@0ktA;nu>;=0Ho zlyUoZ?SnUDHQ~ZU?zBhYFUg!?yEH#I2vw(YWlE=8Vy`x&1*j64e*#Az^y&OGECm-; zDtI3%oFngY#L$edi*4WoC6mRGZDpDga}(;^$HSle%%7a&U}tB4=skU++n_w| zf?3>O{{wgO#nRqBNxFKqudh1T^r1Bsg=%?3j+$gYl^X9g{ zwY9Z!D_^|L#z%U*<>fqn#GSx82!Myl0zmtcDw|j`NEJtnmLnfi?XS-@-@Si7;T#W4 zZ_jF{sTD(%`va`lnx?N_*r!$Ai%faRNB$KrV?;P?B*Z;`SHz#2@F3vLXCJIPRBThb z=4)S{qTCvuEW7Gcf7^R59PH^xaiV9J<`%t3HX4q)607NWrt-JbO1XaXm9O*Uh}(n3 zsmaHI+0Se`(Z>eTjka@Cc=$K{zO_@`nmMBtD|1$K_aFQHzW8-~{LtCO-}Q8rH#z-R zT#0WUAua)ouML#x*uWh!QxGbCTGs5zti$!nnRoLJ^2JMC0&*G4%R{XDW+6PEzJc_g(Hw^?*8td#d4&LIeMnK(Rn3IBhFBE}X3b-8{Qb%6Lt1AHOud4-69xpGM3N$+P~)Z@iTYaS{>~q@=$i z;B;Kjw+0J@o`Upv4K3S9C#0;r3XQv?hJlQh!6$5P$?f5h5jz~N-`C&4hDk*sqg$W& zr;k;574hETm63Vcd^SOp+e~z3j-SRGWu@&CM8c2s)8F#N3*X@2QYe?^e#9TDH>B-&Lnp7mTV{In@ z8ZwEb$aF;#_NC7+{JFM@9}j_kJ}&*5;x4a$=ve--C-k2lpRGlb??+yQJYOX5BU5zi z@#;tLD`6+N1qC$=JikinUJ79#u3bV<&i(k9X;9CJhd#F|kOIN##|+(Xdt2Ke2%su@ zdYfQ3uLTL~q+z>sZ~7kl8o3PpaQGXn!Q!E%cS+X}zPAO?fuHs7YxVCl=_qWN#@c;p z?#%03a$A2@E{i>T9IytmWAfs>HAaC}TUl9se&)Fe_g)9!;iR#hfqa{RmJ;ATwjr*l zNW$;;tjUkw=WNJJRXCSsr>#^N_Q@=0rA|#xPkxOOC$MGDiP!q0?rc{B0k5=TNjZAmM*ubcd$ zKE@NwCtqFNMR-cT(!);iGx`9b+YwOz*JN4R_P-oyB;3rlC2)3lPH?8!Mft2#lW8B;V4iFDj z_1dK5WLsoml&Q7RwT`&9Nm~Vr>3XyRl7-bxf+jc;WFd%J&=AYj0Plsh zjZGkAXQXdykXO*Dayb$WGLWkIckK8IM85#u4+CavfCR~BkgTrBd5WQFb!^B2zd~bL zY`+fT3ii_!3TN*I(ym0 zSs`4be|((qU3TO*pzD45zFK&JZjx{dK!O0I+>8`G>mYC)lB}2n4#4R#7Xpf!dzyMn<9+^S`-b3%Bs;|$09+_z9 zG6hTMS*ZT}{^d`?8l~$tN9brEw7@qv_HTx95uw8kae>IiSVI^>vaKk92on(IP_5ye zf(&TBiA!D;k}NvQEH!!bd7z>YFvR}?JAj14D<*1cYR45NgZpT=R;|px_*)Ol)|w1` zO*y)DJ;;EXT3SkOZp9@Zh|349Ii_H%-@Mrf#`@o+>A`~%a2jLCtf(V12%S04{{35F z|5U?Il8ze`B!<*R^+(XP*kY=Z<8dQI{KGO&Pe9+H!emSc2@c-)a1SPRD)4lq6M6kw z*J=G8G(IpQp!r`j@hw0OlO6EGdoCbthad~Vk6b-eZ14^Ok)VhOiAJ)=fB*QB%!w0% zNmeB=%>fz>mj54(7b4;3o*on%DwxP#y|bM_UZet%7Yye^Mt**D!Oh@n|5s^kx0_r= z0fPz02NOuyhoz;VSFUU%yujJB8^Oy!@9{!d)F?A; z-I|HfPXGo6x=_No5y8&Dfam<_7J!i@8YC8Ld;5_8$%N4uAjQLq>5iBV-|r7>_3$ z<7Tk_){+6fc65Tq#M1H#NwU|ZU3tM%)~jpWqe5D|efu_H1XKVu&_Fvi*zAR{c_@s< zZN&Kevh62a7bx9n@KN97KK8 z9WdC14pUj~Wq^#F+yE|)@2aD-g+b0Qf2;s z<>i0R7tD z(-Hel&=3v}_`BCSRIerL4qrPRaq%Pw=(M!8Da*Cp22KuwQja}Vk2CSl zjdW)(@z?w}9S)(~U1pJOj|+6FIN1B3>Rc6)mA%;V_~)nfw6s+a%7`Q8^^N_TaF|N? zHqzmP3~bk%1lvKBp*H0(({NExBEq`5Own7h;o1qusiIh=u`uDxKI*A3NK?tD$6`S+ zEcdd#eG36i_^t%-U&nJkh`r2>9UaND{d=^xVBx@Y;LQU8@cz2FyVv9ANv4EX>M@@Q zSS{f=3}uhip1gVor3e&o>-C>E-@;H7h)#f%2*NEr^P8dIyePufg&cKaWa$Kx7)l-5 zo$cpk?H!0#N+?`RsZ>A~J#+I91rex-Dvf;L&?6gNqo$z&a+CvqQj)2^4?#$DP>mS& z>PTkcDzUCk1#%zQPf8znn5&7qgWuqKrH|!uA7N3yV%E0`dPETbPPgk`I#kiU7ryBe zY8#CRH0xkZKR2V(k>Lp@h|@kU$kz+wXfmY!VA&Ljv6uJ2W36 zyNqj|x>7p9{JOqA6iZFljuCn67~pM+2UBUOQ~QpE#YG!fzwk1LYz~bAETA>imq;cY zSJ*}m{8<>lO^D*iuPlzF!+Kx}a4)wzZ|#5g3$h1OD|&KtduU-iu7%1N6BC1P?1Lp6 ze*iN88&GR5hsIfQWypJ!t`=Q&iNm*K>o42Q#iNXq8u`WlF9(~{ncEy1qW^si2L-p` zhLLg3b6+aol%HQBG9M9VVRZ_F7SI-x=>hU|H=Iczv89fL#UF!Q4S&x>mt8EZ>IP&a z@V@~UB`qxw`KLmeeFWe;s)o_YvB7%pXEX*9fchx-_`ECl9(vsLfxj!$KW@s6BQ~gU zyZJg{@j77FK~YqdpI_tH7{@~nD7ctO4EW$Cq;M0fr1mDpcOqtC*RMOA%e@vfq!0Fs z;8?Ivx^gk*$Vuxktpi!P_U&N9hAu=7G4sIDh2L)N^%~enUU~3N;qQTG<}A2mB!tZj z82!+zA$GPSYZ}t7TbIG1GTc8 zh6XSwPNYL{T6R4$Wau*(!-akLq6XOHap|L2ZfcYw^4b&O1>JO1$9`DYUtZm1*!Sz# zCb(MxIMQ#isVk-VEjbaQm)J%j?7+zkc31!&DFVL6w?(gkONojyDq{^1J`O%4Anj`C z6{L*-N`3ucKqwa-DksF|1q|dkYW(l(6BE{=J~?`*ljWaxFh;a=bgWUP1#drzTKF({ zanRn`ITpW2MMN{25<_boY9>KD9X0U9t|k)=TiFjtvH&rjzQL zo7YoJh0{lqQd7;I%#VKi)(3csSC8NuWrtEFY0zwl5VuRgr@&^q$^;Bv*C3NI`{bBF zVd`^x=Eo$%PLDyox^lOxeOOMyl@nXET%&nezgY+Glc2 zHy3>V@2@FxD+0<;$jeAzMtKPabysd1J#qLT!jP>nsiTe4K^}wBrOAP;rvi5&RMmQ>}SMj9B~jC3$Sfi3NNlR=LPXIbckiP%02mozSTDmcr&6e z&gD?ci>SAS{f+`y4*p#cxdsMq+W&#vQD-x`Ux>LZnl2E*=HcE4IHlPShw(XV+^|8G z$v#jq0OPBpe*y6v&z!mX$I8-jBV=26_$C%t=&ghR|J%f4=dhM&e#ossM*4?ZA3X~G zA6W1E@gv%%_OjCa-y+WugF}|fpQ_RWxOpW7xsU2A2^j}Uzt8In|N7awKus;XJI|

>8izt!`;($TvHB=k8rvy|e?> zfPavJCKo2{z)D1Hvi-bY0|Mo4EChPZDx?>vb4Xp+7UajIPC3Dv3Y}09LvG`ObT~6X zL&L)uq}(165kdBQT)O8a(h+%a3-*GzcVIGwsTtz?+6lMTs7@cio{26298N{S+oB-Pa#wDI)P!~;@%zq` z$BEtr3n*zZ;YNf#4EGz}X$4i}#RB(kn;L33jh||8P(unQWCBRs0Z8?zvphKaL9k5Q zkxJz@h@}uJL>(C&9bG^s!VqgLuqLMSGtk6t+jC+AjFtE=*i36@=S`Sk>W=?#+`%FH zeQq|`5KMh#;GALCUS??}_Tg1c=u;7z=X}rhVW6XquqwB1iJtf4b0?>p@F-v_tI+P_ z{A%G9NM#pB86Q|Y4T4Ju2)GJefrUL@hZ5vE@UqD7qyV67E$Q{SxWwfbSWng1l-#1h-2)XAa(Md1Pd^>mW&r>4hbu@8t#Vf&8#%1 z(Jw;+f<@i+&x$wMEWEZ19ZM~^^bcPZUqA=J%V++jVy1uBFJBJ*1nL8N@&Jey@@gkV z?f^7!Xs#M^R(-1mRqqNn7@mCqROz=5d0UC+3r2y6om-Ep!4X1uKjPmb>PpNCgZeJg zzLk`gHUZr^YQ*2b8{XXJewaxDqGA3f9IVJk$a0{&b2zRTC_yy0|NG~+Jt#?PimQjm zj=ziV{Qr$zIf{a&*^H5mDezX4c~baZ@0O@|@726nvl zAdbmE1bo*7tZcwOVo8|op}ZoLp`ux$$(|xx`&7|Ef0C61Vs~jF)rv<`*tG9u_k04vPLFHW;-5$z%tlDh|K(}VwxW8 z$|Pt%J5U)V0}`iR*VL>6KcI%-#k^b~RIsX-FV_$uca}zr^h-E+{jfA0J>R$m=V%at z1U9c*LUstylREuf159px-@iv9VUzxZ0IG3E-Tk)6(4f9al_NYd0G##7lVNyW!u_Yl z4=f^}6@6dJ>JSYPTGNRoS2-m zcH%8Y)D*mr>-lpn*oR?1QZ@2(K$ld%D2>0W`}FD9gK^j}-M=4fewF;Pm`#R57uIP$>uxMBO=a59bs%GR(MH0>eC z=UOLDe}pzjKohU&3nHvpenVt75K2DQj5a#+MHG?BmV4WUEcu`2NPhmnDS0i*kIyUHIX zCnaqTR0Oj_9j~d1&`>AG6F~WD$IojLc4P!CsUz*qg?C@*n3L1f^(b+W<=bNr!{mKL z@IOtVXC$+nr)0`(dwX-EhC@;7qTvULh_Zmdie(c$-r4K%IhQt2Q4t>7aJ;PW=zR3x z7|klGrq{2p3gzQfGUVf%OyB%EfAagPzyOpkt?~U}vhFWOHVb?#_)rm!cmI4H8x!u$ z9W-?{qj^y^o~ab^MTg3q=i{ie z^=rXKN>09@$$dWJBm)=7C^FJV!YCak8e^@XXvmrz7 zKyr6y2dwaBIJ|2h0R)brlU1_Gu7%GZi})Wf=p0}F?VCq|L)`2|vrjLyt4FXu5zh{F ze+nDSU2A43|Nm5V=3zOmUH`vBBq2i-LNZjQ%tM4w#wcYjB}1YRB?_T3lrp7^A(5Fh zP$ZJVCK;khRM;U@WGH$+EBpB!$L~G%AMb;@ulu^r>pa)`PV2#(Cod~dMZ-y=_Ho<*jvGY5m0qnz34ylbI;__OY zk6d};Nc7YpMwWdOqHN6O@?(|4PU!?omT`Dr`0P_l5m z*{_D%J9F@>D>@n(&9ZmaYt{9~9!kf_#i?P4*Y}b2vB=SJ@|qmy*YOjtta}sqYTM-U zb=7ZwztYYNG`;yQe%2M^vWFi>M){x8y*e}JaI#IDp}zRdPbsfgX(%_l0cDmWyh%fs6m9xjm%ipZgfq_%g{%CAJw`hD5 zqy*q=u!kqV#3HTpSzlGy4x;eJt(>`<`#}xl5TWGt{g{u8L8xDene+WPXXlG0` zWRe=d+QSH(SfVPa8qR!))%3%)KLC9WHNr@p$}r^1-&X&cqtW$8bIj}wD7wMAtsITe z!#&E#xJgb;^Xx{+k+tg8tM~PDemnZdX=#7nQLdueBwluuPUWW&r725Kr_fD=@jU@OO1+-mR(OuPQ554B9 zS-X5D1od`&)E<7o_v0)e201`iu8cORxG;uyt^eZ8ZgMsp?fUaC>7~DpwEv&4B!{jxBS0GE*#J+N&rnN5GMQ89u3y{zo-um~Ok%7dyMGlC}=0n$HC?L|rmD;WAvbb&{lb zJlZ%6bc_$J2c4_Bc%(+D$F)Xzj=&60vp;=O(SfwP33=N*vnCIyO&j%Q&B8DCZ_akE zA<-rDE&3n3d4ZbI+?kQ$>xt&_9TH+VmxUB zG_nH+R}@~}c%hbP+>i%$V6o~{@|;sCoR*Ux#G(k+($tLpl}qcpC$wliu>=z45K(gY z)v)THjX4(UQKI?A_E}vvFE`$;0XT)Kv)N%fGbOM`Kaze-65ND3p$-^#sI+3<>{Y4* zE^k4x^)^TY`W5$0*FNKJq12`YiS(ikXqGo$oe`nc=0WKdic{53Nzqtpz_9FIiIK-RB0V zX5BE@6tMO7{_6Ml_jli&6X(Z zp;(WDw&<5!N2_jd$=J<|8m$9fwTwrP+SAThPj4+@g2g#;J?}nv!2dd};<&c%_!&7d z-=>?zzJ8m(G||t+^S*xn^=B_HE{*+e8=t_2A;x$*kl@S{;F+sGoQvC&X*=yw4{Z$~ zf^ERHg@942iyi({{S581uo0zaR-x%fMRkJXfsl}skK93(m0!KNp}!Em{!pT6vdxN} zx#_K!UekG2x7-pC7c@{>XZ#FSi%J3!_o7O+n(1j{TDc{;DxF$$lA4-Y^IywZ1t)Lk z<{tF=nqkJ`dv+wa&#OR}&E|(+-_JI+b-wO&;?qMcm};7a-uP$P-f1suc_ZP}i^-do z^6-nbq&sZF(&^i%C1rGbX{wNhRzb9>E;f7CM7c+89W*0H?w4O$QMvQlnESw}A0Hs29Zr8=dBgf`8 z^~+ropQ~(nBwt4cRpW|x2e-Ua!>JdV}liHWxQ&IbPTZU$eA@tdWm zt*GCifv(Sv{?E?^KdZZWR#=N_o1oaA$38xIe%g8QFl`(Exa$-0LPP7hd}bU>^z`B$C5yNTjmMr=8n+Bj7XhDmU}srY zlHE^?qcm>e+Vw$^LXn@iy1Lql|KI$;^md7=Gvin4n~WT3fxB=|1%>n$Rsz2Djp^^E z9CGlpp9w;{NmR~2@SFh#79Tf+j5i?ATz`pg_tkuPWzcvvT3HEY$8Fqnk9VK(cmBW8XAf2_;Dgc_qrdi>#;l4X3GOcOG3xHnGVrlb?9LcSuD;x zdy{u}1Fd9-!Gl{3`~9^UG=Zbp;@yW1Y18&v%$hacXxsrdRSO6VaSnYOH07mN`9!E= z-v4Hzp%LkZG@#_MS*_@mcc1FD^k%(wNu2N^QJ09h{KW0tu#CV=9mI#DSUr44v%kOB zo|ngrl+}MOWOm&R8bPnW9F2%bqeYu77spbb)7|(^@6%oR(rB4Bff=KsRY$xo`OnwC zzZ^-=>X`b$kOP5j#q3Rj@r;|(?6FgAB<^C}&o?czEfm+fRdto;ah%fq76=pt35cx^z7}cQaV8ZC$sl>rTLM(q{Q9bo_(2&!U3fE~I6?1V zEKAeZ|C!Rg!S#-wwc=*K_@5TQ(=)i;spP~(X9}(?Elo*HmE5!F!sz#$@hT9=bcvH^ z&sHM>|3CYf(+ZmzCp^zYMh?frW@2I@G0_B9h_`zBjdep=xh6)3XO|FC!DeQ;JMUDh zsb8LX^Q=2hI;Z%O=AUbMdRCNKIS)=f(i; zJ6^H&%FwZ_PE_m5iywCqw3F?6BjE`2g@I5gt+@hI0I@}>R-g?AkeWJTP9+(&jWACjtp7s8(R`=jEB37fI#=D=@3 z|7tjH>O0Z%h%!K~j}4|KF2ukZoJ$DLOY`ry$hMHLIgz+XJ9Zec4DcLF2D(~fBI4r; z$tc1?atA;rBW~#&F(XJu$ln|S&s;6tiIh6<*W&C!b0#iy!V3X5Bwph`Revv!FLv1$ z*7x}Q*AQk{@&Rr+Bb`XAprgc1{fgEkx=_!`+&u6~S`|qmIHAp0GrT9{*aO`R3zQIrwYL!23*2O!owYcuDmdqBrK?$_%)I0HD6ZNjevAQ=t(ZRdz;fhRDhhF zjk~`2V=Gh_t~xPnkc5ByPc(SvkrA67ot)I}-sjDOLXI|4$K2Zlms8L0OZmia z^7vDA)x!vo0|-lnaJbuiCdnUI$r=b(O3oy1KkP!@(0&91oSGc~1~XcQYTUp7H&CMY z5XbhDGN_BYq9YTi*_zdi}aqut&6OYoGDc>`=v z`)H{gI9lg8C`V zUVPyu8o#QO9C$w#I6FV+v1|&5b}gZ6fu!XvD;d}w3e~R-u^Y=?4L>N!fO}j{5fMNw zKp)bqWj@v8@?I^lJu_~E_kL7|oj=?#3l?k=XC(~>|8eW=?d%i)_==wXvCy$U+=$0O zgWvZAG~T7mX;<)jgk$UPQSUZ$c0|`Um5{C6T^*E|NF!9+{WlGF66zsqO7yIoSG^ z>be0J*H*gY47z;T%g;iIu6iTIFJEdd-0+t4G7`PH#L%`&xs@0rk7`FJ15DY3nmz?zCQNN0c#leL~RuB}@K! z_(dz_37QjV%4JX`vdLkWoe(r{$r7O!GXgVe@qv#>R}i5~)rUFNO?#e6H4G?|4jxtz zxkq&P7<=U#2L=-@j|6xVxdz!5nkQ0UXZG1UyVewqiyQMNZrAKpi>Km8NcgyQ(epw- z1)%RLmc%@w;>c~4*wn|$!^0!#?%n#;znia%ZHi@)R#_7JwSrQ2$9L-3 zF*@Fv`oq5@;XS#!K2TPbdAnpwtVvSTVb(y~c@Sfmj*{r)tv|K)p731VWCxg3VmnRE z7w=R)|Xe7e;=Ve63N)m09G9D!R+xl(OeSv_@HHJK15{t9o1mFuZt^JPn7Z*olC!SZVS z^`X&cJ7~!t?}k;Ma&0$g-}%DwpQ5HtNbvN3cYe-|M#f=1%iZ_{5sB1o*s6|+3arOA zC>BP?C@0;Y`x(2hi3*cDHuAmh+j6JLE&>5C8tF@EBCI zj1su__dd%DMk$gw%ax9#uB#4b+hCs|+(d4d6$h_;GGb0Gsxa2R>m7~7cN;|G1}jLN zFDSVpTJ!p4?3xPk57%wfW(x}iz>9MzYU$EA(UlYb6zw1{LA_a5mgz`gOibUA=;%G} zZA>4+nci->)PLJHL#L;a{o~&GSv578nX3@?+vhLeDQqnU72g~6+)fa|h-fam*KuB( zWNipEYDPA>W7@tu`@jkl1oYH&^ysah%l>RZehd>!Y&yai9w@$10^HN56N4WG$4FkjvRnyWTjPe*xk&f<6!E-FxGO3ro*b?m0O1NS8(rc(@VTN7)Ru z6iGNjFsu6$-xKAw?3wzsgGzRon+h77&sCr3qIaUm{edviqit>5LYpaFyEX>Yrpvdx zc06XGr&+czz}?%vdFB# zA$3xFvqcm+FiGj@lbr4~$bO+^tlJ-p4_s}#IcvXZU}QCw6Umi`2HiL0Dte7XzxzQK z6az}ooCR9?jEkaBfBbf1hkGZ^_H!EDVduST3s+A`X$GLY&pj(Uy8wMl3M_9xjqrRl&oe-&Iw&{WmFkp6pAS-Ox|(44Q6TMnoS8)Xn+uMwOJthD%=#xdMVf zbH$?zCNlXymz1MhYnKSUvHxM4pW3-}&wA)GjSZchfejs1*UkI)$DbZ%TlNi1wEKxM zE-z2!R~g6oouh?$T2MUOBIILYf$Ps-zqTa=bZj~Mbw$sh7`J^3RFYRpT8&n(xpkJi zx^@r|b75pOltwLYs{E#8e$?_Ny}`U>BpI=sdXSbJ)o3mVc^!jH?PP2-l{!N4Qo;+x z3q-eP&2?Bmk4hxXr4RczIXqy`tDlSD^)P8KI$ewAoUdtdztW; z3564DVq*Kp*mc+}>vGQM-RLqN2kMOewr*U&J(0NfGBj)tb&8lQVr&bq8NTO+&OW+( z$Qa7&W{QHjFZPvOQG3w11^Fv=6wsp0haphP?6AV~`~0DAOG@$I{K@W#*A1xW zf*a!#ryZRf8Ks=hP$zdaeH5(4(d0Y6B3gL_q-)JY6!SPYfNVw2${&LS2J)^xH!Npg z$v}EQ%kS%}f3H&Ph(G_x>%cT(qI)0^){#{EGY)PcF>cn&_XNJzZlM)Zirb(mM_;Mg z8q@tnTE@>~GO9ktcXxa^JTBVd+6vTkFJ3?A$ta?x`Imu|es%AQ zs>A`Y75SVyo_5(1IQe^w8-kA-6||qMxoLL_OKWWX__;xgEwNZX=Ahy~N`hH=CsxO$ zVjEBz?vyq5)5307KiZkaEp8aM*BK8=LXP`hY63-c!r}90hrCaW{lo!k<>I8A8Iim) z3*o^Q%H4eoPXii7usnZ^d$D@ijf^fV1Iw0Ut~x8%+VC_bO}Jgh)}KmQE2sRAz$PYK z)cRZfiENi*3G~c4%t_v7KEx%uNF}T~{`s_yMKXf8u`te9QdLmF0J_r9N(&SavyE{7x=S-C-rnp(rlv&~i67DO6@X*}L-q+Wc$f4~<@2 z_Z(nc`LR#-a>%;uiY5OcTQ&2pwuNqUGwk*LAJ7>*b+>6T{4b?uYA#WN6I3I+sfkxV>Ir!H6LYL3S7G0$f3g!DeR z%j7U%;ohV~C3_@UZ#K;u4UN_LSC{Sg1z1W9gHSS14;*P}iJ>Vg7gKdBl2nbKtT;#T z0M*9}GN*wyFJcGJ-l=J-b7gEk#q^rKn1CdH3J>eQ`8U>YUa(@Grd#r%*JrZzv^j4h z%1cW{*iOMx9}OIU$D_d8hOQ|ZUe>K6k|}^|Vx|C<9A4k@MJO|OilJfNIJQK|QQYC$ zbkUMMjYcYCZ~?G;ja`u2Mug^eZrjcSM2vtK-g9=gL(~&?)bGGX!GpwuhyZbtk!1)e z^ldOP5za93#=76sMtaf;WMo)WY4RH-+uFz=+zz)Gkt=_jC9Ts0i!G-~p^=aJu&u~F}&e0MfyIyrvbEb3ts1{~F(}C6sWfFax9Uqf+|9-gksVAMy50Y)`-#^Yi zeAk+~+`SXpKoDpZy^kK%k@Pksx7j-f`=%D!(;T0BxB*Ysc9PB`a}CA>R1IDKr%mm; zYDa3iKCR&1c1a4Lts{uHB)ewDS8yTH>&$g*DkUR~FS3suz-KoJ+PTXlDZ1L^l{>&( zITbPzRxOn@4qR^q?^JYzWk7csL`vQ^{rRNqB%`37XZ&3d1a7AMH1TPE#dXgA3+7=* zSAthg;f-Wc7XV*j<^Fc!jJ(0v3gs7JnA3}xyXO=+h7z0P5X~A3kCKv?Ur&8w&+U+o zNosQbr;+gCdFw6ZI2A7 zq+{0av^bB37H+HltNE(}0YTr43W!4EP+K8C^pLFIFgo&zsbM6n)@&Ru#kV_+-@Bzx zzrKAhW(*ptvg-QQL$65*J7~N%KF@8LJOYHhuoZf}W>gZJ)XO=$uR;vIYzY$wH3{#= zVgwqQM|9TUO|$w%#`NsH-0F3=8(0Dba{(GcGY7eYpju54+TVHqk4wBeOjNNNa-q?a z--1VRR}BVr@Hew6xons`hDkE|LZN{?k_@KG-RRxGC@2-oeoOAuM{EJ)rRneszYx7l z0q`19sbCiGD}Zt7$E-edjMrR^s53PU?G)Hn>}c81ch~_iFz(e4O|cA-oj>1q#HPAi zFVr_YT>?~}#DTbT|Ne1Cq2k&>Wf{z2nG#z01XCg}^Ck|HHQzrqMB+Le2_oeU*A6)KjiJ9)IyOcp0jOFG#BwiZT|eJ zTq+?>zz0rEj3!uabYFo(TB;*hEG%7@B zbPy|%ZF7L%5o0j+E_d}KjyE?dBlW(?9=OPC5#_(BbIR@8cIk(Q3ipBP&zHAZP)vL{ zte)F|-(qBx^cmvNWT=k-%1f3Qf+k}d;6K6bMd0|`eegLV!D;eo#7rq5A#n);?v9P{w#l)o}o}6W*oBSs|qg_DhWvykZ20A zb_1GQ-W_|ch?qGYE(c0W!V}pYqh_Ln7Pj-y`--`P7qYf~tgrrDZO^0r?~^h` zrVG!BAMiEIpDb#am|FBAGFKHVoB($vp-+Lvbm=u{vQ-*4mZ`Sf??xMjEBFtJ2OGnk*hz<%wBHY#y+>O^;rb4|uKmgx1q)OA3>y=^Gq zMZmt`A5E@LGF{p22a}mFMYL2o|oiH*urT#-EkKPXZBN+jJOk5O&{pJ>Dc6?jOB?qW+n?lj) ziBtw(hOp?2tFpMOKVNH+1X2z+@~&*k=wkx%*kjcm03a7c+dUtcuJE?Hs<@gD^4tw%yRk%VS!nMoLQ7* zX=B04OL0coqvf7{X|6I0Tjtx;qX>x(*SCg&Mek*`W`U(Zy|ekBFaxq*{V`V62c&khWgDqx9`$r^UM~l$oR8bqh>6Jv>ma&1zK85bFUoTvo!6*@#BC9 z^{CVMkDrH4_8!yGGQt1W`@LDbHjvoflqDDiI(>TnsXn`<^}0)I`iC?MH{?j<3~a@& z|3$cF#96PM!6c4QThQmQg+YQm0D~8BikSjV{Q?YO8!54!g0&Q|I083Aj_5Fn8Jsp+ zS3eh?j_Ke*Z2mjl_e>R(TnG0@FJJ*+S0Mz)=~p1wU?!!6f!t1M>v8qYN8WRgu~T2v z%3xXkLJGE9A&sCd)vz&0B_+`c9etT#aT5)zxsfU4_rm59D}o334((E~t9uAHl2_J& zKRoBs2lLRrpq3gM_5ihb`j&A^t*rJuE&mUJ=y41kz!4Nx@4q3o;%}9|l$|}8S_lf! z0Q^?be%z|tlQB%vCJLdM@?Jo==Q$0`=lUS%|&;C^4 z{8x;;Bioh^G6(Zq9q*?fz$WS|3W|%z%KxL|Bih?`9~>)l*LAwQEC^11HyIrKwf?s& z9T8ZAx}|(avK;ZLYce~ybdHQ24z>bo;7sw%QS|GO?HJMOdxL8)6W zU+z<)sd`CV7B<{FI=3!SmTahw|d0G+MSS1gc5BQT=-<E4l(;3ow;t>vMgZlqQuYBq~m1I zf;p*2)fiov##HDQHO!@H0;{*EfKfAW)#Am7j2;Fs?y+Xgt-QSJQ1(EnLsL&~wIx{X8QF_zWKCGqo}n z{%D|GiOnm~LkARA97R0|(={Gxz(^xiouc{6-}x@%W=QS2~U_32UI;~wrAdyAtW?MJj%XCE7-pC38@X%5(qq>LmbCAsw4`I7Sb667YA z%p8Nw?g#r5L3QFju4%Te^ld|*=nNh_`2QL4x6PI?7KD%nb3T$5qBTgK8O$9~{ZF@` zcc6wFZ{#+F!%#G$bH#`j z9llblW=LW|CBZ}IT)O@|h`js$**!24HZM1yW`>~*T%jKygAFSWfp1OvXLvFTKmH#q zXFYANNzf?|@8-k{ADM@;etmU|@l9*>+O z8Y6v%y*ue`xBtpq!Z@-bKHzM4fP) zYDSeqO3lG+wfRdPP~9``AhLhUB}x~wdTae8uBan@e|uV+=MGCx<;+dsEk=&qc5Yu2tB zIIICzLi@?n-G+F_hy(&4WAKP6(?c5Rw#309=@v*f6P|@uY4us4*pWEAb*oke>t^pI zBUWj+e})k-bKNhl5of89JAfT<+R@`1Tvj{Ufo~@&=pU#i8R;_K$UG$Py63cWzdkEv zZcb0jAQ$_{$&;dtdW@Hj1E9b^4u-HK@ZKvwt6d>sMX?nN175hzr?<5es<1O?? zw9e~U(!D?UdACGUf|hfRu04GQ2*DqWX(v+DLp)ZQt|c?Vxdl(Z57=TcgTS@+xZ96M zMs9K6GcEJl71x$iJ2mJNXH)Kc_1PtkzKFd;ifkKfjZm-uGyn8d=aMOh8mvQzDTAS8 z0GXKWY*xKK|0c4*x2lG7t!Cn-kXd7Nxt(jsgO#nmhs`LvWiV~e)t&u!<@C#sz3%38 zb(7L6nG5J;xFaSV%y(Z&3L(=_zFu8$o}xM+zbd;dnpT*>Dd9SbqDI*P2k65p=9fXcr$OnAwxS zlxrJ-*#H2DK5fKJk^BDZYw<{ts5`gpmd13@$THS$>^9S47TYOn4}xeg66d2%k=@#0IK2@u8?{{1nM z4&D3K9nRW^CfY{dREUZ;Rpb}cZKy(h-1|l(tgY?5J@L`MkY(OuV%?H|b z(dn6i$mYdaOV)pKkB^p?78V;A7UTd0N-cbg=?KdRGGnqvxY^fd_iOHOR+iY&q2qljedA86p#lNTZtZG56K98?o5KT%%Txof@l|osEQcu!wT38AWZYTp+>9c z03;*mJPJ)8w_-zsb8~L2lP@X_7~YLz0HzcWeh@+;qdFcDkr5b~NeIgqDGPztphH0x zQ&_t)+KP~(JsO_%aV8<^nv&CGuwb1}Cj%4Wv)ce%0k&G3L@a)4I!Xj-wEK5m*c!LX zxN#;dHx#_fS&s~9Qe&{-E$~wT+p?QjXX00od`QZL7@OT8>@nw%xx8X6OB|pmcmb4O zV^@fdr^X3=# zXDNwU3IYgS3L>7`TzV?~orLX)LK8(xw1dHH5{(JdR43c%-+vt?n-JP3RXV!&?^OxD zV2faAo%;et*b#GxGW~I2**&QOBEVgB2M=zAJIahmACz+-mA6X-aFHQLlH*ODB(lcv zgB8!L!KHu&z{W)H0Q4ny2z2`<{VjEQvRoD~8*IR%z%#MPF&^DmzuxQHvuRGAH(P*C z735I80N7dP@z#9ZJ{tsEj zG0mE~1*7g|cmRA)RLeZi>`MvTdCw#)5#(TAOJ^-pZBmc@8!V*Jm04}1Vf^u18fLr} ztXCNh6MEzpp%IiF>eYiF8W5_W;%mpY_UWnbHIyNOy~xBeef=i+5>JuHuhZ?)r5;y= z3BT@_vIb5g{nUUpHe2z>>Sg5;oTy3MZ=~NtTKJN zo?HJz1ow#RZ|>}A)1JRK#tn0Xy&IMdTu#7d3Ovqqrxl0W^Jk?SIDs_L0r^2_J$-MVN{Qra65}EffFKHdn}R;>3SNrZ zI(0Z}B})J(;m4UDvj;TzzV!$%oj~GSw~h_zq4$}9PH1j%eghIG{Hd#@59B3BP~y+G zt;xUk0r@GetD-|t)y^NsoG0zSX_xX#WiOB4j%{9G?izw*r53=^lkEZ8 zMlGmGl9(-(vNv$gZp4OmB8qoNLmsL2hISkrh_pvyJ%`d-#=dTq4W3!r?QkV}u56f;V$KBq;A1qfcm}8wzpMR!jkuid{*A zlpIX=#4fX_1`s-Dlsx=GuN2z$4paxoj>Mmm9!#FvI5O1!<6EdyyK*|7sK-er$5-uk z?b=E18CrycY#}bI+#AaaGc=RKM?AP9T_dOFC^7(`GMQDx{k4w|0mEe{P|$47_diiJ zjaiql&9V>i2;`p>-$c9Hq*beJfIs#k=q!x%VFQv>5Jm(Y(H{Lto@U#JY+@4cC;Se~ z)}z3Y+RxXg6S!-((}>ITyU7u7alwODGO3SKCAx4jVo>IPhPq@Os(rDiZ7;R*IdY>I z)jP_2@56M0pqjSM~GZ& zPV%T`U&qOD-75sKO^&|B&`uf?fun0Q#(}p+o0{HY$cMX{FMxx6=u~nd9i+#>Urz)M zvL&WCJ>_{dLR1+YuGgGsLCq5reK>h*p$;nG@k#)Y_SjF9&f6%A3k6G;Fc&j>3#%b>e6p%vI|UlOnh5}RzO{iwP+Zq>UH6T^PCOUi$3+{yLSnIpP{QS);E#UL#0h8JuPR(s|()^GK zB6x|YJbz!yjn}2s&^A;2YK5uy){}?Qa*o#qByz8abeSE}X9>6iiI#~~@1WpD=ruRr zl62D8%Oq_FGUFloTjtpIvD(_UYzHDq2yvaXR$s)kI(l zaeyj7Rg*FO)Pxv67?6M4W9{BOZpB|a29M)sb*k$C_#D3%Wo5GYLoF4Z_lqF>?Sc1Cm2CKOZXrM5jyTB${e(eR{;CA3i(mD+%I5J3%;vH*fOM z%a@xa`6rfQYgsY-kSpEj+D}h1Gf^!Xi6Ka1>hZzN<)@0@zrQ76Jye^8FrUHS*VZ2H z>1jlvGw|&BrE8TF0cORtMB@CI6pd4jKaziuDWESt~yvXpQM`|AAj9If9djya}%MHEAcBNs}*8_1V)q38EEnwebf4|cx zb_`lF?xNkN&t(^fQqFrR9d$sVlTTZ1s(Dy^Q5juFnfRjn5YN+R z&-%bA*w3$xl(wT%Olg#`CZ~3mGwgwb>d=kTH-tqx6Q!z>W zIkd^jpoae?L(_G3u+!vM%}o3zLT$D4`l1ArM_OOr2mPy9U~=6)uota#cQ-LJtHytQMVEl!hWIX-K8-(75;U1H4;|$-z272M<=-)|EFrk30dQb- zpmrt6jOtK4L6G-@=S`U+x8TlgnueOxd?@lH|49l+ceEdHpOBDShZJ>cTHm|n+CP7% z>1`Wwgs%u4-jK6)!v4Nv|3temiHd?jR5>7kN&Yd5FhwaT*mUgzj&e(yZXLC*b`$>; zhH;LBM(2`cL^Ow3T22_+Ho|G8_=N^Ww0Sct>uP?JvF|s^^Y2XNd1B=!BN^9#o;bOoYCUL5CSAwpCJ;=a6^TpvwU3~O91&QP!>qEljugSmEG$qjW3tdFKU>r z^`G*vYRTE9&$0`!hxhz}m%SG?V(_!!d^yiouZ@JBAHjW5O+B0eIcoV^N zq`-7|yw>#c{HDE?NEy*3bPc)$k5ekB31d%&+ovl*UrRwHsq?g-k`jq@pKz)LPjikX zb(yzuEF-n(qXCYV(K$o$+fhf(4(az=3pXo8h75WKZ#s9{g=&%m|yoJ`>d0 zPelm+@bFO4v5U1Ob4&f3eV;M&z4TMU=Z*9uk_rdIZ{{^|&8;r;YF$ABAiBql<(?SLCv zW&QrGAEw?A&A%^|*j!XNEUJBm_U(|wu#Fjie@?6#9lq1mz;jUGP&SXUvg#valoU|y zvVN?z48wq4*Q)i z9EgJ7%A97^sdMLw!Dbe8G8a}TeO*1&jcE8=q-yPhJ#6NmvdK8;7gyL`w5s*md3jU3 z#(apCs05g5EOi;qp2$38kS<6sf|Rtl@g19TBr#Oc;80zY4m%}_7If8T?mu@A|9Cd^ z7rQDcV=_w-#Ufq&`lJ<+<2S1}lq40v!WtN&WrgAkwj2gxU7oh;K1hm^!sO`1oC@@h zj19V%oVc)>aGGmXGuOIRLeQ{8gm0w&(f*ddsN0O`(_aAu(lOt`;81tfW6iY%?R_F2 z!DIkxmrv)&dL9v=qUfxy67TH4b*sJPjm1F@-d|5N)%EzrYyOwxM{h2+Ka74R3Me;y-Cd!fX$k)w}pT7O3Pqx3fH@S)q{DNKxj}cJ#JngZb5a^O{2*y;+=A zF@#q`@4DY$RTUMvRAS!GuN=k(;ADW4n-Y&rZFWfS_mVun`byU_it7mz0s#6OP{u;p&RpC1l9vv*3X~Bb`^P&zO(hALxdDpqw6S^Rh z4#ie@Q<9X8*`7vC1bbS1JFQL7m%x05LN#jas1a7ZvF4`OWj6@jd2Cy?rm}27-Hdm~WW0rPZzn`}*usKDYHP@YMT$CiRXabs2@j zz*t-Hg|c(m?D|?kn1};VT7Ao|`GMXHz6ZA@?IG1oU!mAJ#)s&RpC-CsJeZFEjX804 z@Z{9}8Y2BAR@7?nMt`cGes5mC)(PtoI-Y-p_duQYC&~8gr46J+)=~ha4m$7hFVN1F z9`>UInHUpc8F6{*`(=|G>&$D82T&p>4`y=A@NrKUYtO{dE`u|XW0$wu%>L0{Uw`X| z-lk-gwPq1V+Q%=)n+>*uzdx=wSl244YPV}p_zEFGKHCvLD~eiu(* zq_hj&6M9U`Gohq0C0nNKkVE{Z7CQ4LbO-Tw&WiHg9ZMu&aF0F%38QVrUh)DZi(j&f_v&>ir<6wM<=%o{F;V3(-ES5@xtDn+xcm0VE4C%WOTtrKH*a%h z!CVF|THSi+g>xR}6DK~{@iHVsDkTgX)2VD9|CT3z37zVlcjbd-+>c#oXpr-~JVnz% zc&+#oOK?>=d|a<|Rx>neS+MWQIVFu3_4VV99~-Yhyi;X+w|mmZtOMm$yuI>|?E zlKEc$+`G38nbM?)k`^2+W@poEA$_|Mv{roa0q4!n?Hw4U_WZjjAIiEcud4pHP0EL_ zY9Ed+X&t#XYr8?`&L^HfIC`Ut2LYDUYJT7AHX6BbVV7B;uR<7kvI`w~c(k802R<#M zO8xpH2wo6@27;UO9Y8m^_yKsah@g1P2uA;kWS}w;vy9 z9({+_Ys|-vmM-#iBSx1qB4&;Ax+~-(k~j&17(Vn(qy?>!gyce^K+*dh2k| zOTIB|moM>VUaEGz*Aa3nGb~;tt~EkIBiiS#UC%VJ*PzIEN zqrQ&s_E?_LpVFZ9NtlLYv=c)!i0!P;$et#14#^SQWR!rj->wRJuMB96wbrRb#n{gbJW2+BC|`0n$xYSj4c zxf6&~M4OFL-52%s&9pQ}tDK`!^Y`xC_cjHV#Lx5lb@``wg@oc|3elAwOHHni+2w~u)li<)yzosl%&nB_&m1Ohm+Et5F2AhpNkvHO{!Sp`CW~f6EU;xW+14Sq` zyRaBk0`&mMe}%@5uW7DWafNfLf3+0R4v4e>-k~(<+?TCTe&_E`@f65A@z1ZtfatQc zxR-K&==Cm91cAH-tZw*q^E4mz2&rTVBthf&3U-$Ynt-D^psV0*=`@NX^_oi~N}I<} zl+kdGNN*K;R*)l`=5G<@oZd(h@d5eIXiY9%b z29;k9mRb*j@8NrzCm*$$>_2$Jw$g?pOylTNr(I!BOhLVg`#jxk>%o8Dzc)bvDs)tBP>TPs60?Z>kogFf$^Knak zZYzMkOvi@Y5|$2{9R!?@^`((}KEI--6Tbw=DV+s%2@t@B$$E9LGLh*cY5|$ni`H7Q zU;yWyziA{!K%(iX*`f4+#lHcH_$49f1o6&hkXXM!HPdN-rQ(D0XbM8_(?7!2N7`*JYm8>2zE?a_1FOI z7#Ac3>C*L;v&pIdoqJTH8V4l zDr_UzYt?*vdnHgUVR{Juogh`s-@R z0K+}6qyW-jPXq|`9o3<@ybb%&&xW_TM2f<)Nm@~p+UOn%cD*L+CtxvSyw&|v%33<` zlMNNNRlRn1E6T~qAqb^?bsfdE>7N)!(+&VcCWgnveDl#BTSu-V&1b}vuIeUfbw=`_. + ``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 -- GitLab From 2167d8d98c1b6200a9a67cefe3ba8d7bee06393f Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 2 Jun 2020 11:37:11 -0500 Subject: [PATCH 22/33] Finished adding firedrake mesh design to the docs --- doc/interop.rst | 133 +++++++++++++++++++++++++----------------------- 1 file changed, 68 insertions(+), 65 deletions(-) diff --git a/doc/interop.rst b/doc/interop.rst index f621b49b..41ad3b63 100644 --- a/doc/interop.rst +++ b/doc/interop.rst @@ -19,73 +19,76 @@ FInAT Firedrake --------- +TODO include this + + +Implementation Details +---------------------- + Firedrake Function Space Design ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -.. note:: - *(If you are just an end-user, you probably don't need to know this)* - - 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 +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 -- GitLab From 5d9f5d8eeb6652a15ac3722b9e48e5d7c8e1a776 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 2 Jun 2020 13:09:44 -0500 Subject: [PATCH 23/33] Some bug fixes (mostly from fd2mm conversion) revealed by testing --- meshmode/interop/firedrake/__init__.py | 3 ++- meshmode/interop/firedrake/function_space.py | 6 +++--- .../firedrake/function_space_coordless.py | 16 ++++++++++------ .../firedrake/function_space_shared_data.py | 8 ++++---- meshmode/interop/firedrake/mesh_geometry.py | 4 ++-- 5 files changed, 21 insertions(+), 16 deletions(-) diff --git a/meshmode/interop/firedrake/__init__.py b/meshmode/interop/firedrake/__init__.py index 47c56e6a..33daeb39 100644 --- a/meshmode/interop/firedrake/__init__.py +++ b/meshmode/interop/firedrake/__init__.py @@ -163,7 +163,8 @@ def import_firedrake_function_space(cl_ctx, fdrake_fspace, mesh_importer): from meshmode.interop.firedrake.function_space_coordless import \ FiredrakeFunctionSpaceImporter mesh_importer.init(cl_ctx) - finat_elt_importer = FinatLagrangeElementImporter(fdrake_fspace.finat_elt) + finat_elt_importer = \ + FinatLagrangeElementImporter(fdrake_fspace.finat_element) topological_importer = FiredrakeFunctionSpaceImporter(fdrake_fspace, mesh_importer, finat_elt_importer) diff --git a/meshmode/interop/firedrake/function_space.py b/meshmode/interop/firedrake/function_space.py index 998bed3b..85b3d6c2 100644 --- a/meshmode/interop/firedrake/function_space.py +++ b/meshmode/interop/firedrake/function_space.py @@ -88,7 +88,7 @@ class FiredrakeWithGeometryImporter(ExternalImportHandler): mesh_order = mesh_importer.data.coordinates.\ function_space().finat_element.degree - if mesh_order > self.data.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.") @@ -98,7 +98,7 @@ class FiredrakeWithGeometryImporter(ExternalImportHandler): self._resampling_mat_mm2fd = None def __getattr__(self, attr): - return getattr(self._topology_a, attr) + return getattr(self._topology_importer, attr) def mesh_importer(self): return self._mesh_importer @@ -118,7 +118,7 @@ class FiredrakeWithGeometryImporter(ExternalImportHandler): if self._resampling_mat_fd2mm is None: element_grp = self.discretization().groups[0] self._resampling_mat_fd2mm = \ - self.finat_element_a.make_resampling_matrix(element_grp) + self.finat_element_importer.make_resampling_matrix(element_grp) self._resampling_mat_mm2fd = np.linalg.inv(self._resampling_mat_fd2mm) diff --git a/meshmode/interop/firedrake/function_space_coordless.py b/meshmode/interop/firedrake/function_space_coordless.py index 3d785bc9..e6fdec6c 100644 --- a/meshmode/interop/firedrake/function_space_coordless.py +++ b/meshmode/interop/firedrake/function_space_coordless.py @@ -82,6 +82,10 @@ class FiredrakeFunctionSpaceImporter(ExternalImportHandler): 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 @@ -95,7 +99,7 @@ class FiredrakeFunctionSpaceImporter(ExternalImportHandler): raise TypeError(":param:`mesh_importer` must be either *None* " "or of type :class:`meshmode.interop.firedrake." "FiredrakeMeshTopologyImporter`") - if not function_space.mesh() == mesh_importer.data: + 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.") @@ -114,8 +118,7 @@ class FiredrakeFunctionSpaceImporter(ExternalImportHandler): # }}} - # We want to ignore any geometry and then finish initialization - function_space = function_space.topological + # finish initialization super(FiredrakeFunctionSpaceImporter, self).__init__(function_space) self._mesh_importer = mesh_importer @@ -170,6 +173,10 @@ class FiredrakeCoordinatelessFunctionImporter(ExternalImportHandler): 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 @@ -189,9 +196,6 @@ class FiredrakeCoordinatelessFunctionImporter(ExternalImportHandler): "attribute and ``function.function_space()`` " "must be identical.") - function = function.topological - function_space_importer = function_space_importer.topological_importer - super(FiredrakeCoordinatelessFunctionImporter, self).__init__(function) self._function_space_importer = function_space_importer diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index 1682ffdd..9bd03f27 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -216,15 +216,15 @@ class FiredrakeFunctionSpaceDataImporter(ExternalImportHandler): Note that :param:`mesh_importer` and :param:`finat_element_importer` are used for checking """ - if mesh_importer.topological_a == mesh_importer: + if mesh_importer.topological_importer == mesh_importer: raise TypeError(":param:`mesh_importer` is a " "FiredrakeMeshTopologyImporter, but " " must be a FiredrakeMeshGeometryImporter") - importer = (mesh_importer.importer(), finat_element_importer.importer()) + importer = (mesh_importer.data, finat_element_importer.data) super(FiredrakeFunctionSpaceDataImporter, self).__init__(importer) - self._fspace_data = FunctionSpaceData(mesh_importer.importer(), - finat_element_importer.importer()) + self._fspace_data = FunctionSpaceData(mesh_importer.data, + finat_element_importer.data) self._cl_ctx = cl_ctx self._mesh_importer = mesh_importer diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py index 735bbc6a..216cb7b1 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -91,7 +91,7 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): fspace_importer = coordinates_importer.function_space_importer() topology_importer = fspace_importer.mesh_importer() - if topology_importer != mesh.topology: + if topology_importer.data != mesh.topology: raise ValueError("Topology :param:`coordinates` lives on must be " "the same " "topology that :param:`mesh` lives on") @@ -138,7 +138,7 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): coordinates_fs, coordinates_fs_importer, self) - f = Function(fspace_importer.data, val=self._coordinates_a.data) + f = Function(fspace_importer.data, val=self._coordinates_importer.data) self._coordinates_function_importer = \ FiredrakeFunctionImporter(f, fspace_importer) -- GitLab From c6099406caa241b8a2b39c16701dc7aa3e1d0322 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 2 Jun 2020 13:10:07 -0500 Subject: [PATCH 24/33] Added an initial test which makes sure conversion is idempotent --- test/test_firedrake_interop.py | 121 +++++++++++++++++++++++++++++++++ 1 file changed, 121 insertions(+) create mode 100644 test/test_firedrake_interop.py diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py new file mode 100644 index 00000000..9661a548 --- /dev/null +++ b/test/test_firedrake_interop.py @@ -0,0 +1,121 @@ +__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 pyopencl as cl +import numpy as np + +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, sin, exp, pi, as_vector) + + +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=[1, 2, 3], ids=["P^1", "P^2", "P^3"]) +def fdrake_degree(request): + return request.param + + +# {{{ 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) + + print(np.sum(np.abs(fdrake_function.dat.data - fdrake_function_copy.dat.data))) + print(type(fdrake_function.dat.data)) + 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_getter, 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_getter() + fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_fspace) + check_idempotency(fdrake_connection, fdrake_unique) + + +def test_vector_idempotency(ctx_getter, 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_getter() + fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_vfspace) + check_idempotency(fdrake_connection, fdrake_unique) + +# }}} -- GitLab From 30a19e717711f695f90bf86207d4a19983b32cf5 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 2 Jun 2020 13:14:05 -0500 Subject: [PATCH 25/33] Some flake8 fixes, ctx_getter -> ctx_factory --- test/test_firedrake_interop.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py index 9661a548..c1f94b39 100644 --- a/test/test_firedrake_interop.py +++ b/test/test_firedrake_interop.py @@ -20,7 +20,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import pyopencl as cl import numpy as np from pyopencl.tools import ( # noqa @@ -40,7 +39,7 @@ firedrake = pytest.importorskip("firedrake") from firedrake import ( UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh, FunctionSpace, VectorFunctionSpace, Function, - SpatialCoordinate, sin, exp, pi, as_vector) + SpatialCoordinate) CLOSE_ATOL = 10**-12 @@ -87,7 +86,7 @@ def check_idempotency(fdrake_connection, fdrake_function): np.testing.assert_allclose(mm_field_copy, mm_field, atol=CLOSE_ATOL) -def test_scalar_idempotency(ctx_getter, fdrake_mesh, fdrake_degree): +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 """ @@ -98,12 +97,12 @@ def test_scalar_idempotency(ctx_getter, fdrake_mesh, fdrake_degree): fdrake_unique.dat.data[:] = np.arange(fdrake_unique.dat.data.shape[0]) # test idempotency - cl_ctx = ctx_getter() + cl_ctx = ctx_factory() fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_fspace) check_idempotency(fdrake_connection, fdrake_unique) -def test_vector_idempotency(ctx_getter, fdrake_mesh, fdrake_degree): +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 """ @@ -114,7 +113,7 @@ def test_vector_idempotency(ctx_getter, fdrake_mesh, fdrake_degree): fdrake_unique = Function(fdrake_vfspace).interpolate(xx) # test idempotency - cl_ctx = ctx_getter() + cl_ctx = ctx_factory() fdrake_connection = FromFiredrakeConnection(cl_ctx, fdrake_vfspace) check_idempotency(fdrake_connection, fdrake_unique) -- GitLab From d2871c40a9f0fe36cb6a1e34a049d52e6a282184 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 2 Jun 2020 13:47:46 -0500 Subject: [PATCH 26/33] Removed dependency on abc --- meshmode/interop/__init__.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py index 089e4471..1421458c 100644 --- a/meshmode/interop/__init__.py +++ b/meshmode/interop/__init__.py @@ -20,8 +20,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from abc import ABC - __doc__ = """ Development Interface --------------------- @@ -33,7 +31,7 @@ Development Interface # {{{ Generic, most abstract class for transporting meshmode <-> external -class ExternalDataHandler(ABC): +class ExternalDataHandler: """ A data handler takes data from meshmode and facilitates its use in another package or the reverse: takes data from another package -- GitLab From 7af97f693866c5838037230daf8dd6763ff1e608 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Mon, 22 Jun 2020 12:56:16 -0500 Subject: [PATCH 27/33] Added a test to check some basic consistency --- test/test_firedrake_interop.py | 48 ++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py index c1f94b39..eec00304 100644 --- a/test/test_firedrake_interop.py +++ b/test/test_firedrake_interop.py @@ -62,6 +62,54 @@ def fdrake_degree(request): return request.param +# {{{ Basic conversion tests + +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. 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 + +# }}} + + # {{{ Idempotency tests fd->mm->fd and (fd->)mm->fd->mm for connection def check_idempotency(fdrake_connection, fdrake_function): -- GitLab From 6a7a21e7897977b043331bbd54b654d1e43e5395 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Mon, 22 Jun 2020 12:57:07 -0500 Subject: [PATCH 28/33] Fixed the unit node coordinates routine --- meshmode/interop/FInAT/lagrange_element.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/meshmode/interop/FInAT/lagrange_element.py b/meshmode/interop/FInAT/lagrange_element.py index 6a80ed9f..858bbf01 100644 --- a/meshmode/interop/FInAT/lagrange_element.py +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -59,11 +59,16 @@ class FinatLagrangeElementImporter(ExternalImportHandler): # Check types from finat.fiat_elements import DiscontinuousLagrange, Lagrange - if not isinstance(finat_element, (Lagrange, DiscontinuousLagrange)): + 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` or" - " `finat.fiat_elements.DiscontinuousLagrange`", - " not type `%s`" % type(finat_element)) + " `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" @@ -90,6 +95,9 @@ class FinatLagrangeElementImporter(ExternalImportHandler): node_nr_to_coords = {} unit_vertex_indices = [] + # FIXME : This works, but is very ad-hoc. It is probably better + # to get permission from the FInAT people to reach into + # the fiat element and get the nodes explicitly # Get unit nodes for dim, element_nrs in six.iteritems( self.data.entity_support_dofs()): @@ -100,7 +108,7 @@ class FinatLagrangeElementImporter(ExternalImportHandler): # Record any new nodes i = 0 for node_nr in node_list: - if node_nr not in node_nr_to_coords: + if node_nr not in node_nr_to_coords and i < len(pts_on_element): node_nr_to_coords[node_nr] = pts_on_element[i] i += 1 # If is a vertex, store the index -- GitLab From 46d9fca3b43a7fcbe210ae3f6147765d54e93b26 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Mon, 22 Jun 2020 14:56:34 -0500 Subject: [PATCH 29/33] flake8 + fixed careless bug by preventing memoized functions from both calling each other --- meshmode/interop/FInAT/lagrange_element.py | 3 +- .../firedrake/function_space_shared_data.py | 1 + meshmode/interop/firedrake/mesh_geometry.py | 80 ++++++++++++------- 3 files changed, 54 insertions(+), 30 deletions(-) diff --git a/meshmode/interop/FInAT/lagrange_element.py b/meshmode/interop/FInAT/lagrange_element.py index 858bbf01..b3d4c50f 100644 --- a/meshmode/interop/FInAT/lagrange_element.py +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -108,7 +108,8 @@ class FinatLagrangeElementImporter(ExternalImportHandler): # Record any new nodes i = 0 for node_nr in node_list: - if node_nr not in node_nr_to_coords and i < len(pts_on_element): + if node_nr not in node_nr_to_coords and \ + i < len(pts_on_element): node_nr_to_coords[node_nr] = pts_on_element[i] i += 1 # If is a vertex, store the index diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index 9bd03f27..52c0e6d2 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -146,6 +146,7 @@ def reordering_array(mesh_importer, key, fspace_data): (order.shape[0]//nunit_nodes, nunit_nodes) + order.shape[1:]) flip_mat = finat_element_importer.flip_matrix() + print("ORIENT=", mesh_importer._orient) reorder_nodes(mesh_importer.orientations(), new_order, flip_mat, unflip=firedrake_to_meshmode) new_order = new_order.flatten() diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py index 216cb7b1..8ce88f26 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -255,50 +255,72 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): return self._nodes - def group(self): + def _get_unflipped_group(self): """ Return an instance of :class:meshmode.mesh.SimplexElementGroup` - corresponding to the mesh :attr:`data` + 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.) """ - if self._group is None: + # 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 - from meshmode.mesh.processing import flip_simplex_element_group fspace_importer = self.coordinates_importer.function_space_importer() finat_element_importer = fspace_importer.finat_element_importer - # IMPORTANT that set :attr:`_group` because - # :meth:`orientations` may call :meth:`group` - self._group = SimplexElementGroup( + 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()) - self._group = flip_simplex_element_group(self.vertices(), self._group, + 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*. + 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 @@ -311,9 +333,9 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): from meshmode.mesh.processing import \ find_volume_mesh_element_group_orientation - orient = \ - find_volume_mesh_element_group_orientation(self.vertices(), - self.group()) + 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 @@ -339,8 +361,8 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): 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. + NOTE : This should be guaranteed by previous checks, + but is here anyway in case of future development. """ assert self._orient is not None -- GitLab From 910303e3e89c29ddeae95b9113d5626060e9e62d Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 23 Jun 2020 10:44:14 -0500 Subject: [PATCH 30/33] Added some function transfer tests --- test/test_firedrake_interop.py | 75 +++++++++++++++++++++++++++++++--- 1 file changed, 70 insertions(+), 5 deletions(-) diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py index eec00304..f3337c4a 100644 --- a/test/test_firedrake_interop.py +++ b/test/test_firedrake_interop.py @@ -21,6 +21,7 @@ THE SOFTWARE. """ import numpy as np +import pyopencl as cl from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl @@ -39,7 +40,7 @@ firedrake = pytest.importorskip("firedrake") from firedrake import ( UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh, FunctionSpace, VectorFunctionSpace, Function, - SpatialCoordinate) + SpatialCoordinate, Constant) CLOSE_ATOL = 10**-12 @@ -57,17 +58,22 @@ def fdrake_mesh(request): 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 tests +# {{{ 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. Also ensure the + 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 """ @@ -110,6 +116,67 @@ def test_discretization_consistency(ctx_factory, fdrake_mesh, fdrake_degree): # }}} +# {{{ 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): @@ -123,8 +190,6 @@ def check_idempotency(fdrake_connection, fdrake_function): fdrake_function_copy = Function(fdrake_fspace) fdrake_connection.from_meshmode(mm_field, fdrake_function_copy) - print(np.sum(np.abs(fdrake_function.dat.data - fdrake_function_copy.dat.data))) - print(type(fdrake_function.dat.data)) np.testing.assert_allclose(fdrake_function_copy.dat.data, fdrake_function.dat.data, atol=CLOSE_ATOL) -- GitLab From 41d1ea91785d6da61a7a9c8b817d466caed22f1a Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 23 Jun 2020 10:46:06 -0500 Subject: [PATCH 31/33] Redid cell and coordinate handling a bit to account for 1D bug/get correct unit nodes --- meshmode/interop/FInAT/lagrange_element.py | 79 +++++++++---------- meshmode/interop/fiat/simplex_cell.py | 30 +++---- .../firedrake/function_space_shared_data.py | 1 - meshmode/interop/firedrake/mesh_geometry.py | 8 ++ 4 files changed, 55 insertions(+), 63 deletions(-) diff --git a/meshmode/interop/FInAT/lagrange_element.py b/meshmode/interop/FInAT/lagrange_element.py index b3d4c50f..4b3d977c 100644 --- a/meshmode/interop/FInAT/lagrange_element.py +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -22,7 +22,6 @@ THE SOFTWARE. import numpy as np import numpy.linalg as la -import six from meshmode.interop import ExternalImportHandler from meshmode.interop.fiat import FIATSimplexCellImporter @@ -44,6 +43,8 @@ class FinatLagrangeElementImporter(ExternalImportHandler): :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"``) @@ -51,6 +52,7 @@ class FinatLagrangeElementImporter(ExternalImportHandler): :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 @@ -91,39 +93,33 @@ class FinatLagrangeElementImporter(ExternalImportHandler): 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 - node_nr_to_coords = {} - unit_vertex_indices = [] - - # FIXME : This works, but is very ad-hoc. It is probably better - # to get permission from the FInAT people to reach into - # the fiat element and get the nodes explicitly - # Get unit nodes - for dim, element_nrs in six.iteritems( - self.data.entity_support_dofs()): - for element_nr, node_list in six.iteritems(element_nrs): - # Get the nodes on the element (in meshmode reference coords) - pts_on_element = self.cell_importer.make_points( - dim, element_nr, self.data.degree) - # Record any new nodes - i = 0 - for node_nr in node_list: - if node_nr not in node_nr_to_coords and \ - i < len(pts_on_element): - node_nr_to_coords[node_nr] = pts_on_element[i] - i += 1 - # If is a vertex, store the index - if dim == 0: - unit_vertex_indices.append(node_nr) - - # store vertex indices - self._unit_vertex_indices = np.array(sorted(unit_vertex_indices)) - - # Convert unit_nodes to array, then change to (dim, nunit_nodes) - # from (nunit_nodes, dim) - unit_nodes = np.array([node_nr_to_coords[i] for i in - range(len(node_nr_to_coords))]) - self._unit_nodes = unit_nodes.T.copy() + + # 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) # }}} @@ -135,9 +131,10 @@ class FinatLagrangeElementImporter(ExternalImportHandler): def unit_vertex_indices(self): """ - :return: An array of shape *(dim+1,)* of indices + :return: A numpy integer array of indices so that *self.unit_nodes()[self.unit_vertex_indices()]* - are the vertices of the reference element. + 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 @@ -203,12 +200,12 @@ class FinatLagrangeElementImporter(ExternalImportHandler): 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) + :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), \ diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py index cf0c1622..c6c0ca2c 100644 --- a/meshmode/interop/fiat/simplex_cell.py +++ b/meshmode/interop/fiat/simplex_cell.py @@ -105,30 +105,18 @@ class FIATSimplexCellImporter(ExternalImportHandler): self._mat, self._shift = get_affine_mapping(reference_vertices, self._unit_vertices) - def make_points(self, dim, entity_id, order): + def affinely_map_firedrake_to_meshmode(self, points): """ - Constructs a lattice of points on the *entity_id*th facet - of dimension *dim*. - - Args are exactly as in - :meth:`fiat.FIAT.reference_element.Cell.make_points` - (see `FIAT docs `_), - but the unit nodes are (affinely) mapped to :mod:`modepy` + Map points on the firedrake reference simplex to + :mod:`modepy` `unit coordinates `_. - :arg dim: Dimension of the facet we are constructing points on. - :arg entity_id: identifier to determine which facet of - dimension *dim* to construct the points on. - :arg order: Number of points to include in each direction. - - :return: an *np.array* of shape *(dim, npoints)* holding the - coordinates of each of the ver + :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 """ - points = self.data.make_points(dim, entity_id, order) - if not points: - return points - points = np.array(points) - # Points is (npoints, dim) so have to transpose - return (np.matmul(self._mat, points.T) + self._shift[:, np.newaxis]).T + return np.matmul(self._mat, points) + self._shift[:, np.newaxis] # }}} diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index 52c0e6d2..9bd03f27 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -146,7 +146,6 @@ def reordering_array(mesh_importer, key, fspace_data): (order.shape[0]//nunit_nodes, nunit_nodes) + order.shape[1:]) flip_mat = finat_element_importer.flip_matrix() - print("ORIENT=", mesh_importer._orient) reorder_nodes(mesh_importer.orientations(), new_order, flip_mat, unflip=firedrake_to_meshmode) new_order = new_order.flatten() diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py index 8ce88f26..da965cbe 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -180,6 +180,14 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): 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() -- GitLab From d2e843a0d3f8fc4a04b63378fbb52ac70f5b3035 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 23 Jun 2020 14:25:52 -0500 Subject: [PATCH 32/33] Some documentation fixes --- .../interop/firedrake/function_space_shared_data.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index 9bd03f27..c7b316a7 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -55,16 +55,16 @@ def cached(f, mesh_importer, key, *args, **kwargs): def reorder_nodes(orient, nodes, flip_matrix, unflip=False): """ - :param orient: An array of shape (nelements) of orientations, + 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 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 - - flips :param:`nodes` """ # reorder nodes (Code adapted from # meshmode.mesh.processing.flip_simplex_element_group) @@ -106,7 +106,7 @@ def reorder_nodes(orient, nodes, flip_matrix, unflip=False): @cached def reordering_array(mesh_importer, key, fspace_data): """ - :param key: A tuple (finat_element_anlog, firedrake_to_meshmode) + :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 -- GitLab From 7dad1fcbe5ac3a9ff5bda350384e093313a5ba31 Mon Sep 17 00:00:00 2001 From: benSepanski Date: Tue, 23 Jun 2020 14:26:17 -0500 Subject: [PATCH 33/33] Removed unused/unsafe __eq__/__hash__ from handlers --- meshmode/interop/__init__.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py index 1421458c..2c1bf072 100644 --- a/meshmode/interop/__init__.py +++ b/meshmode/interop/__init__.py @@ -46,17 +46,6 @@ class ExternalDataHandler: def __init__(self, data): self.data = data - def __hash__(self): - return hash((type(self), self.data)) - - def __eq__(self, other): - return isinstance(other, type(self)) and \ - isinstance(self, type(other)) and \ - self.data == other.data - - def __neq__(self, other): - return not self.__eq__(other) - # }}} -- GitLab