diff --git a/meshmode/interop/firedrake/connection.py b/meshmode/interop/firedrake/connection.py
index b9cefaec80fbd83b3992c5fac7d67912a3164d96..619773305dfc3e39e3e4071f2ac72f66abf90b7d 100644
--- a/meshmode/interop/firedrake/connection.py
+++ b/meshmode/interop/firedrake/connection.py
@@ -26,6 +26,9 @@ __doc__ = """
import numpy as np
import numpy.linalg as la
+import six
+
+from modepy import resampling_matrix
from meshmode.interop.firedrake.mesh import import_firedrake_mesh
from meshmode.interop.firedrake.reference_cell import (
@@ -37,8 +40,6 @@ from meshmode.discretization.poly_element import \
InterpolatoryQuadratureSimplexGroupFactory
from meshmode.discretization import Discretization
-from modepy import resampling_matrix
-
def _reorder_nodes(orient, nodes, flip_matrix, unflip=False):
"""
@@ -130,19 +131,45 @@ class FromFiredrakeConnection:
new_nodes=fd_unit_nodes,
old_nodes=mm_unit_nodes)
- # handle reordering fd->mm
+ # 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.
+ #
+ # One should note there is something a bit more subtle going on
+ # in the continuous case. All meshmode discretizations use
+ # are discontinuous, so nodes are associated with elements(cells)
+ # not vertices. In a continuous firedrake space, some nodes
+ # are shared between multiple cells. In particular, while the
+ # below "reordering" is indeed a permutation if the firedrake space
+ # is discontinuous, if the firedrake space is continuous then
+ # some firedrake nodes correspond to nodes on multiple meshmode
+ # elements, i.e. those nodes appear multiple times
+ # in the "reordering" array
flip_mat = get_simplex_element_flip_matrix(ufl_elt.degree(),
fd_unit_nodes)
fd_cell_node_list = fdrake_fspace.cell_node_list
_reorder_nodes(orient, fd_cell_node_list, flip_mat, unflip=False)
self._reordering_arr_fd2mm = fd_cell_node_list.flatten()
- # handle reordering mm->fd (this only works in the discontinuous
- # case)
- nnodes = self.to_discr.nnodes
- mm_cell_node_list = self.to_discr.groups[0].view(np.arange(nnodes))
- _reorder_nodes(orient, mm_cell_node_list, flip_mat, unflip=True)
- self._reordering_arr_mm2fd = mm_cell_node_list.flatten()
+ # Now handle the possibility of duplicate nodes
+ unique_fd_nodes, counts = np.unique(self._reordering_arr_fd2mm,
+ return_counts=True)
+ # self._duplicate_nodes
+ # maps firedrake nodes associated to more than 1 meshmode node
+ # to all associated meshmode nodes.
+ self._duplicate_nodes = {}
+ if ufl_elt.family() == 'Discontinuous Lagrange':
+ assert np.all(counts == 1), \
+ "This error should never happen, some nodes in a firedrake " \
+ "discontinuous space were duplicated. Contact the developer "
+ else:
+ dup_fd_nodes = set(unique_fd_nodes[counts > 1])
+ for mm_inode, fd_inode in enumerate(self._reordering_arr_fd2mm):
+ if fd_inode in dup_fd_nodes:
+ self._duplicate_nodes.setdefault(fd_inode, [])
+ self._duplicate_nodes[fd_inode].append(mm_inode)
# Store things that we need for *from_fspace*
self._ufl_element = ufl_elt
@@ -228,7 +255,9 @@ class FromFiredrakeConnection:
np.matmul(out_view, self._resampling_mat_fd2mm.T, out=out_view)
return out
- def from_meshmode(self, mm_field, out=None):
+ def from_meshmode(self, mm_field, out=None,
+ assert_fdrake_discontinuous=True,
+ continuity_tolerance=None):
"""
transport meshmode field from :attr:`to_discr` into an
appropriate firedrake function space.
@@ -238,15 +267,26 @@ class FromFiredrakeConnection:
:param out: If *None* then ignored, otherwise a :mod:`firedrake`
function of the right function space for the transported data
to be stored in.
+ :param assert_fdrake_discontinuous: If *True*,
+ disallows conversion to a continuous firedrake function space
+ (i.e. this function checks that ``self.from_fspace()`` is
+ discontinuous and raises a *ValueError* otherwise)
+ :param continuity_tolerance: If converting to a continuous firedrake
+ function space (i.e. if ``self.from_fspace()`` is continuous),
+ assert that at any two meshmode nodes corresponding to the
+ same firedrake node (meshmode is a discontinuous space, so this
+ situation will almost certainly happen), the function being transported
+ has values at most :param:`continuity_tolerance` distance
+ apart. If *None*, no checks are performed.
:return: a :mod:`firedrake` :class:`Function` holding the transported
data.
"""
- if self._ufl_element.family() == 'Lagrange':
- raise ValueError("Cannot convert functions from discontinuous "
- " space (meshmode) to continuous firedrake "
- " space (reference element family %s)."
- % type(self._ufl_element.family()))
+ if self._ufl_element.family() == 'Lagrange' \
+ and assert_fdrake_discontinuous:
+ raise ValueError("Trying to convert to continuous function space "
+ " with :param:`assert_fdrake_discontinuous` set "
+ " to *True*")
# make sure out is a firedrake function in an appropriate
# function space
if out is not None:
@@ -273,16 +313,47 @@ class FromFiredrakeConnection:
if len(out.dat.data.shape) == 1 and len(mm_field.shape) > 1:
mm_field = mm_field.reshape(mm_field.shape[1])
- # resample from nodes
+ # resample from nodes on reordered view. Have to do this in
+ # a bit of a roundabout way to be careful about duplicated
+ # firedrake nodes.
by_cell_field_view = self.to_discr.groups[0].view(mm_field)
- by_cell_out_view = self.to_discr.groups[0].view(out.dat.data.T)
- np.matmul(by_cell_field_view, self._resampling_mat_mm2fd.T,
- out=by_cell_out_view)
-
- # reorder data
if len(out.dat.data.shape) == 1:
- out.dat.data[:] = out.dat.data[self._reordering_arr_mm2fd]
+ reordered_outdata = out.dat.data[self._reordering_arr_fd2mm]
else:
- out.dat.data[:] = out.dat.data[self._reordering_arr_mm2fd, :]
+ reordered_outdata = out.dat.data.T[:, self._reordering_arr_fd2mm]
+ by_cell_reordered_view = self.to_discr.groups[0].view(reordered_outdata)
+ np.matmul(by_cell_field_view, self._resampling_mat_mm2fd.T,
+ out=by_cell_reordered_view)
+ out.dat.data[self._reordering_arr_fd2mm] = reordered_outdata.T
+
+ # Continuity checks
+ if self._ufl_element.family() == 'Lagrange' \
+ and continuity_tolerance is not None:
+ assert isinstance(continuity_tolerance, float)
+ assert continuity_tolerance >= 0
+ # Check each firedrake node which has been duplicated
+ # that all associated values are within the continuity
+ # tolerance
+ for fd_inode, duplicated_mm_nodes in \
+ six.iteritems(self._duplicate_nodes):
+ mm_inode = duplicated_mm_nodes[0]
+ # Make sure to compare using reordered_outdata not mm_field,
+ # because two meshmode nodes associated to the same firedrake
+ # nodes may have been resampled to distinct nodes on different
+ # elements. reordered_outdata has undone that resampling.
+ for dup_mm_inode in duplicated_mm_nodes[1:]:
+ if len(reordered_outdata.shape) > 1:
+ dist = la.norm(reordered_outdata[:, mm_inode]
+ - reordered_outdata[:, dup_mm_inode])
+ else:
+ dist = la.norm(reordered_outdata[mm_inode]
+ - reordered_outdata[dup_mm_inode])
+ if dist >= continuity_tolerance:
+ raise ValueError("Meshmode nodes %s and %s represent "
+ "the same firedrake node %s, but "
+ ":param:`mm_field`'s values are "
+ " %s > %s apart)"
+ % (mm_inode, dup_mm_inode, fd_inode,
+ dist, continuity_tolerance))
return out
diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py
index 470cd724dcc732a75d3c50f7c1a24d3dcfa62fb1..d9060b9bfeadf04d5e19a850893797ae5622f33b 100644
--- a/meshmode/interop/firedrake/mesh.py
+++ b/meshmode/interop/firedrake/mesh.py
@@ -267,25 +267,13 @@ def _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices,
Return the orientations of the mesh elements:
:param fdrake_mesh: A :mod:`firedrake` instance of :class:`MeshGeometry`
-
:param unflipped_group: A :class:`SimplexElementGroup` instance with
(potentially) some negatively oriented elements.
-
:param vertices: The vertex coordinates as a numpy array of shape
*(ambient_dim, nvertices)*
-
- :param normals: _Only_ used if :param:`mesh` is a 1-surface
- embedded in 2-space. In this case,
- - If *None* then
- all elements are assumed to be positively oriented.
- - Else, should be a list/array whose *i*th entry
- is the normal for the *i*th element (*i*th
- in :param:`mesh`*.coordinate.function_space()*'s
- :attribute:`cell_node_list`)
-
- :param no_normals_warn: If *True*, raises a warning
- if :param:`mesh` is a 1-surface embedded in 2-space
- and :param:`normals` is *None*.
+ :param normals: As described in the kwargs of :func:`import_firedrake_mesh`
+ :param no_normals_warn: As described in the kwargs of
+ :func:`import_firedrake_mesh`
:return: A numpy array, the *i*th element is > 0 if the *ith* element
is positively oriented, < 0 if negatively oriented.
@@ -337,7 +325,7 @@ def _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices,
# {{{ Mesh conversion
-def import_firedrake_mesh(fdrake_mesh):
+def import_firedrake_mesh(fdrake_mesh, **kwargs):
"""
Create a :mod:`meshmode` :class:`Mesh` from a :mod:`firedrake`
:class:`MeshGeometry` with the same cells/elements, vertices, nodes,
@@ -357,7 +345,14 @@ def import_firedrake_mesh(fdrake_mesh):
1, 2, or 3 and have co-dimension of 0 or 1.
It must use a simplex as a reference element.
- Finally, its ``coordinates`` attribute must have a function
+ 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 the keyword arguments below.
+
+ 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
@@ -366,6 +361,19 @@ def import_firedrake_mesh(fdrake_mesh):
verify this by looking at the ``[0]`` entry of
``fdrake_mesh.coordinates.function_space().finat_element.entity_dofs()``.
+ :Keyword Arguments:
+ * *normals*: _Only_ used if :param:`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 :param:`mesh`*.coordinate.function_space()*'s
+ :attribute:`cell_node_list`)
+ * *no_normals_warn: If *True* (the default), raises a warning
+ if :param:`fdrake_mesh` is a 1-surface embedded in 2-space
+ and :param:`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)
@@ -449,8 +457,11 @@ def import_firedrake_mesh(fdrake_mesh):
vertices = np.array([vertices[i] for i in range(len(vertices))]).T
# Use the vertices to compute the orientations and flip the group
- # FIXME : Allow for passing in normals/no normals warn
- orient = _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices)
+ normals = kwargs.get('normals', None)
+ no_normals_warn = kwargs.get('no_normals_warn', True)
+ orient = _get_firedrake_orientations(fdrake_mesh, unflipped_group, vertices,
+ normals=normals,
+ no_normals_warn=no_normals_warn)
from meshmode.mesh.processing import flip_simplex_element_group
group = flip_simplex_element_group(vertices, unflipped_group, orient < 0)
diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py
index 31e4a4ba8484424f859042c5a6a7a6cb2de7959f..7074b58e2fbe35b942a0e8857c8459f636ff9527 100644
--- a/test/test_firedrake_interop.py
+++ b/test/test_firedrake_interop.py
@@ -181,7 +181,7 @@ def test_function_transfer(ctx_factory,
def check_idempotency(fdrake_connection, fdrake_function):
"""
- Make sure fd->mm->fd and mm->fd->mm are identity for DG spaces
+ Make sure fd->mm->fd and mm->fd->mm are identity
"""
vdim = None
if len(fdrake_function.dat.data.shape) > 1:
@@ -191,7 +191,9 @@ def check_idempotency(fdrake_connection, fdrake_function):
# Test for idempotency fd->mm->fd
mm_field = fdrake_connection.from_firedrake(fdrake_function)
fdrake_function_copy = Function(fdrake_fspace)
- fdrake_connection.from_meshmode(mm_field, fdrake_function_copy)
+ fdrake_connection.from_meshmode(mm_field, fdrake_function_copy,
+ assert_fdrake_discontinuous=False,
+ continuity_tolerance=1e-8)
np.testing.assert_allclose(fdrake_function_copy.dat.data,
fdrake_function.dat.data,
@@ -202,11 +204,12 @@ def check_idempotency(fdrake_connection, fdrake_function):
np.testing.assert_allclose(mm_field_copy, mm_field, atol=CLOSE_ATOL)
-def test_scalar_idempotency(ctx_factory, fdrake_mesh, fdrake_degree):
+def test_scalar_idempotency(ctx_factory, fdrake_mesh,
+ fdrake_family, fdrake_degree):
"""
- Make sure fd->mm->fd and mm->fd->mm are identity for scalar DG spaces
+ Make sure fd->mm->fd and mm->fd->mm are identity for scalar spaces
"""
- fdrake_fspace = FunctionSpace(fdrake_mesh, 'DG', fdrake_degree)
+ fdrake_fspace = FunctionSpace(fdrake_mesh, fdrake_family, fdrake_degree)
# Make a function with unique values at each node
fdrake_unique = Function(fdrake_fspace)
@@ -218,11 +221,12 @@ def test_scalar_idempotency(ctx_factory, fdrake_mesh, fdrake_degree):
check_idempotency(fdrake_connection, fdrake_unique)
-def test_vector_idempotency(ctx_factory, fdrake_mesh, fdrake_degree):
+def test_vector_idempotency(ctx_factory, fdrake_mesh,
+ fdrake_family, fdrake_degree):
"""
- Make sure fd->mm->fd and mm->fd->mm are identity for vector DG spaces
+ Make sure fd->mm->fd and mm->fd->mm are identity for vector spaces
"""
- fdrake_vfspace = VectorFunctionSpace(fdrake_mesh, 'DG', fdrake_degree)
+ fdrake_vfspace = VectorFunctionSpace(fdrake_mesh, fdrake_family, fdrake_degree)
# Make a function with unique values at each node
xx = SpatialCoordinate(fdrake_vfspace.mesh())