diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index fd2920800884fb605761f5390a4a6aa733af9697..94eb5faef2e01194d822ed61ff31f7d4bd466fbe 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -33,6 +33,69 @@ jobs: curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh . ./build-and-test-py-project-within-miniconda.sh + firedrake: + name: Pytest Firedrake + runs-on: ubuntu-latest + container: + image: 'firedrakeproject/firedrake' + steps: + - uses: actions/checkout@v1 + # This Dependencies setup is the same setup used for running firedrake + # examples, but had to be copied since github doesn't support yaml + # anchors. Changes made here should also be made there + - name: "Dependencies" + run: | + sudo apt update + sudo apt upgrade -y + sudo apt install time + sudo apt install -y pocl-opencl-icd ocl-icd-opencl-dev + sudo chown -R $(whoami) /github/home + + . /home/firedrake/firedrake/bin/activate + grep -v loopy requirements.txt > /tmp/myreq.txt + pip install -r /tmp/myreq.txt + pip install --force-reinstall git+https://github.com/benSepanski/loopy.git@firedrake-usable_for_potentials + # The Firedrake container is based on Py3.6 as of 2020-10-10, which + # doesn't have dataclasses. + pip install pytest dataclasses + pip install . + - name: "Test" + run: | + . /home/firedrake/firedrake/bin/activate + cd test + python -m pytest --tb=native -rxsw test_firedrake_interop.py + + firedrake_examples: + name: Examples Firedrake + runs-on: ubuntu-latest + container: + image: 'firedrakeproject/firedrake' + steps: + - uses: actions/checkout@v1 + # This Dependencies setup is the same setup used for running firedrake + # tests, but had to be copied since github doesn't support yaml + # anchors. Changes made here should also be made there + - name: "Dependencies" + run: | + sudo apt update + sudo apt upgrade -y + sudo apt install time + sudo apt install -y pocl-opencl-icd ocl-icd-opencl-dev + sudo chown -R $(whoami) /github/home + + . /home/firedrake/firedrake/bin/activate + grep -v loopy requirements.txt > /tmp/myreq.txt + pip install -r /tmp/myreq.txt + pip install --force-reinstall git+https://github.com/benSepanski/loopy.git@firedrake-usable_for_potentials + # The Firedrake container is based on Py3.6 as of 2020-10-10, which + # doesn't have dataclasses. + pip install pytest dataclasses + pip install . + - name: "Examples" + run: | + . /home/firedrake/firedrake/bin/activate + . ./.github/workflows/run_firedrake_examples.sh + examples3: name: Examples Conda Py3 runs-on: ubuntu-latest diff --git a/.github/workflows/run_firedrake_examples.sh b/.github/workflows/run_firedrake_examples.sh new file mode 100644 index 0000000000000000000000000000000000000000..d5f59f7aa6ad63601c6643da07dc6d9013388376 --- /dev/null +++ b/.github/workflows/run_firedrake_examples.sh @@ -0,0 +1,9 @@ +cd examples +for i in $(find . -name '*firedrake*.py' -exec grep -q __main__ '{}' \; -print ); do + echo "-----------------------------------------------------------------------" + echo "RUNNING $i" + echo "-----------------------------------------------------------------------" + dn=$(dirname "$i") + bn=$(basename "$i") + (cd $dn; time python3 "$bn") +done diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index aa2c970918b911a5273d4ce4a5b38f76f648febb..715fded35dd1e8d7a417deb0c0afdeb33707acb3 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -49,6 +49,43 @@ Python 3 POCL Examples: except: - tags +Python 3 POCL Firedrake: + tags: + - "docker-runner" + image: "firedrakeproject/firedrake" + script: + - &firedrake_setup | + sudo apt update + sudo apt upgrade -y + sudo apt install time + sudo apt install -y pocl-opencl-icd ocl-icd-opencl-dev + source ~/firedrake/bin/activate + grep -v loopy requirements.txt > myreq.txt + pip install -r myreq.txt + pip install --force-reinstall git+https://github.com/benSepanski/loopy.git@firedrake-usable_for_potentials + # The Firedrake container is based on Py3.6 as of 2020-10-10, which + # doesn't have dataclasses. + pip install pytest dataclasses + pip install . + # run tests + - cd test + - python -m pytest --tb=native --junitxml=pytest.xml -rxsw test_firedrake_interop.py + artifacts: + reports: + junit: test/pytest.xml + +Python 3 POCL Firedrake Examples: + tags: + - "docker-runner" + image: "firedrakeproject/firedrake" + script: + - *firedrake_setup + # run examples + - . ./.github/workflows/run_firedrake_examples.sh + artifacts: + reports: + junit: test/pytest.xml + Python 3 Conda: script: | CONDA_ENVIRONMENT=.test-conda-env-py3.yml diff --git a/doc/conf.py b/doc/conf.py index 738dcfc0eb42e9cef3564696ca4a828d382d6e10..cb25d8fe88defe6e882a7e6c71206e80c4e6be22 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -33,6 +33,7 @@ extensions = [ 'sphinx.ext.doctest', 'sphinx.ext.intersphinx', 'sphinx.ext.mathjax', + 'sphinx.ext.graphviz', ] # Add any paths that contain templates here, relative to this directory. @@ -281,5 +282,8 @@ intersphinx_mapping = { 'https://documen.tician.de/meshpy': None, 'https://documen.tician.de/modepy': None, 'https://documen.tician.de/loopy': None, + 'https://firedrakeproject.org/': None, 'https://tisaac.gitlab.io/recursivenodes/': None, + 'https://fenics.readthedocs.io/projects/fiat/en/latest/': None, + 'https://finat.github.io/FInAT/': None, } diff --git a/doc/index.rst b/doc/index.rst index a5e39233edd23cf2ae9fc7a3ec6757d5fffd8fe1..9caf6446ed33cfb3dcc7fa6d6d0f4782426efbbe 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -10,6 +10,7 @@ Contents: array discretization connection + interop misc Indices and tables diff --git a/doc/interop.rst b/doc/interop.rst new file mode 100644 index 0000000000000000000000000000000000000000..c56665e5032bc6d919df8ae401ed8a85bb4b22d9 --- /dev/null +++ b/doc/interop.rst @@ -0,0 +1,255 @@ +Interoperability with Other Discretization Packages +=================================================== + +Functionality in this subpackage helps import and export data to/from other +pieces of software, typically PDE solvers. + + +Firedrake +--------- + +Function Spaces/Discretizations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Users wishing to interact with :mod:`meshmode` from :mod:`firedrake` +will create a +:class:`~meshmode.interop.firedrake.connection.FiredrakeConnection` +using :func:`~meshmode.interop.firedrake.connection.build_connection_from_firedrake`, +while users wishing +to interact with :mod:`firedrake` from :mod:`meshmode` will use +will create a +:class:`~meshmode.interop.firedrake.connection.FiredrakeConnection` +using :func:`~meshmode.interop.firedrake.connection.build_connection_to_firedrake`. +the :class:`~meshmode.interop.firedrake.connection.ToFiredrakeConnection` class. +It is not recommended to create a +:class:`~meshmode.interop.firedrake.connection.FiredrakeConnection` directly. + +.. automodule:: meshmode.interop.firedrake.connection + +Meshes +^^^^^^ + +.. automodule:: meshmode.interop.firedrake.mesh + +Reference Cells +^^^^^^^^^^^^^^^ + +.. automodule:: meshmode.interop.firedrake.reference_cell + + +Implementation Details +^^^^^^^^^^^^^^^^^^^^^^ + +Converting between :mod:`firedrake` and :mod:`meshmode` is in general +straightforward. Some language is different: + +* In a mesh, a :mod:`meshmode` "element" is a :mod:`firedrake` "cell" +* A :class:`~meshmode.discretization.Discretization` is a :mod:`firedrake` + :class:`~firedrake.functionspaceimpl.WithGeometry`, usually + created by calling the function :func:`~firedrake.functionspace.FunctionSpace` + and referred to as a "function space" +* In a mesh, any vertices, faces, cells, etc. are :mod:`firedrake` + "entities" (see `the PETSc documentation on dmplex `__ + for more info on how topological mesh information is stored + in :mod:`firedrake`). + +Other than carefully tabulating how and which vertices/faces +correspond to other vertices/faces/cells, there are two main difficulties. + +1. :mod:`meshmode` has discontinuous polynomial function spaces + which may use different unit nodes than :mod:`firedrake`. +2. :mod:`meshmode` requires that all mesh elements be positively oriented, + :mod:`firedrake` does not. Meanwhile, when :mod:`firedrake` creates + a mesh, it changes the element ordering and the local vertex ordering. + +(1.) is easily handled by insisting that the :mod:`firedrake` +:class:`~firedrake.functionspaceimpl.WithGeometry` uses polynomial elements +and that the group of the :class:`~meshmode.discretization.Discretization` +being converted is a +:class:`~meshmode.discretization.poly_element.InterpolatoryQuadratureSimplexElementGroup` +of the same order. Then, on each element, the function space being +represented is the same in :mod:`firedrake` and :mod:`meshmode`. +We may simply resample to one system or another's unit nodes. + +To handle (2.), +once we associate a :mod:`meshmode` +element to the correct :mod:`firedrake` cell, we have something +like this picture: + +.. graphviz:: + + digraph{ + // created with graphviz2.38 dot + // NODES + + mmNodes [label="Meshmode\nnodes"]; + mmRef [label="Meshmode\nunit nodes"]; + fdRef [label="Firedrake\nunit nodes"]; + fdNodes [label="Firedrake\nnodes"]; + + // EDGES + + mmRef -> mmNodes [label=" f "]; + fdRef -> fdNodes [label=" g "]; + } + +(Assume we have already +ensured that :mod:`meshmode` and :mod:`firedrake` use the +same reference element by mapping :mod:`firedrake`'s reference +element onto :mod:`meshmode`'s). +If :math:`f=g`, then we can resample function values from +one node set to the other. However, if :mod:`firedrake` +has reordered the vertices or if we flipped their order to +ensure :mod:`meshmode` has positively-oriented elements, +there is some map :math:`A` applied to the reference element +which implements this permutation of barycentric coordinates. +In this case, :math:`f=g\circ A`. Now, we have a connected diagram: + +.. graphviz:: + + digraph{ + // created with graphviz2.38 dot + // NODES + + mmNodes [label="Meshmode\nnodes"]; + mmRef [label="Meshmode\nunit nodes"]; + fdRef [label="Firedrake\nunit nodes"]; + fdRef2 [label="Firedrake\nunit nodes"]; + fdNodes [label="Firedrake\nnodes"]; + + // EDGES + + {rank=same; mmRef; fdRef;} + {rank=same; mmNodes; fdNodes;} + mmRef -> fdRef [label="Resampling", dir="both"]; + mmRef -> mmNodes [label=" f "]; + fdRef -> fdRef2 [label=" A "]; + fdRef2 -> fdNodes [label=" g "]; + } + +In short, once we reorder the :mod:`firedrake` nodes so +that the mapping from the :mod:`meshmode` and :mod:`firedrake` +reference elements are the same, we can resample function values +at nodes from one set of unit nodes to another (and then undo +the reordering if converting function values +from :mod:`meshmode` to :mod:`firedrake`). The +information for this whole reordering process is +stored in +:attr:`~meshmode.interop.firedrake.connection.FiredrakeConnection.mm2fd_node_mapping`, +an array which associates each :mod:`meshmode` node +to the :mod:`firedrake` node found by tracing the +above diagram (i.e. it stores +:math:`g\circ A\circ \text{Resampling} \circ f^{-1}`). + +For Developers: Firedrake Function Space Design Crash Course +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In Firedrake, meshes and function spaces have a close relationship. +In particular, this is due to some structure described in this +`Firedrake pull request `_. +If you wish to develop on / add to the implementation of conversion +between :mod:`meshmode` and :mod:`firedrake`, you will need +to understand their design style. Below is a crash course. + +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 :class:`~firedrake.mesh.MeshTopology` + which holds information about connectivity + and other topological properties, but nothing about geometry/coordinates + etc. + +(3) A class :class:`~firedrake.functionspaceimpl.FunctionSpace` + created from a :mod:`finat` element and a + :class:`~firedrake.mesh.MeshTopology` which allows us to + define functions mapping the nodes (defined by the + :mod:`finat` element) of each element in the + :class:`~firedrake.mesh.MeshTopology` to some values. + Note that the function :func:`~firedrake.functionspace.FunctionSpace` + in the firedrake API is used to create objects of class + :class:`~firedrake.functionspaceimpl.FunctionSpace` + and :class:`~firedrake.functionspaceimpl.WithGeometry` (see + (6)). + +(4) A :class:`~firedrake.function.CoordinatelessFunction` + (in the sense that its *domain* has no coordinates) + which is a function in a + :class:`~firedrake.functionspaceimpl.FunctionSpace`. + +(5) A :class:`~firedrake.mesh.MeshGeometry` created from a + :class:`~firedrake.functionspaceimpl.FunctionSpace` + and a :class:`~firedrake.function.CoordinatelessFunction` + in that :class:`~firedrake.functionspaceimpl.FunctionSpace` + which maps each dof to its geometric coordinates. + +(6) A :class:`~firedrake.functionspaceimpl.WithGeometry` which is a + :class:`~firedrake.functionspaceimpl.FunctionSpace` together + with a :class:`~firedrake.mesh.MeshGeometry`. + This is the object returned + usually returned to the user by a call + to the :mod:`firedrake` function + :func:`~firedrake.functionspace.FunctionSpace`. + +(7) A :class:`~firedrake.function.Function` is defined on a + :class:`~firedrake.functionspaceimpl.WithGeometry`. + +Thus, by the coordinates of a mesh geometry we mean + +(a) On the hidden back-end: a :class:`~firedrake.function.CoordinatelessFunction` + *f* on some function space defined only on the mesh topology. +(b) On the front-end: A :class:`~firedrake.function.Function` + with the values of *f* but defined + on a :class:`~firedrake.functionspaceimpl.WithGeometry` + created from the :class:`~firedrake.functionspaceimpl.FunctionSpace` + *f* lives in and the :class:`~firedrake.mesh.MeshGeometry` *f* defines. + +Basically, it's this picture (where :math:`a\to b` if :math:`b` depends on :math:`a`) + +.. warning:: + + In general, the :class:`~firedrake.functionspaceimpl.FunctionSpace` + of the coordinates function + of a :class:`~firedrake.functionspaceimpl.WithGeometry` may not be the same + :class:`~firedrake.functionspaceimpl.FunctionSpace` + as for functions which live in the + :class:`~firedrake.functionspaceimpl.WithGeometry`. + This picture + only shows how the class definitions depend on each other. + + +.. graphviz:: + + digraph{ + // created with graphviz2.38 dot + // 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/examples/from_firedrake.py b/examples/from_firedrake.py new file mode 100644 index 0000000000000000000000000000000000000000..624d1dc4f03969bccf307be512418b2e33aa0bf7 --- /dev/null +++ b/examples/from_firedrake.py @@ -0,0 +1,109 @@ +__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 + + +# This example provides a brief template for bringing information in +# from firedrake and makes some plots to help you better understand +# what a FromBoundaryFiredrakeConnection does +def main(): + # If can't import firedrake, do nothing + # + # filename MUST include "firedrake" (i.e. match *firedrake*.py) in order + # to be run during CI + try: + import firedrake # noqa : F401 + except ImportError: + return 0 + + from meshmode.interop.firedrake import build_connection_from_firedrake + from firedrake import ( + UnitSquareMesh, FunctionSpace, SpatialCoordinate, Function, cos + ) + + # Create a firedrake mesh and interpolate cos(x+y) onto it + fd_mesh = UnitSquareMesh(10, 10) + fd_fspace = FunctionSpace(fd_mesh, 'DG', 2) + spatial_coord = SpatialCoordinate(fd_mesh) + fd_fntn = Function(fd_fspace).interpolate(cos(sum(spatial_coord))) + + # Make connections + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + from meshmode.array_context import PyOpenCLArrayContext + actx = PyOpenCLArrayContext(queue) + + fd_connection = build_connection_from_firedrake(actx, fd_fspace) + fd_bdy_connection = \ + build_connection_from_firedrake(actx, + fd_fspace, + restrict_to_boundary='on_boundary') + + # Plot the meshmode meshes that the connections connect to + import matplotlib.pyplot as plt + from meshmode.mesh.visualization import draw_2d_mesh + fig, (ax1, ax2) = plt.subplots(1, 2) + ax1.set_title("FiredrakeConnection") + plt.sca(ax1) + draw_2d_mesh(fd_connection.discr.mesh, + draw_vertex_numbers=False, + draw_element_numbers=False, + set_bounding_box=True) + ax2.set_title("FiredrakeConnection 'on_boundary'") + plt.sca(ax2) + draw_2d_mesh(fd_bdy_connection.discr.mesh, + draw_vertex_numbers=False, + draw_element_numbers=False, + set_bounding_box=True) + plt.show() + + # Plot fd_fntn using unrestricted FiredrakeConnection + from meshmode.discretization.visualization import make_visualizer + discr = fd_connection.discr + vis = make_visualizer(actx, discr, discr.groups[0].order+3) + field = fd_connection.from_firedrake(fd_fntn, actx=actx) + + fig = plt.figure() + ax1 = fig.add_subplot(1, 2, 1, projection='3d') + ax1.set_title("cos(x+y) in\nFiredrakeConnection") + vis.show_scalar_in_matplotlib_3d(field, do_show=False) + + # Now repeat using FiredrakeConnection restricted to 'on_boundary' + bdy_discr = fd_bdy_connection.discr + bdy_vis = make_visualizer(actx, bdy_discr, bdy_discr.groups[0].order+3) + bdy_field = fd_bdy_connection.from_firedrake(fd_fntn, actx=actx) + + ax2 = fig.add_subplot(1, 2, 2, projection='3d') + plt.sca(ax2) + ax2.set_title("cos(x+y) in\nFiredrakeConnection 'on_boundary'") + bdy_vis.show_scalar_in_matplotlib_3d(bdy_field, do_show=False) + + import matplotlib.cm as cm + fig.colorbar(cm.ScalarMappable()) + plt.show() + + +if __name__ == "__main__": + main() + +# vim: foldmethod=marker diff --git a/examples/to_firedrake.py b/examples/to_firedrake.py new file mode 100644 index 0000000000000000000000000000000000000000..30e86bcb376c0cacd2b227a912a919afdc7f2309 --- /dev/null +++ b/examples/to_firedrake.py @@ -0,0 +1,132 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +import numpy as np +import pyopencl as cl + + +# Nb: Some of the initial setup was adapted from meshmode/examplse/simple-dg.py +# written by Andreas Klockner: +# https://gitlab.tiker.net/inducer/meshmode/-/blob/7826fa5e13854bf1dae425b4226865acc10ee01f/examples/simple-dg.py # noqa : E501 +def main(): + # If can't import firedrake, do nothing + # + # filename MUST include "firedrake" (i.e. match *firedrake*.py) in order + # to be run during CI + try: + import firedrake # noqa : F401 + except ImportError: + return 0 + + # For this example, imagine we wish to solve the Laplace equation + # on a meshmode mesh with some given Dirichlet boundary conditions, + # and decide to use firedrake. + # + # To verify this is working, we use a solution to the wave equation + # to get our boundary conditions + + # {{{ First we make a discretization in meshmode and get our bcs + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + from meshmode.array_context import PyOpenCLArrayContext + actx = PyOpenCLArrayContext(queue) + + nel_1d = 16 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5, -0.5), + b=(0.5, 0.5), + n=(nel_1d, nel_1d)) + + order = 3 + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + group_factory = InterpolatoryQuadratureSimplexGroupFactory(order=order) + discr = Discretization(actx, mesh, group_factory) + + # Get our solution: we will use + # Real(e^z) = Real(e^{x+iy}) + # = e^x Real(e^{iy}) + # = e^x cos(y) + nodes = discr.nodes() + from meshmode.dof_array import thaw + for i in range(len(nodes)): + nodes[i] = thaw(actx, nodes[i]) + # First index is dimension + candidate_sol = actx.np.exp(nodes[0]) * actx.np.cos(nodes[1]) + + # }}} + + # {{{ Now send candidate_sol into firedrake and use it for boundary conditions + + from meshmode.interop.firedrake import build_connection_to_firedrake + fd_connection = build_connection_to_firedrake(discr, group_nr=0) + # convert candidate_sol to firedrake + fd_candidate_sol = fd_connection.from_meshmode(candidate_sol) + # get the firedrake function space + fd_fspace = fd_connection.firedrake_fspace() + + # set up dirichlet laplace problem in fd and solve + from firedrake import ( + FunctionSpace, TrialFunction, TestFunction, Function, inner, grad, dx, + Constant, project, DirichletBC, solve) + + # because it's easier to write down the variational problem, + # we're going to project from our "DG" space + # into a continuous one. + cfd_fspace = FunctionSpace(fd_fspace.mesh(), 'CG', order) + u = TrialFunction(cfd_fspace) + v = TestFunction(cfd_fspace) + sol = Function(cfd_fspace) + + a = inner(grad(u), grad(v)) * dx + rhs = Constant(0.0) * v * dx + bc_value = project(fd_candidate_sol, cfd_fspace) + bc = DirichletBC(cfd_fspace, bc_value, 'on_boundary') + params = {'ksp_monitor': None} + solve(a == rhs, sol, bcs=[bc], solver_parameters=params) + + # project back into our "DG" space + sol = project(sol, fd_fspace) + + # }}} + + # {{{ Take the solution from firedrake and compare it to candidate_sol + + true_sol = fd_connection.from_firedrake(sol, actx=actx) + # pull back into numpy + true_sol = actx.to_numpy(true_sol[0]) + candidate_sol = actx.to_numpy(candidate_sol[0]) + print("l^2 difference between candidate solution and firedrake solution=", + np.linalg.norm(true_sol - candidate_sol)) + + # }}} + + +if __name__ == "__main__": + main() + +# vim: foldmethod=marker diff --git a/meshmode/array_context.py b/meshmode/array_context.py index cff1690c6e1f9a91a774e224ddcf1882acd61abe..84dd81013d1d00d8342192f70cb040f6ad6b3b1e 100644 --- a/meshmode/array_context.py +++ b/meshmode/array_context.py @@ -276,9 +276,13 @@ class PyOpenCLArrayContext(ArrayContext): def call_loopy(self, program, **kwargs): program = self.transform_loopy_program(program) - if not ( - program.options.return_dict - and program.options.no_numpy): + + # accommodate loopy with and without kernel callables + try: + options = program.options + except AttributeError: + options = program.root_kernel.options + if not (options.return_dict and options.no_numpy): raise ValueError("Loopy program passed to call_loopy must " "have return_dict and no_numpy options set. " "Did you use meshmode.array_context.make_loopy_program " @@ -300,7 +304,11 @@ class PyOpenCLArrayContext(ArrayContext): def transform_loopy_program(self, program): # FIXME: This could be much smarter. import loopy as lp - all_inames = program.all_inames() + # accommodate loopy with and without kernel callables + try: + all_inames = program.all_inames() + except AttributeError: + all_inames = program.root_kernel.all_inames() inner_iname = None if "iel" not in all_inames and "i0" in all_inames: diff --git a/meshmode/interop/__init__.py b/meshmode/interop/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..ba465860a4545ee4aeeaec84c9df953efd10db76 --- /dev/null +++ b/meshmode/interop/__init__.py @@ -0,0 +1,21 @@ +__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. +""" diff --git a/meshmode/interop/firedrake/__init__.py b/meshmode/interop/firedrake/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..727e2c67c939f513b8ef2cf3293aa167982955b6 --- /dev/null +++ b/meshmode/interop/firedrake/__init__.py @@ -0,0 +1,33 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +from meshmode.interop.firedrake.connection import ( + build_connection_from_firedrake, build_connection_to_firedrake, + FiredrakeConnection) +from meshmode.interop.firedrake.mesh import ( + import_firedrake_mesh, export_mesh_to_firedrake) + +__all__ = ["build_connection_from_firedrake", "build_connection_to_firedrake", + "FiredrakeConnection", "import_firedrake_mesh", + "export_mesh_to_firedrake", + ] diff --git a/meshmode/interop/firedrake/connection.py b/meshmode/interop/firedrake/connection.py new file mode 100644 index 0000000000000000000000000000000000000000..4c61113b07ee0d10907b5df4139193f144e72688 --- /dev/null +++ b/meshmode/interop/firedrake/connection.py @@ -0,0 +1,803 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +__doc__ = """ +.. autoclass:: FiredrakeConnection +.. autofunction:: build_connection_to_firedrake +.. autofunction:: build_connection_from_firedrake +""" + +import numpy as np +import numpy.linalg as la + +from modepy import resampling_matrix + +from meshmode.interop.firedrake.mesh import ( + import_firedrake_mesh, export_mesh_to_firedrake) +from meshmode.interop.firedrake.reference_cell import ( + get_affine_reference_simplex_mapping, get_finat_element_unit_nodes) + +from meshmode.mesh.processing import get_simplex_element_flip_matrix + +from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendGroupFactory, + PolynomialRecursiveNodesGroupFactory, + ElementGroupFactory) +from meshmode.discretization import ( + Discretization, InterpolatoryElementGroupBase) + +from pytools import memoize_method + + +def _reorder_nodes(orient, nodes, flip_matrix, unflip=False): + """ + Return a flipped copy of *nodes* according to *orient*. + + :arg orient: An array of shape *(nelements)* of orientations, + >0 for positive, <0 for negative + :arg nodes: a *(nelements, nunit_nodes)* shaped array of nodes + :arg flip_matrix: The matrix used to flip each negatively-oriented + element + :arg unflip: If *True*, use transpose of *flip_matrix* to + flip negatively-oriented elements + """ + # reorder nodes (Code adapted from + # meshmode.mesh.processing.flip_simplex_element_group) + + # ( round to int bc applying on integers) + flip_mat = np.rint(flip_matrix) + if unflip: + flip_mat = flip_mat.T + + # flipping twice should be identity + assert la.norm( + np.dot(flip_mat, flip_mat) + - np.eye(len(flip_mat))) < 1e-13 + + # flip nodes that need to be flipped + flipped_nodes = np.copy(nodes) + flipped_nodes[orient < 0] = np.einsum( + "ij,ej->ei", + flip_mat, nodes[orient < 0]) + + return flipped_nodes + + +# {{{ Connection between a fd function space and mm discretization + +class FiredrakeConnection: + r""" + A connection between one group of + a meshmode discretization and a firedrake "DG" + function space. + + Users should instantiate this using + :func:`build_connection_to_firedrake` + or :func:`build_connection_from_firedrake` + + .. attribute:: discr + + A :class:`meshmode.discretization.Discretization`. + + .. attribute:: group_nr + + The group number identifying which element group of + :attr:`discr` is being connected to a firedrake function space + + .. attribute:: mm2fd_node_mapping + + Letting *element_grp = self.discr.groups[self.group_nr]*, + *mm2fd_node_mapping* is a numpy array of shape + *(element_grp.nelements, element_grp.nunit_dofs)* + whose *(i, j)*\ th entry is the :mod:`firedrake` node + index associated to the *j*\ th degree of freedom of the + *i*\ th element in *element_grp*. + + :attr:`mm2fd_node_mapping` must encode an embedding + into the :mod:`firedrake` mesh, i.e. no two :mod:`meshmode` nodes + may be associated to the same :mod:`firedrake` node + + Degrees of freedom should be associated so that + the implicit mapping from a reference element to element *i* + which maps meshmode unit dofs *0,1,...,n-1* to actual + dofs *(i, 0), (i, 1), ..., (i, n-1)* + is the same mapping which maps firedrake unit dofs + *0,1,...,n-1* to firedrake dofs + *mm2fd_node_mapping[i,0], mm2fd_node_mapping[i,1],..., + mm2fd_node_mapping[i,n-1]*. + + (See :mod:`meshmode.interop.firedrake.reference_cell` to + obtain firedrake unit dofs on the meshmode reference cell) + + .. automethod:: __init__ + """ + def __init__(self, discr, fdrake_fspace, mm2fd_node_mapping, group_nr=None): + """ + :param discr: A :class:`meshmode.discretization.Discretization` + :param fdrake_fspace: A + :class:`firedrake.functionspaceimpl.WithGeometry`. + Must use ufl family ``'Discontinuous Lagrange'``. + :param mm2fd_node_mapping: Used as attribute :attr:`mm2fd_node_mapping`. + A 2-D numpy integer array with the same dtype as + ``fdrake_fspace.cell_node_list.dtype`` + :param group_nr: The index of the group in *discr* which is + being connected to *fdrake_fspace*. The group must be a + :class:`~meshmode.discretization.InterpolatoryElementGroupBase` + of the same topological dimension as *fdrake_fspace*. + If *discr* has only one group, *group_nr=None* may be supplied. + + :raises TypeError: If any input arguments are of the wrong type, + if the designated group is of the wrong type, + or if *fdrake_fspace* is of the wrong family. + :raises ValueError: If + *mm2fd_node_mapping* is of the wrong shape + or dtype, if *group_nr* is an invalid index, or + if *group_nr* is *None* when *discr* has more than one group. + """ + # {{{ Validate input + if not isinstance(discr, Discretization): + raise TypeError("'discr' must be of type " + "meshmode.discretization.Discretization, " + "not '%s'`." % type(discr)) + from firedrake.functionspaceimpl import WithGeometry + if not isinstance(fdrake_fspace, WithGeometry): + raise TypeError("'fdrake_fspace' must be of type " + "firedrake.functionspaceimpl.WithGeometry, " + "not '%s'." % type(fdrake_fspace)) + if not isinstance(mm2fd_node_mapping, np.ndarray): + raise TypeError("'mm2fd_node_mapping' must be of type " + "numpy.ndarray, " + "not '%s'." % type(mm2fd_node_mapping)) + if not isinstance(group_nr, int) and group_nr is not None: + raise TypeError("'group_nr' must be of type int or be " + "*None*, not of type '%s'." % type(group_nr)) + # Convert group_nr to an integer if *None* + if group_nr is None: + if len(discr.groups) != 1: + raise ValueError("'group_nr' is *None* but 'discr' " + "has '%s' != 1 groups." % len(discr.groups)) + group_nr = 0 + # store element_grp as variable for convenience + element_grp = discr.groups[group_nr] + + if group_nr < 0 or group_nr >= len(discr.groups): + raise ValueError("'group_nr' has value '%s', which an invalid " + "index into list 'discr.groups' of length '%s'." + % (group_nr, len(discr.groups))) + if not isinstance(element_grp, InterpolatoryElementGroupBase): + raise TypeError("'discr.groups[group_nr]' must be of type " + "InterpolatoryElementGroupBase" + ", not '%s'." % type(element_grp)) + if fdrake_fspace.ufl_element().family() != 'Discontinuous Lagrange': + raise TypeError("'fdrake_fspace.ufl_element().family()' must be" + "'Discontinuous Lagrange', not " + f"'{fdrake_fspace.ufl_element().family()}'") + if mm2fd_node_mapping.shape != (element_grp.nelements, + element_grp.nunit_dofs): + raise ValueError("'mm2fd_node_mapping' must be of shape ", + "(%s,), not '%s'" + % ((discr.groups[group_nr].ndofs,), + mm2fd_node_mapping.shape)) + if mm2fd_node_mapping.dtype != fdrake_fspace.cell_node_list.dtype: + raise ValueError("'mm2fd_node_mapping' must have dtype " + "%s, not '%s'" % (fdrake_fspace.cell_node_list.dtype, + mm2fd_node_mapping.dtype)) + if np.size(np.unique(mm2fd_node_mapping)) != np.size(mm2fd_node_mapping): + raise ValueError("'mm2fd_node_mapping' must have unique entries; " + "no two meshmode nodes may be associated to the " + "same Firedrake node") + # }}} + + # Get meshmode unit nodes + mm_unit_nodes = element_grp.unit_nodes + # get firedrake unit nodes and map onto meshmode reference element + tdim = fdrake_fspace.mesh().topological_dimension() + fd_ref_cell_to_mm = get_affine_reference_simplex_mapping(tdim, True) + fd_unit_nodes = get_finat_element_unit_nodes(fdrake_fspace.finat_element) + fd_unit_nodes = fd_ref_cell_to_mm(fd_unit_nodes) + + # compute and store resampling matrices + self._resampling_mat_fd2mm = resampling_matrix(element_grp.basis(), + new_nodes=mm_unit_nodes, + old_nodes=fd_unit_nodes) + self._resampling_mat_mm2fd = resampling_matrix(element_grp.basis(), + new_nodes=fd_unit_nodes, + old_nodes=mm_unit_nodes) + + # Store input + self.discr = discr + self.group_nr = group_nr + self.mm2fd_node_mapping = mm2fd_node_mapping + self._mesh_geometry = fdrake_fspace.mesh() + self._ufl_element = fdrake_fspace.ufl_element() + + @memoize_method + def firedrake_fspace(self, shape=None): + """ + Return a firedrake function space that + *self.discr.groups[self.group_nr]* is connected to + of the appropriate vector dimension + + :arg shape: Either *None*, in which case a function space which maps + to scalar values is returned, a positive integer *n*, + in which case a function space which maps into *\\R^n* + is returned, or a tuple of integers defining + the shape of values in a tensor function space, + in which case a tensor function space is returned + :return: A :class:`firedrake.functionspaceimpl.WithGeometry` which + corresponds to *self.discr.groups[self.group_nr]* of the appropriate + vector dimension + + :raises TypeError: If *shape* is of the wrong type + """ + if shape is None: + from firedrake import FunctionSpace + return FunctionSpace(self._mesh_geometry, + self._ufl_element.family(), + degree=self._ufl_element.degree()) + elif isinstance(shape, int): + from firedrake import VectorFunctionSpace + return VectorFunctionSpace(self._mesh_geometry, + self._ufl_element.family(), + degree=self._ufl_element.degree(), + dim=shape) + elif isinstance(shape, tuple): + from firedrake import TensorFunctionSpace + return TensorFunctionSpace(self._mesh_geometry, + self._ufl_element.family(), + degree=self._ufl_element.degree(), + shape=shape) + else: + raise TypeError("'shape' must be *None*, an integer, " + " or a tuple of integers, not of type '%s'." + % type(shape)) + + def _validate_function(self, function, function_name, + shape=None, dtype=None): + """ + Handy helper function to validate that a firedrake function + is convertible (or can be the recipient of a conversion). + Raises error messages if wrong types, shape, dtype found + etc. + """ + # Validate that *function* is convertible + from firedrake.function import Function + if not isinstance(function, Function): + raise TypeError(f"'{function_name} must be a firedrake Function" + f" but is of unexpected type '{type(function)}'") + ufl_elt = function.function_space().ufl_element() + if ufl_elt.family() != self._ufl_element.family(): + raise ValueError(f"'{function_name}.function_space().ufl_element()" + f".family()' must be {self._ufl_element.family()}" + f", not '{type(ufl_elt.family())}'") + if ufl_elt.degree() != self._ufl_element.degree(): + raise ValueError(f"'{function_name}.function_space().ufl_element()" + f".degree()' must be {self._ufl_element.degree()}" + f", not '{ufl_elt.degree()}'") + if function.function_space().mesh() is not self._mesh_geometry: + raise ValueError(f"'{function_name}.function_space().mesh()' must" + " be the same mesh used by this connection") + if dtype is not None and function.dat.data.dtype != dtype: + raise ValueError(f"'{function_name}.dat.dtype' must be " + f"{dtype}, not '{function.dat.data.dtype}'") + if shape is not None and function.function_space().shape != shape: + raise ValueError("'{function_name}.function_space().shape' must be" + " {shape}, not '{function.function_space().shape}" + "'") + + def _validate_field(self, field, field_name, shape=None, dtype=None): + """ + Handy helper function to validate that a meshmode field + is convertible (or can be the recipient of a conversion). + Raises error messages if wrong types, shapes, dtypes found + etc. + """ + from meshmode.dof_array import DOFArray + element_grp = self.discr.groups[self.group_nr] + group_shape = (element_grp.nelements, element_grp.nunit_dofs) + + def check_dof_array(arr, arr_name): + if not isinstance(arr, DOFArray): + raise TypeError(f"'{arr_name}' must be of type " + f"meshmode.dof_array.DOFArray, not " + f"{type(arr)}") + if arr.array_context is None: + raise ValueError(f"'{arr_name}' must have a non-*None* " + "array_context") + if arr.shape != tuple([len(self.discr.groups)]): + raise ValueError(f"'{arr_name}' shape must be " + f"{tuple([len(self.discr.groups)])}, not " + f"'{arr.shape}'") + if arr[self.group_nr].shape != group_shape: + raise ValueError(f"'{arr_name}[{self.group_nr}].shape' must be" + f" {group_shape}, not " + f"'{arr[self.group_nr].shape}'") + if dtype is not None and arr.entry_dtype != dtype: + raise ValueError(f"'{arr_name}.entry_dtype' must be {dtype}," + f" not '{arr.entry_dtype}'") + + if isinstance(field, DOFArray): + if shape is not None and shape != (): + raise ValueError("shape != () and '%s' is of type DOFArray" + " instead of np.ndarray." % field_name) + check_dof_array(field, field_name) + elif isinstance(field, np.ndarray) and field.dtype == np.object: + if shape is not None and field.shape != shape: + raise ValueError(f"'{field_name}.shape' must be {shape}, not " + f"'{field.shape}'") + for multi_index, arr in np.ndenumerate(field): + arr_name = "%s[%s]" % (field_name, multi_index) + try: + check_dof_array(arr, arr_name) + except TypeError as type_err: + msg = type_err.args[0] + prefix = f"'{field_name}' is a numpy array of shape " \ + f"{field.shape}, which is interpreted as a mapping" \ + f" into a space of sahpe {field.shape}. For each " \ + r" multi-index *mi*, the *mi*\ th coordinate values " \ + f" of '{field_name}' should be represented as a " \ + f"DOFArray stored in '{field_name}[mi]'. If you are " \ + " not trying to represent a mapping into a space of " \ + f" shape {field.shape}, look at the documentation " \ + " for FiredrakeConnection.from_meshmode or " \ + "FiredrakeConnection.from_firedrake to see how " \ + "fields in a discretization are represented." + raise TypeError(prefix + "\n" + msg) + else: + raise TypeError("'field' must be of type DOFArray or a numpy object " + "array of those, not '%s'." % type(field)) + + def from_firedrake(self, function, out=None, actx=None): + """ + Transport firedrake function onto :attr:`discr` + + :arg function: A :class:`firedrake.function.Function` to transfer onto + :attr:`discr`. Its function space must have + the same family, degree, and mesh as ``self.from_fspace()``. + :arg out: Either + + 1. A :class:`~meshmode.dof_array.DOFArray` + 2. A :class:`numpy.ndarray` object array, each of whose + entries is a :class:`~meshmode.dof_array.DOFArray` + 3. *None* + + In the case of (1.), *function* must be in a + scalar function space + (i.e. `function.function_space().shape == (,)`). + In the case of (2.), the shape of *out* must match + `function.function_space().shape`. + + In either case, each :class:`~meshmode.dof_array.DOFArray` + must be a :class:`~meshmode.dof_array.DOFArray` + defined on :attr:`discr` (as described in + the documentation for :class:`~meshmode.dof_array.DOFArray`). + Also, each :class:`~meshmode.dof_array.DOFArray`'s *entry_dtype* must + match the *function.dat.data.dtype*, and be of shape + *(nelements, nunit_dofs)*. + + In case (3.), an array is created satisfying + the above requirements. + + The data in *function* is transported to :attr:`discr` + and stored in *out*, which is then returned. + :arg actx: + * If *out* is *None*, then *actx* is a + :class:`~meshmode.array_context.ArrayContext` on which + to create the :class:`~meshmode.dof_array.DOFArray` + * If *out* is not *None*, *actx* must be *None* or *out*'s + :attr:`~meshmode.dof_array.DOFArray.array_context`. + + :return: *out*, with the converted data from *function* stored in it + """ + self._validate_function(function, "function") + # get function data and shape of values + function_data = function.dat.data + fspace_shape = function.function_space().shape + + # Handle :arg:`out` + if out is not None: + self._validate_field(out, "out", fspace_shape, function_data.dtype) + # If out is supplied, check type, shape, and dtype + if actx not in (None, out.array_context): + raise ValueError("If 'out' is not *None*, 'actx' must be" + " *None* or 'out.array_context'") + else: + # If 'out' is not supplied, create it + from meshmode.array_context import ArrayContext + if not isinstance(actx, ArrayContext): + raise TypeError("If 'out' is *None*, 'actx' must be of type " + "ArrayContext, not '%s'." % type(actx)) + if fspace_shape == (): + out = self.discr.empty(actx, dtype=function_data.dtype) + else: + out = np.ndarray(fspace_shape, dtype=np.object) + for multi_index in np.ndindex(fspace_shape): + out[multi_index] = \ + self.discr.empty(actx, dtype=function_data.dtype) + + def reorder_and_resample(dof_array, fd_data): + # put the firedrake data in meshmode order and then resample, + # storing in dof_array + dof_array[self.group_nr].set( + np.einsum("ij,kj->ik", + fd_data[self.mm2fd_node_mapping], + self._resampling_mat_fd2mm) + ) + + # If scalar, just reorder and resample out + if fspace_shape == (): + reorder_and_resample(out, function_data) + else: + # firedrake drops extra dimensions + if len(function_data.shape) != 1 + len(fspace_shape): + shape = (function_data.shape[0],) + fspace_shape + function_data = function_data.reshape(shape) + # otherwise, have to grab each dofarray and the corresponding + # data from *function_data* + for multi_index in np.ndindex(fspace_shape): + dof_array = out[multi_index] + index = (np.s_[:],) + multi_index + fd_data = function_data[index] + reorder_and_resample(dof_array, fd_data) + + return out + + def from_meshmode(self, mm_field, out=None): + r""" + Transport meshmode field from :attr:`discr` into an + appropriate firedrake function space. + + If *out* is *None*, values at any firedrake + nodes associated to NO meshmode nodes are zeroed out. + If *out* is supplied, values at nodes associated to NO meshmode nodes + are not modified. + + :arg mm_field: Either + + * A :class:`~meshmode.dof_array.DOFArray` representing + a field of shape *tuple()* on :attr:`discr` + * A :class:`numpy.ndarray` of dtype "object" with + entries of class :class:`~meshmode.dof_array.DOFArray` + representing a field of shape *mm_field.shape* + on :attr:`discr` + + See :class:`~meshmode.dof.DOFArray` for further requirements. + The :attr:`group_nr` entry of each + :class:`~meshmode.dof_array.DOFArray` + must be of shape *(nelements, nunit_dofs)* and + the *element_dtype* must match that used for + :class:`firedrake.function.Function`\ s + + :arg out: If *None* then ignored, otherwise a + :class:`firedrake.function.Function` + of the right function space for the transported data + to be stored in. The shape of its function space must + match the shape of *mm_field* + + :return: a :class:`firedrake.function.Function` holding the transported + data (*out*, if *out* was not *None*) + """ + # All firedrake functions are the same dtype + dtype = self.firedrake_fspace().mesh().coordinates.dat.data.dtype + self._validate_field(mm_field, "mm_field", dtype=dtype) + + # get the shape of mm_field + from meshmode.dof_array import DOFArray + if not isinstance(mm_field, DOFArray): + fspace_shape = mm_field.shape + else: + fspace_shape = () + + # make sure out is a firedrake function in an appropriate + # function space + if out is not None: + self._validate_function(out, "out", fspace_shape, dtype) + else: + from firedrake.function import Function + # Translate shape so that don't always get a TensorFunctionSpace, + # but instead get FunctionSpace or VectorFunctionSpace when + # reasonable + shape = fspace_shape + if shape == (): + shape = None + elif len(shape) == 1: + shape = shape[0] + # make a function filled with zeros + out = Function(self.firedrake_fspace(shape)) + out.dat.data[:] = 0.0 + + out_data = out.dat.data + # Handle firedrake dropping dimensions + if len(out.dat.data.shape) != 1 + len(fspace_shape): + shape = (out.dat.data.shape[0],) + fspace_shape + out_data = out_data.reshape(shape) + + def resample_and_reorder(fd_data, dof_array): + # pull data into numpy + dof_np = dof_array.array_context.to_numpy(dof_array[self.group_nr]) + # resample the data and store in firedrake ordering + # store resampled data in firedrake ordering + fd_data[self.mm2fd_node_mapping] = \ + np.einsum("ij,kj->ik", dof_np, self._resampling_mat_mm2fd) + + # If scalar, just reorder and resample out + if fspace_shape == (): + resample_and_reorder(out_data, mm_field) + else: + # otherwise, have to grab each dofarray and the corresponding + # data from *function_data* + for multi_index in np.ndindex(fspace_shape): + # have to be careful to take view and not copy + index = (np.s_[:],) + multi_index + fd_data = out_data[index] + dof_array = mm_field[multi_index] + resample_and_reorder(fd_data, dof_array) + + return out + +# }}} + + +# {{{ Create connection from firedrake into meshmode + + +def _get_cells_to_use(fdrake_mesh, bdy_id): + """ + Return the cell indices of 'fdrake_mesh' which have at least one vertex + coinciding with a facet which is marked with firedrake marker + 'bdy_id'. + + If 'bdy_id' is *None*, returns *None* + + Separated into a function for testing purposes + + :param fdrake_mesh: A mesh as in + :func:`~meshmode.interop.firedrake.mesh.import_firedrake_mesh` + :param bdy_id: As the argument 'restrict_to_boundary' in + :func:`build_connection_from_firedrake` + """ + if bdy_id is None: + return None + + cfspace = fdrake_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) + + from pyop2.datatypes import IntType + return np.nonzero(cell_is_near_bdy)[0].astype(IntType) + + +def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None, + restrict_to_boundary=None): + + """ + Create a :class:`FiredrakeConnection` from a :mod:`firedrake` + ``"DG"`` function space by creates a corresponding + meshmode discretization and facilitating + transfer of functions to and from :mod:`firedrake`. + + :arg actx: A :class:`~meshmode.array_context.ArrayContext` + used to instantiate :attr:`FiredrakeConnection.discr`. + :arg fdrake_fspace: A :mod:`firedrake` ``"DG"`` + function space (of class + :class:`~firedrake.functionspaceimpl.WithGeometry`) built on + a mesh which is importable by + :func:`~meshmode.interop.firedrake.mesh.import_firedrake_mesh`. + :arg grp_factory: (optional) If not *None*, should be + a :class:`~meshmode.discretization.poly_element.ElementGroupFactory` + whose group class is a subclass of + :class:`~meshmode.discretization.InterpolatoryElementGroupBase`. + If *None*, and :mod:`recursivenodes` can be imported, + a :class:`~meshmode.discretization.poly_element.\ +PolynomialRecursiveNodesGroupFactory` with ``'lgl'`` nodes is used. + Note that :mod:`recursivenodes` may not be importable + as it uses :func:`math.comb`, which is new in Python 3.8. + In the case that :mod:`recursivenodes` cannot be successfully + imported, a :class:`~meshmode.discretization.poly_element.\ +PolynomialWarpAndBlendGroupFactory` is used. + :arg restrict_to_boundary: (optional) + If not *None*, then must be a valid boundary marker for + ``fdrake_fspace.mesh()``. In this case, creates a + :class:`~meshmode.discretization.Discretization` on a submesh + of ``fdrake_fspace.mesh()`` created from the cells with at least + one vertex on a facet marked with the marker + *restrict_to_boundary*. + """ + # Ensure fdrake_fspace is a function space with appropriate reference + # element. + from firedrake.functionspaceimpl import WithGeometry + if not isinstance(fdrake_fspace, WithGeometry): + raise TypeError("'fdrake_fspace' must be of firedrake type " + "WithGeometry, not '%s'." + % type(fdrake_fspace)) + ufl_elt = fdrake_fspace.ufl_element() + + if ufl_elt.family() != 'Discontinuous Lagrange': + raise ValueError("the 'fdrake_fspace.ufl_element().family()' of " + "must be be " + "'Discontinuous Lagrange', not '%s'." + % ufl_elt.family()) + # Make sure grp_factory is the right type if provided, and + # uses an interpolatory class. + if grp_factory is not None: + if not isinstance(grp_factory, ElementGroupFactory): + raise TypeError("'grp_factory' must inherit from " + "meshmode.discretization.ElementGroupFactory," + "but is instead of type " + "'%s'." % type(grp_factory)) + if not issubclass(grp_factory.group_class, + InterpolatoryElementGroupBase): + raise TypeError("'grp_factory.group_class' must inherit from" + "meshmode.discretization." + "InterpolatoryElementGroupBase, but" + " is instead of type '%s'" + % type(grp_factory.group_class)) + # If not provided, make one + else: + degree = ufl_elt.degree() + try: + # recursivenodes is only importable in Python 3.8 since + # it uses :func:`math.comb`, so need to check if it can + # be imported + import recursivenodes # noqa : F401 + family = 'lgl' # L-G-Legendre + grp_factory = PolynomialRecursiveNodesGroupFactory(degree, family) + except ImportError: + # If cannot be imported, uses warp-and-blend nodes + grp_factory = PolynomialWarpAndBlendGroupFactory(degree) + if restrict_to_boundary is not None: + uniq_markers = fdrake_fspace.mesh().exterior_facets.unique_markers + allowable_bdy_ids = list(uniq_markers) + ["on_boundary"] + if restrict_to_boundary not in allowable_bdy_ids: + raise ValueError("'restrict_to_boundary' must be one of" + " the following allowable boundary ids: " + f"{allowable_bdy_ids}, not " + f"'{restrict_to_boundary}'") + + # If only converting a portion of the mesh near the boundary, get + # *cells_to_use* as described in + # :func:`meshmode.interop.firedrake.mesh.import_firedrake_mesh` + cells_to_use = _get_cells_to_use(fdrake_fspace.mesh(), + restrict_to_boundary) + + # Create to_discr + mm_mesh, orient = import_firedrake_mesh(fdrake_fspace.mesh(), + cells_to_use=cells_to_use) + to_discr = Discretization(actx, mm_mesh, grp_factory) + + # get firedrake unit nodes and map onto meshmode reference element + group = to_discr.groups[0] + fd_ref_cell_to_mm = get_affine_reference_simplex_mapping(group.dim, + True) + fd_unit_nodes = get_finat_element_unit_nodes(fdrake_fspace.finat_element) + fd_unit_nodes = fd_ref_cell_to_mm(fd_unit_nodes) + # Flipping negative elements corresponds to reordering the nodes. + # We handle reordering by storing the permutation explicitly as + # a numpy array + + # Get the reordering fd->mm. + flip_mat = get_simplex_element_flip_matrix(ufl_elt.degree(), + fd_unit_nodes) + fd_cell_node_list = fdrake_fspace.cell_node_list + if cells_to_use is not None: + fd_cell_node_list = fd_cell_node_list[cells_to_use] + # flip fd_cell_node_list + flipped_cell_node_list = _reorder_nodes(orient, + fd_cell_node_list, + flip_mat, + unflip=False) + + assert np.size(np.unique(flipped_cell_node_list)) == \ + np.size(flipped_cell_node_list), \ + "A firedrake node in a 'DG' space got duplicated" + + return FiredrakeConnection(to_discr, + fdrake_fspace, + flipped_cell_node_list) + +# }}} + + +# {{{ Create connection to firedrake from meshmode + + +def build_connection_to_firedrake(discr, group_nr=None, comm=None): + """ + Create a connection from a meshmode discretization + into firedrake. Create a corresponding "DG" function + space and allow for conversion back and forth + by resampling at the nodes. + + :param discr: A :class:`~meshmode.discretization.Discretization` + to intialize the connection with + :param group_nr: The group number of the discretization to convert. + If *None* there must be only one group. The selected group + must be of type + :class:`~meshmode.discretization.poly_element.\ +InterpolatoryQuadratureSimplexElementGroup`. + + :param comm: Communicator to build a dmplex object on for the created + firedrake mesh + """ + if group_nr is None: + if len(discr.groups) != 1: + raise ValueError("'group_nr' is *None*, but 'discr' has '%s' " + "!= 1 groups." % len(discr.groups)) + group_nr = 0 + el_group = discr.groups[group_nr] + + from firedrake.functionspace import FunctionSpace + fd_mesh, fd_cell_order, perm2cells = \ + export_mesh_to_firedrake(discr.mesh, group_nr, comm) + fspace = FunctionSpace(fd_mesh, 'DG', el_group.order) + # get firedrake unit nodes and map onto meshmode reference element + dim = fspace.mesh().topological_dimension() + fd_ref_cell_to_mm = get_affine_reference_simplex_mapping(dim, True) + fd_unit_nodes = get_finat_element_unit_nodes(fspace.finat_element) + fd_unit_nodes = fd_ref_cell_to_mm(fd_unit_nodes) + + # **_cell_node holds the node nrs in shape *(ncells, nunit_nodes)* + fd_cell_node = fspace.cell_node_list + + # To get the meshmode to firedrake node assocation, we need to handle + # local vertex reordering and cell reordering. + from pyop2.datatypes import IntType + mm2fd_node_mapping = np.ndarray((el_group.nelements, el_group.nunit_dofs), + dtype=IntType) + for perm, cells in perm2cells.items(): + # reordering_arr[i] should be the fd node corresponding to meshmode + # node i + # + # The jth meshmode cell corresponds to the fd_cell_order[j]th + # firedrake cell. If *nodeperm* is the permutation of local nodes + # applied to the *j*\ th meshmode cell, the firedrake node + # fd_cell_node[fd_cell_order[j]][k] corresponds to the + # mm_cell_node[j, nodeperm[k]]th meshmode node. + # + # Note that the permutation on the unit nodes may not be the + # same as the permutation on the barycentric coordinates (*perm*). + # Importantly, the permutation is derived from getting a flip + # matrix from the Firedrake unit nodes, not necessarily the meshmode + # unit nodes + # + flip_mat = get_simplex_element_flip_matrix(el_group.order, + fd_unit_nodes, + np.argsort(perm)) + flip_mat = np.rint(flip_mat).astype(IntType) + fd_permuted_cell_node = np.matmul(fd_cell_node[fd_cell_order[cells]], + flip_mat.T) + mm2fd_node_mapping[cells] = fd_permuted_cell_node + + assert np.size(np.unique(mm2fd_node_mapping)) == \ + np.size(mm2fd_node_mapping), \ + "A firedrake node in a 'DG' space got duplicated" + return FiredrakeConnection(discr, + fspace, + mm2fd_node_mapping, + group_nr=group_nr) + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py new file mode 100644 index 0000000000000000000000000000000000000000..db2f0a2fc8998c801c9c1257dc939c6916580a5b --- /dev/null +++ b/meshmode/interop/firedrake/mesh.py @@ -0,0 +1,947 @@ +arg = "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 logging +import numpy as np + +from modepy import resampling_matrix, simplex_best_available_basis + +from meshmode.mesh import (BTAG_ALL, BTAG_REALLY_ALL, BTAG_INDUCED_BOUNDARY, + FacialAdjacencyGroup, Mesh, NodalAdjacency, SimplexElementGroup) +from meshmode.interop.firedrake.reference_cell import ( + get_affine_reference_simplex_mapping, get_finat_element_unit_nodes) + +from pytools import ProcessLogger + +__doc__ = """ +.. autofunction:: import_firedrake_mesh +.. autofunction:: export_mesh_to_firedrake +""" + + +logger = logging.getLogger(__name__) + + +# {{{ functions to extract information from Mesh Topology + + +def _get_firedrake_nodal_info(fdrake_mesh_topology, cells_to_use=None): + """ + Get nodal adjacency and vertex indices corresponding + to a firedrake mesh topology. Note that as we do not use + geometric information, there is no guarantee that elements + have a positive orientation. + + The elements (in firedrake lingo, the cells) + are guaranteed to have the same numbering in :mod:`meshmode` + as :mod:`firedrdake` + + :arg fdrake_mesh_topology: A :mod:`firedrake` instance of class + :class:`MeshTopology` or :class:`MeshGeometry`. + + :arg cells_to_use: Ignored if *None*. Otherwise, assumed to be + a numpy array holding the firedrake cell indices to use. + Any cells not on this array are ignored. + Cells used are assumed to appear exactly once in the array. + The cell's element index in the nodal adjacency will be its + index in *cells_to_use*. + + :return: Returns *vertex_indices* as a numpy array of shape + *(nelements, ref_element.nvertices)* (as described by + the ``vertex_indices`` attribute of a :class:`MeshElementGroup`) + and a :class:`NodalAdjacency` constructed from + *fdrake_mesh_topology* + as a tuple *(vertex_indices, nodal_adjacency)*. + + Even if *cells_to_use* is not *None*, the vertex indices + are still the global firedrake vertex indices. + """ + top = fdrake_mesh_topology.topology + + # If you don't understand dmplex, look at the PETSc reference + # here: https://cse.buffalo.edu/~knepley/classes/caam519/CSBook.pdf + # used to get topology info + # FIXME : not sure how to get around the private access + top_dm = top._topology_dm + + # Get range of dmplex ids for cells, facets, and vertices + f_start, f_end = top_dm.getHeightStratum(1) + v_start, v_end = top_dm.getDepthStratum(0) + + # FIXME : not sure how to get around the private accesses + # Maps dmplex vert id -> firedrake vert index + vert_id_dmp_to_fd = top._vertex_numbering.getOffset + + # We will fill in the values of vertex indices as we go + if cells_to_use is None: + num_cells = top.num_cells() + else: + num_cells = np.size(cells_to_use) + from pyop2.datatypes import IntType + vertex_indices = -np.ones((num_cells, top.ufl_cell().num_vertices()), + dtype=IntType) + # This will map fd cell ndx (or its new index as dictated by + # *cells_to_use* if *cells_to_use* + # is not *None*) + # -> list of fd cell indices which share a vertex + cell_to_nodal_neighbors = {} + # This will map dmplex facet id -> list of adjacent + # (fd cell ndx, firedrake local fac num) + facet_to_cells = {} + # This will map dmplex vert id -> list of fd cell + # indices which touch this vertex, + # Primarily used to construct cell_to_nodal_neighbors + vert_to_cells = {} + + # Loop through each cell (cell closure is all the dmplex ids for any + # verts, faces, etc. associated with the cell) + cell_closure = top.cell_closure + if cells_to_use is not None: + cell_closure = cell_closure[cells_to_use, :] + for fd_cell_ndx, closure_dmp_ids in enumerate(cell_closure): + # Store the vertex indices + dmp_verts = closure_dmp_ids[np.logical_and(v_start <= closure_dmp_ids, + closure_dmp_ids < v_end)] + vertex_indices[fd_cell_ndx][:] = np.array([ + vert_id_dmp_to_fd(dmp_vert) for dmp_vert in dmp_verts]) + + # Record this cell as touching the facet and remember its local + # facet number (the index at which it appears) + dmp_fac_ids = closure_dmp_ids[np.logical_and(f_start <= closure_dmp_ids, + closure_dmp_ids < f_end)] + for loc_fac_nr, dmp_fac_id in enumerate(dmp_fac_ids): + facet_to_cells.setdefault(dmp_fac_id, []).append((fd_cell_ndx, + loc_fac_nr)) + + # Record this vertex as touching the cell, and mark this cell + # as nodally adjacent (in cell_to_nodal_neighbors) to any + # cells already documented as touching this cell + cell_to_nodal_neighbors[fd_cell_ndx] = [] + for dmp_vert_id in dmp_verts: + vert_to_cells.setdefault(dmp_vert_id, []) + for other_cell_ndx in vert_to_cells[dmp_vert_id]: + cell_to_nodal_neighbors[fd_cell_ndx].append(other_cell_ndx) + cell_to_nodal_neighbors[other_cell_ndx].append(fd_cell_ndx) + vert_to_cells[dmp_vert_id].append(fd_cell_ndx) + + # make sure that no -1s remain in vertex_indices (i.e. none are left unset) + assert np.all(vertex_indices >= 0) + + # Next go ahead and compute nodal adjacency by creating + # neighbors and neighbor_starts as specified by :class:`NodalAdjacency` + neighbors = [] + if cells_to_use is None: + num_cells = top.num_cells() + else: + num_cells = np.size(cells_to_use) + neighbors_starts = np.zeros(num_cells + 1, dtype=IntType) + for iel in range(len(cell_to_nodal_neighbors)): + neighbors += cell_to_nodal_neighbors[iel] + neighbors_starts[iel+1] = len(neighbors) + + neighbors = np.array(neighbors, dtype=IntType) + + nodal_adjacency = NodalAdjacency(neighbors_starts=neighbors_starts, + neighbors=neighbors) + + return vertex_indices, nodal_adjacency + + +def _get_firedrake_boundary_tags(fdrake_mesh, tag_induced_boundary=False): + """ + 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`) + + :arg fdrake_mesh: A :mod:`firedrake` :class:`MeshTopology` or + :class:`MeshGeometry` + :arg tag_induced_boundary: If *True*, tag induced boundary with + :class:`~meshmode.mesh.BTAG_INDUCED_BOUNDARY` + + :return: A tuple of boundary tags + """ + bdy_tags = [BTAG_ALL, BTAG_REALLY_ALL] + if tag_induced_boundary: + bdy_tags.append(BTAG_INDUCED_BOUNDARY) + + unique_markers = fdrake_mesh.topology.exterior_facets.unique_markers + if unique_markers is not None: + bdy_tags.extend(unique_markers) + + return bdy_tags + + +def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, + cells_to_use=None): + """ + Return facial_adjacency_groups corresponding to + the given firedrake mesh topology. Note that as we do not + have geometric information, elements may need to be + flipped later. + + :arg fdrake_mesh_topology: A :mod:`firedrake` instance of class + :class:`MeshTopology` or :class:`MeshGeometry`. + :arg cells_to_use: If *None*, then this argument is ignored. + Otherwise, assumed to be a numpy array of unique firedrake + cell ids indicating which cells of the mesh to include, + as well as inducing a new cell index for those cells. + Also, in this case boundary faces are tagged + with :class:`BTAG_INDUCED_BOUNDARY` if they are not a boundary + face in *fdrake_mesh_topology* but become a boundary + because the opposite cell is not in *cells_to_use*. + Boundary faces in *fdrake_mesh_topology* are marked + with :class:`BTAG_ALL`. Both are marked with + :class:`BTAG_REALLY_ALL`. + + :return: A list of maps to :class:`FacialAdjacencyGroup`s as required + by a :mod:`meshmode` :class:`Mesh`. + """ + top = fdrake_mesh_topology.topology + # We only need one group + # for interconnectivity and one for boundary connectivity. + # The tricky part is moving from firedrake local facet numbering + # (ordered lexicographically by the vertex excluded from the face, + # search for "local facet number" in the following paper for + # a reference on this... + # https://spiral.imperial.ac.uk/bitstream/10044/1/28819/2/mlange-firedrake-dmplex-accepted.pdf # noqa : E501 + # ) + # and meshmode's facet ordering: obtained from a simplex element + # group + mm_simp_group = SimplexElementGroup(1, None, None, + dim=top.cell_dimension()) + mm_face_vertex_indices = mm_simp_group.face_vertex_indices() + # map firedrake local face number to meshmode local face number + fd_loc_fac_nr_to_mm = {} + # Figure out which vertex is excluded to get the corresponding + # firedrake local index + all_local_facet_nrs = set(range(top.ufl_cell().num_vertices())) + for mm_local_facet_nr, face in enumerate(mm_face_vertex_indices): + fd_local_facet_nr = all_local_facet_nrs - set(face) + assert len(fd_local_facet_nr) == 1 + (fd_local_facet_nr,) = fd_local_facet_nr # extract id from set({id}) + fd_loc_fac_nr_to_mm[fd_local_facet_nr] = mm_local_facet_nr + + # build a look-up table from firedrake markers to the appropriate values + # in the neighbors array for the external and internal facial adjacency + # groups + bdy_tags = _get_firedrake_boundary_tags( + top, tag_induced_boundary=cells_to_use is not None) + boundary_tag_to_index = {bdy_tag: i for i, bdy_tag in enumerate(bdy_tags)} + marker_to_neighbor_value = {} + from meshmode.mesh import _boundary_tag_bit + # None for no marker + marker_to_neighbor_value[None] = \ + -(_boundary_tag_bit(bdy_tags, boundary_tag_to_index, BTAG_REALLY_ALL) + | _boundary_tag_bit(bdy_tags, boundary_tag_to_index, BTAG_ALL)) + for marker in top.exterior_facets.unique_markers: + marker_to_neighbor_value[marker] = \ + -(_boundary_tag_bit(bdy_tags, boundary_tag_to_index, marker) + | -marker_to_neighbor_value[None]) + + # {{{ build the FacialAdjacencyGroup for internal connectivity + + # Get the firedrake cells associated to each interior facet + int_facet_cell = top.interior_facets.facet_cell + # Get the firedrake local facet numbers and map them to the + # meshmode local facet numbers + int_fac_loc_nr = top.interior_facets.local_facet_dat.data + int_fac_loc_nr = \ + np.array([[fd_loc_fac_nr_to_mm[fac_nr] for fac_nr in fac_nrs] + for fac_nrs in int_fac_loc_nr]) + # elements neighbors element_faces neighbor_faces are as required + # for a :class:`FacialAdjacencyGroup`. + + int_elements = int_facet_cell.flatten() + int_neighbors = np.concatenate((int_facet_cell[:, 1], int_facet_cell[:, 0])) + int_element_faces = int_fac_loc_nr.flatten().astype(Mesh.face_id_dtype) + int_neighbor_faces = np.concatenate((int_fac_loc_nr[:, 1], + int_fac_loc_nr[:, 0])) + int_neighbor_faces = int_neighbor_faces.astype(Mesh.face_id_dtype) + # If only using some of the cells + from pyop2.datatypes import IntType + if cells_to_use is not None: + to_keep = np.isin(int_elements, cells_to_use) + cells_to_use_inv = dict(zip(cells_to_use, + np.arange(np.size(cells_to_use), + dtype=IntType))) + + # Keep the cells that we are using and change old cell index + # to new cell index + int_elements = np.vectorize(cells_to_use_inv.__getitem__)( + int_elements[to_keep]) + int_element_faces = int_element_faces[to_keep] + int_neighbors = int_neighbors[to_keep] + int_neighbor_faces = int_neighbor_faces[to_keep] + # For neighbor cells, change to new cell index or record + # as a new boundary (if the neighbor cell is not being used) + newly_created_exterior_facs = [] + for ndx, icell in enumerate(int_neighbors): + try: + int_neighbors[ndx] = cells_to_use_inv[icell] + except KeyError: + newly_created_exterior_facs.append(ndx) + # Make boolean array: 1 if a newly created exterior facet, 0 if + # remains an interior facet + newly_created_exterior_facs = np.isin(np.arange(np.size(int_elements)), + newly_created_exterior_facs) + new_ext_elements = int_elements[newly_created_exterior_facs] + new_ext_element_faces = int_element_faces[newly_created_exterior_facs] + new_ext_neighbor_tag = -(_boundary_tag_bit(bdy_tags, + boundary_tag_to_index, + BTAG_REALLY_ALL) + | _boundary_tag_bit(bdy_tags, + boundary_tag_to_index, + BTAG_INDUCED_BOUNDARY)) + new_ext_neighbors = np.full(new_ext_elements.shape, + new_ext_neighbor_tag, + dtype=IntType) + new_ext_neighbor_faces = np.full(new_ext_elements.shape, + 0, + dtype=Mesh.face_id_dtype) + # Remove any (previously) interior facets that have become exterior + # facets + remaining_int_facs = np.logical_not(newly_created_exterior_facs) + int_elements = int_elements[remaining_int_facs] + int_element_faces = int_element_faces[remaining_int_facs] + int_neighbors = int_neighbors[remaining_int_facs] + int_neighbor_faces = int_neighbor_faces[remaining_int_facs] + + interconnectivity_grp = FacialAdjacencyGroup(igroup=0, ineighbor_group=0, + elements=int_elements, + neighbors=int_neighbors, + element_faces=int_element_faces, + neighbor_faces=int_neighbor_faces) + + # }}} + + # {{{ build the FacialAdjacencyGroup for boundary faces + + # We can get the elements directly from exterior facets + ext_elements = top.exterior_facets.facet_cell.flatten() + + ext_element_faces = np.array([fd_loc_fac_nr_to_mm[fac_nr] for fac_nr in + top.exterior_facets.local_facet_dat.data], + dtype=Mesh.face_id_dtype) + ext_neighbor_faces = np.zeros(ext_element_faces.shape, + dtype=Mesh.face_id_dtype) + # If only using some of the cells, throw away unused cells and + # move to new cell index + if cells_to_use is not None: + to_keep = np.isin(ext_elements, cells_to_use) + ext_elements = np.vectorize(cells_to_use_inv.__getitem__)( + ext_elements[to_keep]) + ext_element_faces = ext_element_faces[to_keep] + ext_neighbor_faces = ext_neighbor_faces[to_keep] + + # tag the boundary, making sure to record custom tags + # (firedrake "markers") if present + if top.exterior_facets.markers is not None: + ext_neighbors = np.zeros(ext_elements.shape, dtype=IntType) + for ifac, marker in enumerate(top.exterior_facets.markers): + ext_neighbors[ifac] = marker_to_neighbor_value[marker] + else: + ext_neighbors = np.full(ext_elements.shape, + marker_to_neighbor_value[None], + dtype=IntType) + + # If not using all the cells, some interior facets may have become + # exterior facets: + if cells_to_use is not None: + # Record any newly created exterior facets + ext_elements = np.concatenate((ext_elements, new_ext_elements)) + ext_element_faces = np.concatenate((ext_element_faces, + new_ext_element_faces)) + ext_neighbor_faces = np.concatenate((ext_neighbor_faces, + new_ext_neighbor_faces)) + ext_neighbors = np.concatenate((ext_neighbors, new_ext_neighbors)) + + exterior_grp = FacialAdjacencyGroup(igroup=0, ineighbor=None, + elements=ext_elements, + element_faces=ext_element_faces, + neighbors=ext_neighbors, + neighbor_faces=ext_neighbor_faces) + + # }}} + + return [{0: interconnectivity_grp, None: exterior_grp}] + +# }}} + + +# {{{ Orientation computation + +def _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices, + cells_to_use, + normals=None, no_normals_warn=True): + r""" + Return the orientations of the mesh elements: + + :arg fdrake_mesh: As described in :func:`import_firedrake_mesh` + :arg unflipped_group: A :class:`SimplexElementGroup` instance with + (potentially) some negatively oriented elements. + :arg vertices: The vertex coordinates as a numpy array of shape + *(ambient_dim, nvertices)* (the vertices of *unflipped_group*) + :arg normals: As described in :func:`import_firedrake_mesh` + :arg no_normals_warn: As described in :func:`import_firedrake_mesh` + :arg cells_to_use: If *None*, then ignored. Otherwise, a numpy array + of unique firedrake cell indices indicating which cells to use. + + :return: A numpy array, the *i*\ th element is > 0 if the *i*\ th element + is positively oriented, < 0 if negatively oriented. + Mesh must have co-dimension 0 or 1. If *cells_to_use* is not + *None*, then the *i*\ th entry corresponds to the + *cells_to_use[i]*\ th element. + """ + # compute orientations + tdim = fdrake_mesh.topological_dimension() + gdim = fdrake_mesh.geometric_dimension() + + orient = None + if gdim == tdim: + # If the co-dimension is 0, :mod:`meshmode` has a convenient + # function to compute cell orientations + from meshmode.mesh.processing import \ + find_volume_mesh_element_group_orientation + + orient = find_volume_mesh_element_group_orientation(vertices, + unflipped_group) + + elif tdim == 1 and gdim == 2: + # In this case we have a 1-surface embedded in 2-space. + # Firedrake does not provide any convenient way of + # letting the user set cell orientations in this case, so we + # have to ask the user for cell normals directly. + if cells_to_use is None: + num_cells = fdrake_mesh.num_cells() + else: + num_cells = np.size(cells_to_use) + orient = np.ones(num_cells) + if normals: + for i, (normal, vert_indices) in enumerate( + zip(np.array(normals), unflipped_group.vertex_indices)): + edge = vertices[:, vert_indices[1]] - vertices[:, vert_indices[0]] + if np.cross(normal, edge) < 0: + orient[i] = -1.0 + elif 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. + # In this case, we assume the user has called + # :func:`firedrake.mesh.MeshGeometry.init_cell_orientations`, see + # https://www.firedrakeproject.org/variational-problems.html#ensuring-consistent-cell-orientations # noqa : E501 + # for a tutorial on how these are usually initialized. + # + # Unfortunately, *init_cell_orientations* is currently only implemented + # in 3D, so we can't use this in the 1/2 case. + orient = fdrake_mesh.cell_orientations().dat.data + if cells_to_use is not None: + orient = orient[cells_to_use] + r""" + Convert (0 \implies negative, 1 \implies positive) to + (-1 \implies negative, 1 \implies positive) + """ + orient *= 2 + orient -= 1 + # Make sure the mesh fell into one of the cases + # Nb : This should be guaranteed by previous checks, + # but is here anyway in case of future development. + assert orient is not None + return orient + +# }}} + + +# {{{ Mesh importing from firedrake + +def import_firedrake_mesh(fdrake_mesh, cells_to_use=None, + normals=None, no_normals_warn=None): + """ + Create a :class:`meshmode.mesh.Mesh` + from a :class:`firedrake.mesh.MeshGeometry` + with the same cells/elements, vertices, nodes, + mesh order, and facial adjacency. + + The vertex and node coordinates will be the same, as well + as the cell/element ordering. However, :mod:`firedrake` + does not require elements to be positively oriented, + so any negative elements are flipped + as in :func:`meshmode.processing.flip_simplex_element_group`. + + The flipped cells/elements are identified by the returned + *firedrake_orient* array + + :arg fdrake_mesh: :class:`firedrake.mesh.MeshGeometry`. + This mesh **must** be in a space of ambient dimension + 1, 2, or 3 and have co-dimension of 0 or 1. + It must use a simplex as a reference element. + + In the case of a 2-dimensional mesh embedded in 3-space, + the method ``fdrake_mesh.init_cell_orientations`` must + have been called. + + In the case of a 1-dimensional mesh embedded in 2-space, + see parameters *normals* and *no_normals_warn*. + + Finally, the ``coordinates`` attribute must have a function + space whose *finat_element* associates a degree + of freedom with each vertex. In particular, + this means that the vertices of the mesh must have well-defined + coordinates. + For those unfamiliar with :mod:`firedrake`, you can + verify this by looking at + + .. code-block:: python + + coords_fspace = fdrake_mesh.coordinates.function_space() + vertex_entity_dofs = coords_fspace.finat_element.entity_dofs()[0] + for entity, dof_list in vertex_entity_dofs.items(): + assert len(dof_list) > 0 + + :arg cells_to_use: *cells_to_use* is primarily intended for use + internally by :func:`~meshmode.interop.firedrake.connection.\ +build_connection_from_firedrake`. + *cells_to_use* must be either + + 1. *None*, in which case this argument is ignored, or + 2. a numpy array of unique firedrake cell indexes. + + In case (2.), + only cells whose index appears in *cells_to_use* are included + in the resultant mesh, and their index in *cells_to_use* + becomes the element index in the resultant mesh element group. + Any faces or vertices which do not touch a cell in + *cells_to_use* are also ignored. + Note that in this latter case, some faces that are not + boundaries in *fdrake_mesh* may become boundaries in the + returned mesh. These "induced" boundaries are marked with + :class:`~meshmode.interop.firedrake.mesh.BTAG_INDUCED_BOUNDARY` + instead of :class:`~meshmode.mesh.BTAG_ALL`. + + :arg normals: **Only** used if *fdrake_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 *mesh.coordinate.function_space()*'s + :attr:`cell_node_list`) + + :arg no_normals_warn: If *True* (the default), raises a warning + if *fdrake_mesh* is a 1-surface embedded in 2-space + and *normals* is *None*. + + :return: A tuple *(meshmode mesh, firedrake_orient)*. + ``firedrake_orient < 0`` is *True* for any negatively + oriented firedrake cell (which was flipped by meshmode) + and False for any positively oriented firedrake cell + (which was not flipped by meshmode). + """ + # Type validation + from firedrake.mesh import MeshGeometry + if not isinstance(fdrake_mesh, MeshGeometry): + raise TypeError("'fdrake_mesh_topology' must be an instance of " + "firedrake.mesh.MeshGeometry, " + "not '%s'." % type(fdrake_mesh)) + if cells_to_use is not None: + if not isinstance(cells_to_use, np.ndarray): + raise TypeError("'cells_to_use' must be a np.ndarray or " + "*None*") + assert len(cells_to_use.shape) == 1 + assert np.size(np.unique(cells_to_use)) == np.size(cells_to_use), \ + ":arg:`cells_to_use` must have unique entries" + assert np.all(np.logical_and(cells_to_use >= 0, + cells_to_use < fdrake_mesh.num_cells())) + assert fdrake_mesh.ufl_cell().is_simplex(), "Mesh must use simplex cells" + gdim = fdrake_mesh.geometric_dimension() + tdim = fdrake_mesh.topological_dimension() + assert gdim in [1, 2, 3], "Mesh must be in space of ambient dim 1, 2, or 3" + assert gdim - tdim in [0, 1], "Mesh co-dimension must be 0 or 1" + # firedrake meshes are not guaranteed be fully instantiated until + # the .init() method is called. In particular, the coordinates function + # may not be accessible if we do not call init(). If the mesh has + # already been initialized, nothing will change. For more details + # on why we need a second initialization, see + # this pull request: + # https://github.com/firedrakeproject/firedrake/pull/627 + # which details how Firedrake implements a mesh's coordinates + # as a function on that very same mesh + fdrake_mesh.init() + + # Get all the nodal information we can from the topology + bdy_tags = _get_firedrake_boundary_tags( + fdrake_mesh, tag_induced_boundary=cells_to_use is not None) + + with ProcessLogger(logger, "Retrieving vertex indices and computing " + "NodalAdjacency from firedrake mesh"): + vertex_indices, nodal_adjacency = \ + _get_firedrake_nodal_info(fdrake_mesh, cells_to_use=cells_to_use) + + # If only using some cells, vertices may need new indices as many + # will be removed + if cells_to_use is not None: + vert_ndx_new2old = np.unique(vertex_indices.flatten()) + vert_ndx_old2new = dict(zip(vert_ndx_new2old, + np.arange(np.size(vert_ndx_new2old), + dtype=vertex_indices.dtype))) + vertex_indices = \ + np.vectorize(vert_ndx_old2new.__getitem__)(vertex_indices) + + with ProcessLogger(logger, "Building (possibly) unflipped " + "SimplexElementGroup from firedrake unit nodes/nodes"): + + # Grab the mesh reference element and cell dimension + coord_finat_elt = fdrake_mesh.coordinates.function_space().finat_element + cell_dim = fdrake_mesh.cell_dimension() + + # Get finat unit nodes and map them onto the meshmode reference simplex + finat_unit_nodes = get_finat_element_unit_nodes(coord_finat_elt) + fd_ref_to_mm = get_affine_reference_simplex_mapping(cell_dim, True) + finat_unit_nodes = fd_ref_to_mm(finat_unit_nodes) + + # Now grab the nodes + coords = fdrake_mesh.coordinates + cell_node_list = coords.function_space().cell_node_list + if cells_to_use is not None: + cell_node_list = cell_node_list[cells_to_use] + nodes = np.real(coords.dat.data[cell_node_list]) + # Add extra dim in 1D for shape (nelements, nunit_nodes, dim) + if tdim == 1: + nodes = np.reshape(nodes, nodes.shape + (1,)) + # Transpose nodes to have shape (dim, nelements, nunit_nodes) + nodes = np.transpose(nodes, (2, 0, 1)) + + # make a group (possibly with some elements that need to be flipped) + unflipped_group = SimplexElementGroup(coord_finat_elt.degree, + vertex_indices, + nodes, + dim=cell_dim, + unit_nodes=finat_unit_nodes) + + # Next get the vertices (we'll need these for the orientations) + with ProcessLogger(logger, "Obtaining vertex coordinates"): + coord_finat = fdrake_mesh.coordinates.function_space().finat_element + # unit_vertex_indices are the element-local indices of the nodes + # which coincide with the vertices, i.e. for element *i*, + # vertex 0's coordinates would be nodes[i][unit_vertex_indices[0]]. + # This assumes each vertex has some node which coincides with it... + # which is normally fine to assume for firedrake meshes. + unit_vertex_indices = [] + # iterate through the dofs associated to each vertex on the + # reference element + for _, dofs in sorted(coord_finat.entity_dofs()[0].items()): + assert len(dofs) == 1, \ + "The function space of the mesh coordinates must have" \ + " exactly one degree of freedom associated with " \ + " each vertex in order to determine vertex coordinates" + dof, = dofs + unit_vertex_indices.append(dof) + + # Now get the vertex coordinates as *(dim, nvertices)*-shaped array + if cells_to_use is not None: + nvertices = np.size(vert_ndx_new2old) + else: + nvertices = fdrake_mesh.num_vertices() + vertices = np.ndarray((gdim, nvertices), dtype=nodes.dtype) + recorded_verts = set() + for icell, cell_vertex_indices in enumerate(vertex_indices): + for local_vert_id, global_vert_id in enumerate(cell_vertex_indices): + if global_vert_id not in recorded_verts: + recorded_verts.add(global_vert_id) + local_node_nr = unit_vertex_indices[local_vert_id] + vertices[:, global_vert_id] = nodes[:, icell, local_node_nr] + + # Use the vertices to compute the orientations and flip the group + with ProcessLogger(logger, "Computing cell orientations"): + orient = _get_firedrake_orientations(fdrake_mesh, + unflipped_group, + vertices, + cells_to_use=cells_to_use, + normals=normals, + no_normals_warn=no_normals_warn) + + with ProcessLogger(logger, "Flipping group"): + from meshmode.mesh.processing import flip_simplex_element_group + group = flip_simplex_element_group(vertices, unflipped_group, orient < 0) + + # Now, any flipped element had its 0 vertex and 1 vertex exchanged. + # This changes the local facet nr, so we need to create and then + # fix our facial adjacency groups. To do that, we need to figure + # out which local facet numbers switched. + face_vertex_indices = group.face_vertex_indices() + # face indices of the faces not containing vertex 0 and not + # containing vertex 1, respectively + no_zero_face_ndx, no_one_face_ndx = None, None + for iface, face in enumerate(face_vertex_indices): + if 0 not in face: + no_zero_face_ndx = iface + elif 1 not in face: + no_one_face_ndx = iface + + with ProcessLogger(logger, "Building (possibly) unflipped " + "FacialAdjacencyGroups"): + unflipped_facial_adjacency_groups = \ + _get_firedrake_facial_adjacency_groups(fdrake_mesh, + cells_to_use=cells_to_use) + + # applied below to take elements and element_faces + # (or neighbors and neighbor_faces) and flip in any faces that need to + # be flipped. + def flip_local_face_indices(faces, elements): + faces = np.copy(faces) + neg_elements = np.full(elements.shape, False) + # To handle neighbor case, we only need to flip at elements + # who have a neighbor, i.e. where neighbors is not a negative + # bitmask of bdy tags + neg_elements[elements >= 0] = (orient[elements[elements >= 0]] < 0) + no_zero = np.logical_and(neg_elements, faces == no_zero_face_ndx) + no_one = np.logical_and(neg_elements, faces == no_one_face_ndx) + faces[no_zero], faces[no_one] = no_one_face_ndx, no_zero_face_ndx + return faces + + # Create new facial adjacency groups that have been flipped + with ProcessLogger(logger, "Flipping FacialAdjacencyGroups"): + facial_adjacency_groups = [] + for igroup, fagrps in enumerate(unflipped_facial_adjacency_groups): + facial_adjacency_groups.append({}) + for ineighbor_group, fagrp in fagrps.items(): + new_element_faces = flip_local_face_indices(fagrp.element_faces, + fagrp.elements) + new_neighbor_faces = flip_local_face_indices(fagrp.neighbor_faces, + fagrp.neighbors) + new_fagrp = FacialAdjacencyGroup(igroup=igroup, + ineighbor_group=ineighbor_group, + elements=fagrp.elements, + element_faces=new_element_faces, + neighbors=fagrp.neighbors, + neighbor_faces=new_neighbor_faces) + facial_adjacency_groups[igroup][ineighbor_group] = new_fagrp + + return (Mesh(vertices, [group], + boundary_tags=bdy_tags, + nodal_adjacency=nodal_adjacency, + facial_adjacency_groups=facial_adjacency_groups), + orient) + +# }}} + + +# {{{ Mesh exporting to firedrake + +def export_mesh_to_firedrake(mesh, group_nr=None, comm=None): + r""" + Create a firedrake mesh corresponding to one + :class:`~meshmode.mesh.Mesh`'s + :class:`~meshmode.mesh.SimplexElementGroup`. + + :param mesh: A :class:`~meshmode.mesh.Mesh` to convert with + at least one :class:`~meshmode.mesh.SimplexElementGroup`. + 'mesh.is_conforming' must evaluate to *True*. + 'mesh' must have vertices supplied, i.e. + 'mesh.vertices' must not be *None*. + :param group_nr: The group number to be converted into a firedrake + mesh. The corresponding group must be of type + :class:`~meshmode.mesh.SimplexElementGroup`. If *None* and + *mesh* has only one group, that group is used. Otherwise, + a *ValueError* is raised. + :param comm: The communicator to build the dmplex mesh on + + :return: A tuple *(fdrake_mesh, fdrake_cell_ordering, perm2cell)* + where + + * *fdrake_mesh* is a :mod:`firedrake` + :class:`~firedrake.mesh.MeshGeometry` corresponding to + *mesh* + * *fdrake_cell_ordering* is a numpy array whose *i*\ th + element in *mesh* (i.e. the *i*\ th element in + *mesh.groups[group_nr].vertex_indices*) corresponds to the + *fdrake_cell_ordering[i]*\ th :mod:`firedrake` cell + * *perm2cell* is a dictionary, mapping tuples to + 1-D numpy arrays of meshmode element indices. + Each meshmode element index + appears in exactly one of these arrays. The corresponding + tuple describes how firedrake reordered the local vertex + indices on that cell. In particular, if *c* + is in the list *perm2cell[p]* for a tuple *p*, then + the *p[i]*\ th local vertex of the *fdrake_cell_ordering[c]*\ th + firedrake cell corresponds to the *i*\ th local vertex + of the *c*\ th meshmode element. + + .. warning:: + Currently, no custom boundary tags are exported along with the mesh. + :mod:`firedrake` seems to only allow one marker on each facet, whereas + :mod:`meshmode` allows many. + """ + if not isinstance(mesh, Mesh): + raise TypeError("'mesh' must of type meshmode.mesh.Mesh," + " not '%s'." % type(mesh)) + if group_nr is None: + if len(mesh.groups) != 1: + raise ValueError("'group_nr' is *None* but 'mesh' has " + "more than one group.") + group_nr = 0 + if not isinstance(group_nr, int): + raise TypeError("Expecting 'group_nr' to be of type int, not " + f"'{type(group_nr)}'") + if group_nr < 0 or group_nr >= len(mesh.groups): + raise ValueError("'group_nr' is an invalid group index:" + f" '{group_nr}' fails to satisfy " + f"0 <= {group_nr} < {len(mesh.groups)}") + if not isinstance(mesh.groups[group_nr], SimplexElementGroup): + raise TypeError("Expecting 'mesh.groups[group_nr]' to be of type " + "meshmode.mesh.SimplexElementGroup, not " + f"'{type(mesh.groups[group_nr])}'") + if mesh.vertices is None: + raise ValueError("'mesh' has no vertices " + "('mesh.vertices' is *None*)") + if not mesh.is_conforming: + raise ValueError(f"'mesh.is_conforming' is {mesh.is_conforming} " + "instead of *True*. Converting non-conforming " + " meshes to Firedrake is not supported") + + # Get the vertices and vertex indices of the requested group + with ProcessLogger(logger, "Obtaining vertices from selected group"): + group = mesh.groups[group_nr] + fd2mm_indices = np.unique(group.vertex_indices.flatten()) + coords = mesh.vertices[:, fd2mm_indices].T + mm2fd_indices = dict(zip(fd2mm_indices, np.arange(np.size(fd2mm_indices)))) + cells = np.vectorize(mm2fd_indices.__getitem__)(group.vertex_indices) + + # Get a dmplex object and then a mesh topology + with ProcessLogger(logger, "Building dmplex object and MeshTopology"): + if comm is None: + from pyop2.mpi import COMM_WORLD + comm = COMM_WORLD + # FIXME : not sure how to get around the private accesses + import firedrake.mesh as fd_mesh + plex = fd_mesh._from_cell_list(group.dim, cells, coords, comm) + # Nb : One might be tempted to pass reorder=False and thereby save some + # hassle in exchange for forcing firedrake to have slightly + # less efficient caching. Unfortunately, that only prevents + # the cells from being reordered, and does not prevent the + # vertices from being (locally) reordered on each cell... + # the tl;dr is we don't actually save any hassle + top = fd_mesh.Mesh(plex, dim=mesh.ambient_dim) # mesh topology + top.init() + + # Get new element ordering: + with ProcessLogger(logger, "Determining permutations applied" + " to local vertex numbers"): + c_start, c_end = top._topology_dm.getHeightStratum(0) + cell_index_mm2fd = np.vectorize(top._cell_numbering.getOffset)( + np.arange(c_start, c_end)) + v_start, v_end = top._topology_dm.getDepthStratum(0) + + # Firedrake goes crazy reordering local vertex numbers, + # we've got to work to figure out what changes they made. + # + # *perm2cells* will map permutations of local vertex numbers to + # the list of all the meshmode cells + # which firedrake reordered according to that permutation + # + # Permutations on *n* vertices are stored as a tuple + # containing all of the integers *0*, *1*, *2*, ..., *n-1* + # exactly once. A permutation *p* + # represents relabeling the *i*\ th local vertex + # of a meshmode element as the *p[i]*\ th local vertex + # in the corresponding firedrake cell. + # + # *perm2cells[p]* is a list of all the meshmode element indices + # for which *p* represents the reordering applied by firedrake + perm2cells = {} + for mm_cell_id, dmp_ids in enumerate(top.cell_closure[cell_index_mm2fd]): + # look at order of vertices in firedrake cell + vert_dmp_ids = \ + dmp_ids[np.logical_and(v_start <= dmp_ids, dmp_ids < v_end)] + fdrake_order = vert_dmp_ids - v_start + # get original order + mm_order = mesh.groups[group_nr].vertex_indices[mm_cell_id] + # want permutation p so that mm_order[p] = fdrake_order + # To do so, look at permutations acting by composition. + # + # mm_order \circ argsort(mm_order) = + # fdrake_order \circ argsort(fdrake_order) + # so + # mm_order \circ argsort(mm_order) \circ inv(argsort(fdrake_order)) + # = fdrake_order + # + # argsort acts as an inverse, so the desired permutation is: + perm = tuple(np.argsort(mm_order)[np.argsort(np.argsort(fdrake_order))]) + perm2cells.setdefault(perm, []) + perm2cells[perm].append(mm_cell_id) + + # Make perm2cells map to numpy arrays instead of lists + perm2cells = {perm: np.array(cells) + for perm, cells in perm2cells.items()} + + # Now make a coordinates function + with ProcessLogger(logger, "Building firedrake function " + "space for mesh coordinates"): + from firedrake import VectorFunctionSpace, Function + coords_fspace = VectorFunctionSpace(top, 'CG', group.order, + dim=mesh.ambient_dim) + coords = Function(coords_fspace) + + # get firedrake unit nodes and map onto meshmode reference element + fd_ref_cell_to_mm = get_affine_reference_simplex_mapping(group.dim, True) + fd_unit_nodes = get_finat_element_unit_nodes(coords_fspace.finat_element) + fd_unit_nodes = fd_ref_cell_to_mm(fd_unit_nodes) + + basis = simplex_best_available_basis(group.dim, group.order) + resampling_mat = resampling_matrix(basis, + new_nodes=fd_unit_nodes, + old_nodes=group.unit_nodes) + # Store the meshmode data resampled to firedrake unit nodes + # (but still in meshmode order) + resampled_group_nodes = np.matmul(group.nodes, resampling_mat.T) + + # Now put the nodes in the right local order + # nodes is shaped *(ambient dim, nelements, nunit nodes)* + with ProcessLogger(logger, "Storing meshmode mesh coordinates" + " in firedrake nodal order"): + from meshmode.mesh.processing import get_simplex_element_flip_matrix + for perm, cells in perm2cells.items(): + flip_mat = get_simplex_element_flip_matrix(group.order, + fd_unit_nodes, + perm) + flip_mat = np.rint(flip_mat).astype(np.int32) + resampled_group_nodes[:, cells, :] = \ + np.matmul(resampled_group_nodes[:, cells, :], flip_mat.T) + + # store resampled data in right cell ordering + with ProcessLogger(logger, "resampling mesh coordinates to " + "firedrake unit nodes"): + reordered_cell_node_list = coords_fspace.cell_node_list[cell_index_mm2fd] + coords.dat.data[reordered_cell_node_list, :] = \ + resampled_group_nodes.transpose((1, 2, 0)) + + return fd_mesh.Mesh(coords), cell_index_mm2fd, perm2cells + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/interop/firedrake/reference_cell.py b/meshmode/interop/firedrake/reference_cell.py new file mode 100644 index 0000000000000000000000000000000000000000..327f9cd4236503dbdbcddafb8183e9d0241d0600 --- /dev/null +++ b/meshmode/interop/firedrake/reference_cell.py @@ -0,0 +1,166 @@ +__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 + + +__doc__ = """ +.. autofunction:: get_affine_reference_simplex_mapping +.. autofunction:: get_finat_element_unit_nodes +""" + + +# {{{ Map between reference simplices + +def get_affine_reference_simplex_mapping(ambient_dim, firedrake_to_meshmode=True): + """ + Returns a function which takes a numpy array points + on one reference cell and maps each + point to another using a positive affine map. + + :arg ambient_dim: The spatial dimension + :arg firedrake_to_meshmode: If true, the returned function maps from + the firedrake reference element to + meshmode, if false maps from + meshmode to firedrake. More specifically, + :mod:`firedrake` uses the standard :mod:`FIAT` + simplex and :mod:`meshmode` uses + :mod:`modepy`'s + `unit coordinates `_. + :return: A function which takes a numpy array of *n* points with + shape *(dim, n)* on one reference cell and maps + each point to another using a positive affine map. + Note that the returned function performs + no input validation. + """ + # validate input + if not isinstance(ambient_dim, int): + raise TypeError("'ambient_dim' must be an int, not " + f"'{type(ambient_dim)}'") + if ambient_dim < 0: + raise ValueError("'ambient_dim' must be non-negative") + if not isinstance(firedrake_to_meshmode, bool): + raise TypeError("'firedrake_to_meshmode' must be a bool, not " + f"'{type(firedrake_to_meshmode)}'") + + from FIAT.reference_element import ufc_simplex + from modepy.tools import unit_vertices + # Get the unit vertices from each system, + # each stored with shape *(dim, nunit_vertices)* + firedrake_unit_vertices = np.array(ufc_simplex(ambient_dim).vertices).T + modepy_unit_vertices = unit_vertices(ambient_dim).T + + if firedrake_to_meshmode: + from_verts = firedrake_unit_vertices + to_verts = modepy_unit_vertices + else: + from_verts = modepy_unit_vertices + to_verts = firedrake_unit_vertices + + # Compute matrix A and vector b so that A f_i + b -> t_i + # for each "from" vertex f_i and corresponding "to" vertex t_i + assert from_verts.shape == to_verts.shape + dim, nvects = from_verts.shape + + # If we only have one vertex, have A = I and b = to_vert - from_vert + if nvects == 1: + shift = to_verts[:, 0] - from_verts[:, 0] + + def affine_map(points): + return points + shift[:, np.newaxis] + # Otherwise, we have to solve for A and b + else: + # span verts: v1 - v0, v2 - v0, ... + from_span_verts = from_verts[:, 1:] - from_verts[:, 0, np.newaxis] + to_span_verts = to_verts[:, 1:] - to_verts[:, 0, np.newaxis] + # mat maps (fj - f0) -> (tj - t0), our "A" + mat = la.solve(from_span_verts, to_span_verts) + # A f0 + b -> t0 so b = t0 - A f0 + shift = to_verts[:, 0] - np.matmul(mat, from_verts[:, 0]) + + # Explicitly ensure A is positive + if la.det(mat) < 0: + from meshmode.mesh.processing import get_simplex_element_flip_matrix + flip_matrix = get_simplex_element_flip_matrix(1, to_verts) + mat = np.matmul(flip_matrix, mat) + + def affine_map(points): + return np.matmul(mat, points) + shift[:, np.newaxis] + + return affine_map + +# }}} + + +# {{{ Get firedrake unit nodes + +def get_finat_element_unit_nodes(finat_element): + """ + Returns the unit nodes used by the :mod:`finat` element in firedrake's + (equivalently, :mod:`finat`/:mod:`FIAT`'s) reference coordinates + + :arg finat_element: An instance of one of the following :mod:`finat` + elements + + * :class:`finat.fiat_elements.Lagrange` + * :class:`finat.fiat_elements.DiscontinuousLagrange` + * :class:`finat.fiat_elements.CrouzeixRaviart` + * :class:`finat.spectral.GaussLobattoLegendre` + * :class:`finat.spectral.GaussLegendre` + + :return: A numpy array of shape *(dim, nunit_dofs)* holding the unit + nodes used by this element. *dim* is the dimension spanned + by the finat element's reference element + (see its ``cell`` attribute) + """ + from finat.fiat_elements import ( + Lagrange, DiscontinuousLagrange, CrouzeixRaviart) + from finat.spectral import GaussLobattoLegendre, GaussLegendre + from FIAT.reference_element import Simplex + allowed_finat_elts = (Lagrange, DiscontinuousLagrange, CrouzeixRaviart, + GaussLobattoLegendre, GaussLegendre) + if not isinstance(finat_element, allowed_finat_elts): + raise TypeError("'finat_element' is of unexpected type " + f"{type(finat_element)}. 'finat_element' must be of " + "one of the following types: {allowed_finat_elts}") + if not isinstance(finat_element.cell, Simplex): + raise TypeError("Reference element of the finat element MUST be a" + " simplex, i.e. 'finat_element's *cell* attribute must" + " be of type FIAT.reference_element.Simplex, not " + f"'{type(finat_element.cell)}'") + # We insisted 'finat_element._element' was a + # FIAT.finite_element.CiarletElement, + # so the finat_element._element.dual.nodes ought to represent + # nodal dofs + # + # point evaluators is a list of functions *p_0,...,p_{n-1}*. + # *p_i(f)* evaluates function *f* at node *i* (stored as a tuple), + # so to recover node *i* we need to evaluate *p_i* at the identity + # function + point_evaluators = finat_element._element.dual.nodes + unit_nodes = [p(lambda x: x) for p in point_evaluators] + return np.array(unit_nodes).T + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index c5d7ecd8f247f8a8ca93144c38e3aae3ff9e07bd..532e5e392307be4fcefeb3566297b8feb2ab62de 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -53,33 +53,46 @@ Predefined Boundary tags .. autoclass:: BTAG_REALLY_ALL .. autoclass:: BTAG_NO_BOUNDARY .. autoclass:: BTAG_PARTITION +.. autoclass:: BTAG_INDUCED_BOUNDARY """ # {{{ element tags -class BTAG_NONE(object): # noqa +class BTAG_NONE(object): # noqa: N801 """A boundary tag representing an empty boundary or volume.""" pass -class BTAG_ALL(object): # noqa +class BTAG_ALL(object): # noqa: N801 """A boundary tag representing the entire boundary or volume. In the case of the boundary, :class:`BTAG_ALL` does not include rank boundaries, - or, more generally, anything tagged with :class:`BTAG_NO_BOUNDARY`.""" + or, more generally, anything tagged with :class:`BTAG_NO_BOUNDARY`. + + In the case of a mesh representing an element-wise subset of another, + :class:`BTAG_ALL` does not include boundaries induced by taking the subset. + Instead, these boundaries will be tagged with + :class:`BTAG_INDUCED_BOUNDARY`. + """ pass -class BTAG_REALLY_ALL(object): # noqa +class BTAG_REALLY_ALL(object): # noqa: N801 """A boundary tag representing the entire boundary. Unlike :class:`BTAG_ALL`, this includes rank boundaries, - or, more generally, everything tagged with :class:`BTAG_NO_BOUNDARY`.""" + or, more generally, everything tagged with :class:`BTAG_NO_BOUNDARY`. + + In the case of a mesh representing an element-wise subset of another, + this tag includes boundaries induced by taking the subset, or, more generally, + everything tagged with + :class:`BTAG_INDUCED_BOUNDARY` + """ pass -class BTAG_NO_BOUNDARY(object): # noqa +class BTAG_NO_BOUNDARY(object): # noqa: N801 """A boundary tag indicating that this edge should not fall under :class:`BTAG_ALL`. Among other things, this is used to keep rank boundaries out of :class:`BTAG_ALL`. @@ -87,7 +100,7 @@ class BTAG_NO_BOUNDARY(object): # noqa pass -class BTAG_PARTITION(object): # noqa +class BTAG_PARTITION(object): # noqa: N801 """ A boundary tag indicating that this edge is adjacent to an element of another :class:`Mesh`. The partition number of the adjacent mesh @@ -116,8 +129,19 @@ class BTAG_PARTITION(object): # noqa return "<%s(%s)>" % (type(self).__name__, repr(self.part_nr)) +class BTAG_INDUCED_BOUNDARY(BTAG_NO_BOUNDARY): # noqa: N801 + """When a :class:`Mesh` is created as an element-by-element subset of another + (as, for example, when calling + :func:`meshmode.interop.firedrake.build_connection_from_firedrake` + while passing *restrict_to_boundary*), boundaries may arise where there + were none in the original mesh. This boundary tag is used to indicate + such boundaries. + """ + pass + + SYSTEM_TAGS = set([BTAG_NONE, BTAG_ALL, BTAG_REALLY_ALL, BTAG_NO_BOUNDARY, - BTAG_PARTITION]) + BTAG_PARTITION, BTAG_INDUCED_BOUNDARY]) # }}} diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 7e82ea00852ae93b9e175f5ef446da27e11df936..6ca668de9d41659a92fd0e193e6437efc3926469 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -723,9 +723,58 @@ def test_volume_mesh_element_orientations(mesh): # {{{ flips -def flip_simplex_element_group(vertices, grp, grp_flip_flags): + +def get_simplex_element_flip_matrix(order, unit_nodes, permutation=None): + """ + Generate a resampling matrix that corresponds to a + permutation of the barycentric coordinates being applied. + The default permutation is to swap the + first two barycentric coordinates. + + :param order: The order of the function space on the simplex, + (see second argument in + :fun:`modepy.simplex_best_available_basis`) + :param unit_nodes: A np array of unit nodes with shape + *(dim, nunit_nodes)* + :param permutation: Either *None*, or a tuple of shape + storing a permutation: + the *i*th barycentric coordinate gets mapped to + the *permutation[i]*th coordinate. + + :return: A numpy array of shape *(nunit_nodes, nunit_nodes)* + which, when its transpose is right-applied + to the matrix of nodes (shaped *(dim, nunit_nodes)*), + corresponds to the permutation being applied + """ from modepy.tools import barycentric_to_unit, unit_to_barycentric + bary_unit_nodes = unit_to_barycentric(unit_nodes) + + flipped_bary_unit_nodes = bary_unit_nodes.copy() + if permutation is None: + flipped_bary_unit_nodes[0, :] = bary_unit_nodes[1, :] + flipped_bary_unit_nodes[1, :] = bary_unit_nodes[0, :] + else: + flipped_bary_unit_nodes[permutation, :] = bary_unit_nodes + flipped_unit_nodes = barycentric_to_unit(flipped_bary_unit_nodes) + + dim = unit_nodes.shape[0] + flip_matrix = mp.resampling_matrix( + mp.simplex_best_available_basis(dim, order), + flipped_unit_nodes, unit_nodes) + + flip_matrix[np.abs(flip_matrix) < 1e-15] = 0 + + # Flipping twice should be the identity + if permutation is None: + assert la.norm( + np.dot(flip_matrix, flip_matrix) + - np.eye(len(flip_matrix))) < 1e-13 + + return flip_matrix + + +def flip_simplex_element_group(vertices, grp, grp_flip_flags): from meshmode.mesh import SimplexElementGroup if not isinstance(grp, SimplexElementGroup): @@ -740,28 +789,8 @@ def flip_simplex_element_group(vertices, grp, grp_flip_flags): new_vertex_indices[grp_flip_flags, 1] \ = grp.vertex_indices[grp_flip_flags, 0] - # Generate a resampling matrix that corresponds to the - # first two barycentric coordinates being swapped. - - bary_unit_nodes = unit_to_barycentric(grp.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) - - flip_matrix = mp.resampling_matrix( - mp.simplex_best_available_basis(grp.dim, grp.order), - flipped_unit_nodes, grp.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 - # Apply the flip matrix to the nodes. + flip_matrix = get_simplex_element_flip_matrix(grp.order, grp.unit_nodes) new_nodes = grp.nodes.copy() new_nodes[:, grp_flip_flags] = np.einsum( "ij,dej->dei", diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py new file mode 100644 index 0000000000000000000000000000000000000000..17cabade91abba74c831f7fb4d09767400c035c9 --- /dev/null +++ b/test/test_firedrake_interop.py @@ -0,0 +1,688 @@ +__copyright__ = "Copyright (C) 2020 Benjamin Sepanski" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np +import pyopencl as cl +import six + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl + as pytest_generate_tests) + +from meshmode.array_context import PyOpenCLArrayContext + +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) + +from meshmode.dof_array import DOFArray + +from meshmode.mesh import ( + BTAG_ALL, BTAG_REALLY_ALL, BTAG_INDUCED_BOUNDARY, check_bc_coverage + ) + +from meshmode.interop.firedrake import ( + build_connection_from_firedrake, build_connection_to_firedrake, + import_firedrake_mesh) + +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, TensorFunctionSpace, + Function, SpatialCoordinate, as_tensor, sin) + + +CLOSE_ATOL = 1e-12 + + +@pytest.fixture(params=["annulus.msh", + "blob2d-order1-h4e-2.msh", + "blob2d-order1-h6e-2.msh", + "blob2d-order1-h8e-2.msh", + "blob2d-order4-h4e-2.msh", + "blob2d-order4-h6e-2.msh", + "blob2d-order4-h8e-2.msh", + ]) +def mm_mesh(request): + from meshmode.mesh.io import read_gmsh + return read_gmsh(request.param) + + +@pytest.fixture(params=["FiredrakeUnitIntervalMesh", + "FiredrakeUnitSquareMesh", + "FiredrakeUnitSquareMesh-order2", + "FiredrakeUnitCubeMesh", + "annulus.msh", + "blob2d-order1-h4e-2.msh", + "blob2d-order1-h6e-2.msh", + "blob2d-order1-h8e-2.msh", + ]) +def fdrake_mesh(request): + mesh_name = request.param + if mesh_name == "FiredrakeUnitIntervalMesh": + return UnitIntervalMesh(100) + elif mesh_name == "FiredrakeUnitSquareMesh": + return UnitSquareMesh(10, 10) + elif mesh_name == "FiredrakeUnitSquareMesh-order2": + m = UnitSquareMesh(10, 10) + fspace = VectorFunctionSpace(m, 'CG', 2) + coords = Function(fspace).interpolate(SpatialCoordinate(m)) + from firedrake.mesh import Mesh + return Mesh(coords) + elif mesh_name == "FiredrakeUnitCubeMesh": + return UnitCubeMesh(5, 5, 5) + elif mesh_name not in ("annulus.msh", "blob2d-order1-h4e-2.msh", + "blob2d-order1-h6e-2.msh", "blob2d-order1-h8e-2.msh"): + raise ValueError("Unexpected value for request.param") + + # Firedrake can't read in higher order meshes from gmsh, + # so we can only use the order1 blobs + from firedrake import Mesh + fd_mesh = Mesh(mesh_name) + fd_mesh.init() + return fd_mesh + + +@pytest.fixture(params=[1, 4], ids=["P^1", "P^4"]) +def fspace_degree(request): + return request.param + + +# {{{ Basic conversion checks for the function space + +def check_consistency(fdrake_fspace, discr, group_nr=0): + """ + While nodes may change, vertex conversion should be *identical* up to + reordering, ensure this is the case for DG spaces. Also ensure the + meshes have the same basic properties and the function space/discretization + agree across firedrake vs meshmode + """ + # Get the unit vertex indices (in each cell) + fdrake_mesh = fdrake_fspace.mesh() + cfspace = fdrake_mesh.coordinates.function_space() + entity_dofs = cfspace.finat_element.entity_dofs()[0] + fdrake_unit_vert_indices = [] + for _, local_node_nrs in sorted(six.iteritems(entity_dofs)): + assert len(local_node_nrs) == 1 + fdrake_unit_vert_indices.append(local_node_nrs[0]) + + # get the firedrake vertices, in no particular order + fdrake_vert_indices = cfspace.cell_node_list[:, fdrake_unit_vert_indices] + fdrake_vert_indices = np.unique(fdrake_vert_indices) + fdrake_verts = fdrake_mesh.coordinates.dat.data[fdrake_vert_indices, ...] + if fdrake_mesh.geometric_dimension() == 1: + fdrake_verts = fdrake_verts[:, np.newaxis] + + meshmode_verts = discr.mesh.vertices + + # Ensure the meshmode mesh has one group and make sure both + # meshes agree on some basic properties + assert len(discr.mesh.groups) == 1 + fdrake_mesh_fspace = fdrake_mesh.coordinates.function_space() + fdrake_mesh_order = fdrake_mesh_fspace.finat_element.degree + assert discr.mesh.groups[group_nr].dim == fdrake_mesh.topological_dimension() + assert discr.mesh.groups[group_nr].order == fdrake_mesh_order + assert discr.mesh.groups[group_nr].nelements == fdrake_mesh.num_cells() + assert 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_allclose(lex_sorted_mm_verts, lex_sorted_fdrake_verts.T, + atol=1e-15) + + # Ensure the discretization and the firedrake function space agree on + # some basic properties + finat_elt = fdrake_fspace.finat_element + assert len(discr.groups) == 1 + assert discr.groups[group_nr].order == finat_elt.degree + assert discr.groups[group_nr].nunit_dofs == finat_elt.space_dimension() + assert discr.ndofs == fdrake_fspace.node_count + + +def test_from_fd_consistency(ctx_factory, fdrake_mesh, fspace_degree): + """ + Check basic consistency with a FiredrakeConnection built from firedrake + """ + # make discretization from firedrake + fdrake_fspace = FunctionSpace(fdrake_mesh, 'DG', fspace_degree) + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + fdrake_connection = build_connection_from_firedrake(actx, fdrake_fspace) + discr = fdrake_connection.discr + # Check consistency + check_consistency(fdrake_fspace, discr) + + +def test_to_fd_consistency(ctx_factory, mm_mesh, fspace_degree): + fspace_degree += mm_mesh.groups[0].order + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + factory = InterpolatoryQuadratureSimplexGroupFactory(fspace_degree) + discr = Discretization(actx, mm_mesh, factory) + fdrake_connection = build_connection_to_firedrake(discr) + fdrake_fspace = fdrake_connection.firedrake_fspace() + # Check consistency + check_consistency(fdrake_fspace, discr) + +# }}} + + +# {{{ Now check the FiredrakeConnection consistency when restricted to bdy + +def test_from_boundary_consistency(ctx_factory, + fdrake_mesh, + fspace_degree): + """ + Make basic checks that FiredrakeConnection restricted to cells + near the boundary is not doing + something obviously wrong, + i.e. that the firedrake boundary tags partition the converted meshmode mesh, + that the firedrake boundary tags correspond to the same physical + regions in the converted meshmode mesh as in the original firedrake mesh, + and that each boundary tag is associated to the same number of facets + in the converted meshmode mesh as in the original firedrake mesh. + """ + fdrake_fspace = FunctionSpace(fdrake_mesh, 'DG', fspace_degree) + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + frombdy_conn = \ + build_connection_from_firedrake(actx, + fdrake_fspace, + restrict_to_boundary="on_boundary") + + # Ensure the meshmode mesh has one group and make sure both + # meshes agree on some basic properties + discr = frombdy_conn.discr + assert len(discr.mesh.groups) == 1 + fdrake_mesh_fspace = fdrake_mesh.coordinates.function_space() + fdrake_mesh_order = fdrake_mesh_fspace.finat_element.degree + assert discr.mesh.groups[0].dim == fdrake_mesh.topological_dimension() + assert discr.mesh.groups[0].order == fdrake_mesh_order + + # Get the unit vertex indices (in each cell) + fdrake_mesh = fdrake_fspace.mesh() + cfspace = fdrake_mesh.coordinates.function_space() + entity_dofs = cfspace.finat_element.entity_dofs()[0] + fdrake_unit_vert_indices = [] + for _, local_node_nrs in sorted(six.iteritems(entity_dofs)): + assert len(local_node_nrs) == 1 + fdrake_unit_vert_indices.append(local_node_nrs[0]) + fdrake_unit_vert_indices = np.array(fdrake_unit_vert_indices) + + # only look at cells "near" bdy (with >= 1 vertex on) + from meshmode.interop.firedrake.connection import _get_cells_to_use + cells_near_bdy = _get_cells_to_use(fdrake_mesh, 'on_boundary') + # get the firedrake vertices of cells near the boundary, + # in no particular order + fdrake_vert_indices = \ + cfspace.cell_node_list[cells_near_bdy, + fdrake_unit_vert_indices[:, np.newaxis]] + fdrake_vert_indices = np.unique(fdrake_vert_indices) + fdrake_verts = fdrake_mesh.coordinates.dat.data[fdrake_vert_indices, ...] + if fdrake_mesh.geometric_dimension() == 1: + fdrake_verts = fdrake_verts[:, np.newaxis] + # Get meshmode vertices (shaped like (dim, nverts)) + meshmode_verts = discr.mesh.vertices + + # Ensure that the vertices of firedrake elements on + # the boundary are identical to the resultant meshes' vertices 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_allclose(lex_sorted_mm_verts, lex_sorted_fdrake_verts.T, + atol=CLOSE_ATOL) + + # Ensure the discretization and the firedrake function space reference element + # agree on some basic properties + finat_elt = fdrake_fspace.finat_element + assert len(discr.groups) == 1 + assert discr.groups[0].order == finat_elt.degree + assert discr.groups[0].nunit_dofs == finat_elt.space_dimension() + +# }}} + + +# {{{ Boundary tags checking + +bdy_tests = [(UnitSquareMesh(10, 10), + [1, 2, 3, 4], + [0, 0, 1, 1], + [0.0, 1.0, 0.0, 1.0]), + (UnitCubeMesh(5, 5, 5), + [1, 2, 3, 4, 5, 6], + [0, 0, 1, 1, 2, 2], + [0.0, 1.0, 0.0, 1.0, 0.0, 1.0]), + ] + + +@pytest.mark.parametrize("square_or_cube_mesh,bdy_ids,coord_indices,coord_values", + bdy_tests) +@pytest.mark.parametrize("only_convert_bdy", (True, False)) +def test_bdy_tags(square_or_cube_mesh, bdy_ids, coord_indices, coord_values, + only_convert_bdy): + """ + Make sure the given boundary ids cover the converted mesh. + Make sure that the given coordinate have the given value for the + corresponding boundary tag (see :mod:`firedrake.utility_meshes`'s + documentation to see how the boundary tags for its utility meshes are + defined) + """ + cells_to_use = None + if only_convert_bdy: + from meshmode.interop.firedrake.connection import _get_cells_to_use + cells_to_use = _get_cells_to_use(square_or_cube_mesh, 'on_boundary') + mm_mesh, orient = import_firedrake_mesh(square_or_cube_mesh, + cells_to_use=cells_to_use) + # Ensure meshmode required boundary tags are there + assert set([BTAG_ALL, BTAG_REALLY_ALL]) <= set(mm_mesh.boundary_tags) + # Check disjoint coverage of bdy ids and BTAG_ALL + check_bc_coverage(mm_mesh, [BTAG_ALL]) + check_bc_coverage(mm_mesh, bdy_ids) + + # count number of times the boundary tag appears in the meshmode mesh, + # should be the same as in the firedrake mesh + bdy_id_to_mm_count = {} + # Now make sure we have identified the correct faces + face_vertex_indices = mm_mesh.groups[0].face_vertex_indices() + ext_grp = mm_mesh.facial_adjacency_groups[0][None] + for iel, ifac, bdy_flags in zip( + ext_grp.elements, ext_grp.element_faces, ext_grp.neighbors): + # try: if mm_mesh has boundaries flagged as not boundaries we need to + # skip them + # catch: if mm_mesh does not use BTAG_INDUCED_BOUNDARY flag we get a + # ValueError + try: + # If this facet is flagged as not really a boundary, skip it + if mm_mesh.boundary_tag_bit(BTAG_INDUCED_BOUNDARY) & -bdy_flags: + continue + except ValueError: + pass + + el_vert_indices = mm_mesh.groups[0].vertex_indices[iel] + # numpy nb: have to have comma to use advanced indexing + face_vert_indices = el_vert_indices[face_vertex_indices[ifac], ] + # shape: *(ambient dim, num vertices on face)* + face_verts = mm_mesh.vertices[:, face_vert_indices] + # Figure out which coordinate should have a fixed value, and what + # that value is. Also, count how many times each boundary tag appears + coord_index, val = None, None + for bdy_id_index, bdy_id in enumerate(bdy_ids): + if mm_mesh.boundary_tag_bit(bdy_id) & -bdy_flags: + bdy_id_to_mm_count.setdefault(bdy_id, 0) + bdy_id_to_mm_count[bdy_id] += 1 + coord_index = coord_indices[bdy_id_index] + val = coord_values[bdy_id_index] + break + assert np.max(np.abs(face_verts[coord_index, :] - val)) < CLOSE_ATOL + + # Verify that the number of meshes tagged with a boundary tag + # is the same in meshmode and firedrake for each tag in *bdy_ids* + fdrake_bdy_ids, fdrake_counts = \ + np.unique(square_or_cube_mesh.exterior_facets.markers, return_counts=True) + assert set(fdrake_bdy_ids) == set(bdy_ids) + for bdy_id, fdrake_count in zip(fdrake_bdy_ids, fdrake_counts): + assert fdrake_count == bdy_id_to_mm_count[bdy_id] + +# }}} + + +# TODO : Add test for FiredrakeConnection built from meshmode +# where group_nr != 0 + +# {{{ Double check functions are being transported correctly + + +@pytest.mark.parametrize("fdrake_mesh_name,fdrake_mesh_pars,dim", + [("UnitInterval", [10, 20, 30], 1), + ("UnitSquare", [10, 20, 30], 2), + ("UnitCube", [10, 20, 30], 3), + ("blob2d-order1", ["8e-2", "6e-2", "4e-2"], 2), + ("blob2d-order4", ["8e-2", "6e-2", "4e-2"], 2), + ("warp", [10, 20, 30], 2), + ("warp", [10, 20, 30], 3), + ]) +@pytest.mark.parametrize("only_convert_bdy", [False, True]) +def test_from_fd_transfer(ctx_factory, fspace_degree, + fdrake_mesh_name, fdrake_mesh_pars, dim, + only_convert_bdy): + """ + Make sure creating a function which projects onto + one dimension then transports it is the same + (up to resampling error) as projecting to one + dimension on the transported mesh + """ + # build estimate-of-convergence recorder + from pytools.convergence import EOCRecorder + # (fd -> mm ? True : False, dimension projecting onto) + eoc_recorders = {(True, d): EOCRecorder() for d in range(dim)} + if not only_convert_bdy: + for d in range(dim): + eoc_recorders[(False, d)] = EOCRecorder() + + # make a computing context + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + def get_fdrake_mesh_and_h_from_par(mesh_par): + if fdrake_mesh_name == "UnitInterval": + assert dim == 1 + n = mesh_par + fdrake_mesh = UnitIntervalMesh(n) + h = 1/n + elif fdrake_mesh_name == "UnitSquare": + assert dim == 2 + n = mesh_par + fdrake_mesh = UnitSquareMesh(n, n) + h = 1/n + elif fdrake_mesh_name == "UnitCube": + assert dim == 3 + n = mesh_par + fdrake_mesh = UnitCubeMesh(n, n, n) + h = 1/n + elif fdrake_mesh_name in ("blob2d-order1", "blob2d-order4"): + assert dim == 2 + if fdrake_mesh_name == "blob2d-order1": + from firedrake import Mesh + fdrake_mesh = Mesh("%s-h%s.msh" % (fdrake_mesh_name, mesh_par), + dim=dim) + else: + from meshmode.mesh.io import read_gmsh + from meshmode.interop.firedrake import export_mesh_to_firedrake + mm_mesh = read_gmsh("%s-h%s.msh" % (fdrake_mesh_name, mesh_par), + force_ambient_dim=dim) + fdrake_mesh, _, _ = export_mesh_to_firedrake(mm_mesh) + h = float(mesh_par) + elif fdrake_mesh_name == "warp": + from meshmode.mesh.generation import generate_warped_rect_mesh + from meshmode.interop.firedrake import export_mesh_to_firedrake + mm_mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par) + fdrake_mesh, _, _ = export_mesh_to_firedrake(mm_mesh) + h = 1/mesh_par + else: + raise ValueError("fdrake_mesh_name not recognized") + + return (fdrake_mesh, h) + + # Record error for each refinement of each mesh + for mesh_par in fdrake_mesh_pars: + fdrake_mesh, h = get_fdrake_mesh_and_h_from_par(mesh_par) + # make function space and build connection + fdrake_fspace = FunctionSpace(fdrake_mesh, 'DG', fspace_degree) + if only_convert_bdy: + fdrake_connection = \ + build_connection_from_firedrake(actx, + fdrake_fspace, + restrict_to_boundary='on_boundary') + else: + fdrake_connection = build_connection_from_firedrake(actx, fdrake_fspace) + # get this for making functions in firedrake + spatial_coord = SpatialCoordinate(fdrake_mesh) + + # get nodes in handier format for making meshmode functions + discr = fdrake_connection.discr + # nodes is np array (ambient_dim,) of DOFArray (ngroups,) + # of arrays (nelements, nunit_dofs), we want a single np array + # of shape (ambient_dim, nelements, nunit_dofs) + nodes = discr.nodes() + group_nodes = np.array([actx.to_numpy(dof_arr[0]) for dof_arr in nodes]) + + # Now, for each coordinate d, test transferring the function + # x -> sin(dth component of x) + for d in range(dim): + fdrake_f = Function(fdrake_fspace).interpolate(sin(spatial_coord[d])) + # transport fdrake function and put in numpy + fd2mm_f = fdrake_connection.from_firedrake(fdrake_f, actx=actx) + fd2mm_f = actx.to_numpy(fd2mm_f[0]) + meshmode_f = np.sin(group_nodes[d, :, :]) + + # record fd -> mm error + err = np.max(np.abs(fd2mm_f - meshmode_f)) + eoc_recorders[(True, d)].add_data_point(h, err) + + if not only_convert_bdy: + # now transport mm -> fd + meshmode_f_dofarr = discr.zeros(actx) + meshmode_f_dofarr[0][:] = meshmode_f + mm2fd_f = fdrake_connection.from_meshmode(meshmode_f_dofarr) + # record mm -> fd error + err = np.max(np.abs(fdrake_f.dat.data - mm2fd_f.dat.data)) + eoc_recorders[(False, d)].add_data_point(h, err) + + # assert that order is correct or error is "low enough" + for ((fd2mm, d), eoc_rec) in six.iteritems(eoc_recorders): + print("\nfiredrake -> meshmode: %s\nvector *x* -> *sin(x[%s])*\n" + % (fd2mm, d), eoc_rec) + assert ( + eoc_rec.order_estimate() >= fspace_degree + or eoc_rec.max_error() < 2e-14) + + +@pytest.mark.parametrize("mesh_name,mesh_pars,dim", + [("blob2d-order1", ["8e-2", "6e-2", "4e-2"], 2), + ("blob2d-order4", ["8e-2", "6e-2", "4e-2"], 2), + ("warp", [10, 20, 30], 2), + ("warp", [10, 20, 30], 3), + ]) +def test_to_fd_transfer(ctx_factory, fspace_degree, mesh_name, mesh_pars, dim): + """ + Make sure creating a function which projects onto + one dimension then transports it is the same + (up to resampling error) as projecting to one + dimension on the transported mesh + """ + # build estimate-of-convergence recorder + from pytools.convergence import EOCRecorder + # dimension projecting onto -> EOCRecorder + eoc_recorders = {d: EOCRecorder() for d in range(dim)} + + # make a computing context + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + # Get each of the refinements of the meshmeshes and record + # conversions errors + for mesh_par in mesh_pars: + if mesh_name in ("blob2d-order1", "blob2d-order4"): + assert dim == 2 + from meshmode.mesh.io import read_gmsh + mm_mesh = read_gmsh("%s-h%s.msh" % (mesh_name, mesh_par), + force_ambient_dim=dim) + h = float(mesh_par) + elif mesh_name == "warp": + from meshmode.mesh.generation import generate_warped_rect_mesh + mm_mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par) + h = 1/mesh_par + else: + raise ValueError("mesh_name not recognized") + + # Make discr and connect it to firedrake + factory = InterpolatoryQuadratureSimplexGroupFactory(fspace_degree) + discr = Discretization(actx, mm_mesh, factory) + + fdrake_connection = build_connection_to_firedrake(discr) + fdrake_fspace = fdrake_connection.firedrake_fspace() + spatial_coord = SpatialCoordinate(fdrake_fspace.mesh()) + + # get the group's nodes in a numpy array + nodes = discr.nodes() + group_nodes = np.array([actx.to_numpy(dof_arr[0]) for dof_arr in nodes]) + + for d in range(dim): + meshmode_f = discr.zeros(actx) + meshmode_f[0][:] = group_nodes[d, :, :] + + # connect to firedrake and evaluate expr in firedrake + fdrake_f = Function(fdrake_fspace).interpolate(spatial_coord[d]) + + # transport to firedrake and record error + mm2fd_f = fdrake_connection.from_meshmode(meshmode_f) + + err = np.max(np.abs(fdrake_f.dat.data - mm2fd_f.dat.data)) + eoc_recorders[d].add_data_point(h, err) + + # assert that order is correct or error is "low enough" + for d, eoc_rec in six.iteritems(eoc_recorders): + print("\nvector *x* -> *x[%s]*\n" % d, eoc_rec) + assert ( + eoc_rec.order_estimate() >= fspace_degree + or eoc_rec.max_error() < 2e-14) + +# }}} + + +# {{{ Idempotency tests fd->mm->fd and (fd->)mm->fd->mm for connection + +@pytest.mark.parametrize("fspace_type", ("scalar", "vector", "tensor")) +@pytest.mark.parametrize("only_convert_bdy", (False, True)) +def test_from_fd_idempotency(ctx_factory, + fdrake_mesh, fspace_degree, + fspace_type, only_convert_bdy): + """ + Make sure fd->mm->fd and (fd->)->mm->fd->mm are identity + """ + # Make a function space and a function with unique values at each node + if fspace_type == "scalar": + fdrake_fspace = FunctionSpace(fdrake_mesh, 'DG', fspace_degree) + # Just use the node nr + fdrake_unique = Function(fdrake_fspace) + fdrake_unique.dat.data[:] = np.arange(fdrake_unique.dat.data.shape[0]) + elif fspace_type == "vector": + fdrake_fspace = VectorFunctionSpace(fdrake_mesh, 'DG', fspace_degree) + # use the coordinates + xx = SpatialCoordinate(fdrake_fspace.mesh()) + fdrake_unique = Function(fdrake_fspace).interpolate(xx) + elif fspace_type == "tensor": + fdrake_fspace = TensorFunctionSpace(fdrake_mesh, 'DG', fspace_degree) + # use the coordinates, duplicated into the right tensor shape + xx = SpatialCoordinate(fdrake_fspace.mesh()) + dim = fdrake_fspace.mesh().geometric_dimension() + unique_expr = as_tensor([xx for _ in range(dim)]) + fdrake_unique = Function(fdrake_fspace).interpolate(unique_expr) + + # Make connection + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + # If only converting boundary, first go ahead and do one round of + # fd->mm->fd. This will zero out any degrees of freedom absent in + # the meshmode mesh (because they are not associated to cells + # with >= 1 node on the boundary) + # + # Otherwise, just continue as normal + if only_convert_bdy: + fdrake_connection = \ + build_connection_from_firedrake(actx, + fdrake_fspace, + restrict_to_boundary='on_boundary') + temp = fdrake_connection.from_firedrake(fdrake_unique, actx=actx) + fdrake_unique = fdrake_connection.from_meshmode(temp) + else: + fdrake_connection = build_connection_from_firedrake(actx, fdrake_fspace) + + # Test for idempotency fd->mm->fd + mm_field = fdrake_connection.from_firedrake(fdrake_unique, actx=actx) + fdrake_unique_copy = Function(fdrake_fspace) + fdrake_connection.from_meshmode(mm_field, out=fdrake_unique_copy) + + np.testing.assert_allclose(fdrake_unique_copy.dat.data, + fdrake_unique.dat.data, + atol=CLOSE_ATOL) + + # Test for idempotency (fd->)mm->fd->mm + mm_field_copy = fdrake_connection.from_firedrake(fdrake_unique_copy, + actx=actx) + if fspace_type == "scalar": + np.testing.assert_allclose(actx.to_numpy(mm_field_copy[0]), + actx.to_numpy(mm_field[0]), + atol=CLOSE_ATOL) + else: + for dof_arr_cp, dof_arr in zip(mm_field_copy.flatten(), + mm_field.flatten()): + np.testing.assert_allclose(actx.to_numpy(dof_arr_cp[0]), + actx.to_numpy(dof_arr[0]), + atol=CLOSE_ATOL) + + +def test_to_fd_idempotency(ctx_factory, mm_mesh, fspace_degree): + """ + Make sure mm->fd->mm and (mm->)->fd->mm->fd are identity + """ + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + # make sure degree is higher order than mesh + fspace_degree += mm_mesh.groups[0].order + + # Make a function space and a function with unique values at each node + factory = InterpolatoryQuadratureSimplexGroupFactory(fspace_degree) + discr = Discretization(actx, mm_mesh, factory) + fdrake_connection = build_connection_to_firedrake(discr) + fdrake_mesh = fdrake_connection.firedrake_fspace().mesh() + dtype = fdrake_mesh.coordinates.dat.data.dtype + + mm_unique = discr.zeros(actx, dtype=dtype) + unique_vals = np.arange(np.size(mm_unique[0]), dtype=dtype) + mm_unique[0].set(unique_vals.reshape(mm_unique[0].shape)) + mm_unique_copy = DOFArray.from_list(actx, [mm_unique[0].copy()]) + + # Test for idempotency mm->fd->mm + fdrake_unique = fdrake_connection.from_meshmode(mm_unique) + fdrake_connection.from_firedrake(fdrake_unique, out=mm_unique_copy) + + np.testing.assert_allclose(actx.to_numpy(mm_unique_copy[0]), + actx.to_numpy(mm_unique[0]), + atol=CLOSE_ATOL) + + # Test for idempotency (mm->)fd->mm->fd + fdrake_unique_copy = fdrake_connection.from_meshmode(mm_unique_copy) + np.testing.assert_allclose(fdrake_unique_copy.dat.data, + fdrake_unique.dat.data, + atol=CLOSE_ATOL) + +# }}} + +# vim: foldmethod=marker