Interoperability with Other Discretization Packages

Functionality in this subpackage helps import and export data to/from other pieces of software, typically PDE solvers.

Nodal DG

Provides interoperability with the Matlab/Octave Codes complementing the book “Nodal Discontinuous Galerkin Methods” by Jan Hesthaven and Tim Warburton (Springer, 2008).

class meshmode.interop.nodal_dg.NodalDGContext(path)

Should be used as a context manager to ensure proper cleanup.

__init__(path)

Parameters:

path – The path to the Codes1.1 folder of the nodal DG codes.

set_mesh(mesh: Mesh, order)

Set the mesh information in the nodal DG Octave instance to the one given by mesh.

The mesh must only have a single element group of simplices.

Warning

High-order geometry information is currently silently ignored.

get_discr(actx) → Discretization

Get a discretization with nodes exactly matching the ones used by the nodal-DG code.

The returned discretization contains a new Mesh object constructed from the global Octave state.

push_dof_array(name, ary: DOFArray)

pull_dof_array(actx: ArrayContext, name) → DOFArray

meshmode.interop.nodal_dg.download_nodal_dg_if_not_present(path='nodal-dg')

Download the nodal-DG source code.

Parameters:

path – The destination path.

Firedrake

Function Spaces/Discretizations

Users wishing to interact with meshmode from firedrake will create a FiredrakeConnection using build_connection_from_firedrake(), while users wishing to interact with firedrake from meshmode will use will create a FiredrakeConnection using build_connection_to_firedrake(). It is not recommended to create a FiredrakeConnection directly.

class meshmode.interop.firedrake.connection.FiredrakeConnection(discr, fdrake_fspace, mm2fd_node_mapping, group_nr=None)

A connection between one group of a meshmode discretization and a firedrake “DG” function space.

Users should instantiate this using build_connection_to_firedrake() or build_connection_from_firedrake()

discr

A meshmode.discretization.Discretization.

group_nr

The group number identifying which element group of discr is being connected to a firedrake function space

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 firedrake node index associated to the jth degree of freedom of the ith element in element_grp.

mm2fd_node_mapping must encode an embedding into the firedrake mesh, i.e. no two meshmode nodes may be associated to the same 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 meshmode.interop.firedrake.reference_cell to

obtain firedrake unit dofs on the meshmode reference cell)

__init__(discr, fdrake_fspace, mm2fd_node_mapping, group_nr=None)

Parameters:

• discr – A meshmode.discretization.Discretization

• fdrake_fspace – A firedrake.functionspaceimpl.WithGeometry. Must use ufl family "Discontinuous Lagrange".

• mm2fd_node_mapping – Used as attribute mm2fd_node_mapping. A 2-D numpy integer array with the same dtype as fdrake_fspace.cell_node_list.dtype

• group_nr – The index of the group in discr which is being connected to fdrake_fspace. The group must be a 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.

• 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.

from_meshmode(mm_field, out=None)

Transport meshmode field from 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.

Parameters:

• mm_field –

  Either

  • A DOFArray representing a field of shape tuple() on discr

  • A numpy.ndarray of dtype “object” with entries of class DOFArray representing a field of shape mm_field.shape on discr

  See DOFArray for further requirements. The group_nr entry of each DOFArray must be of shape (nelements, nunit_dofs) and the element_dtype must match that used for firedrake.function.Functions

• out – If None then ignored, otherwise a 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

Returns:

a firedrake.function.Function holding the transported data (out, if out was not None)

from_firedrake(function, out=None, actx=None)

Transport firedrake function onto discr

Parameters:

• function – A firedrake.function.Function to transfer onto discr. Its function space must have the same family, degree, and mesh as self.from_fspace().

• out –

  Either

  1. A DOFArray

  2. A numpy.ndarray object array, each of whose entries is a 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 DOFArray must be a DOFArray defined on discr (as described in the documentation for DOFArray). Also, each 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 discr and stored in out, which is then returned.

• actx –

  • If out is None, then actx is a ArrayContext on which to create the DOFArray

  • If out is not None, actx must be None or out’s array_context.

Returns:

out, with the converted data from function stored in it

firedrake_fspace(shape=None)

Return a firedrake function space that self.discr.groups[self.group_nr] is connected to of the appropriate vector dimension

Parameters:

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

Returns:

A 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

meshmode.interop.firedrake.connection.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.

Parameters:

• discr – A Discretization to initialize the connection with

• 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 InterpolatoryQuadratureSimplexElementGroup.

• comm – Communicator to build a dmplex object on for the created firedrake mesh

meshmode.interop.firedrake.connection.build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None, restrict_to_boundary=None)

Create a FiredrakeConnection from a firedrake "DG" function space by creates a corresponding meshmode discretization and facilitating transfer of functions to and from firedrake.

Parameters:

• actx – A ArrayContext used to instantiate FiredrakeConnection.discr.

• fdrake_fspace – A firedrake "DG" function space (of class WithGeometry) built on a mesh which is importable by import_firedrake_mesh().

• grp_factory – (optional) If not None, should be a ElementGroupFactory whose group class is a subclass of InterpolatoryElementGroupBase. If None, and a default factory is automatically selected.

• restrict_to_boundary –

  (optional) If not None, then must be one of the following:

  • A valid boundary marker for fdrake_fspace.mesh()

  • A tuple of valid boundary markers for fdrake_fspace.mesh()

  • The string "on_boundary" (firedrake equivalent to BTAG_ALL).

  If not None, creates a 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.

Meshes

meshmode.interop.firedrake.mesh.import_firedrake_mesh(fdrake_mesh, cells_to_use=None, normals=None, no_normals_warn=None)

Create a meshmode.mesh.Mesh from a 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, firedrake does not require elements to be positively oriented, so any negative elements are flipped.

The flipped cells/elements are identified by the returned firedrake_orient array

Parameters:

• fdrake_mesh –

  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 firedrake, you can verify this by looking at

  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
  

• cells_to_use –

  cells_to_use is primarily intended for use internally by 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 BTAG_INDUCED_BOUNDARY.

• 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 ith entry is the normal for the ith element (ith in mesh.coordinate.function_space()’s cell_node_list)

• 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.

Returns:

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).

meshmode.interop.firedrake.mesh.export_mesh_to_firedrake(mesh, group_nr=None, comm=None)

Create a firedrake mesh corresponding to one Mesh’s SimplexElementGroup.

Parameters:

• mesh – A Mesh to convert with at least one SimplexElementGroup. ‘mesh.is_conforming’ must evaluate to True. ‘mesh’ must have vertices supplied, i.e. ‘mesh.vertices’ must not be None.

• group_nr – The group number to be converted into a firedrake mesh. The corresponding group must be of type SimplexElementGroup. If None and mesh has only one group, that group is used. Otherwise, a ValueError is raised.

• comm – The communicator to build the dmplex mesh on

Returns:

A tuple (fdrake_mesh, fdrake_cell_ordering, perm2cell) where

• fdrake_mesh is a firedrake firedrake.mesh.MeshGeometry corresponding to mesh

• fdrake_cell_ordering is a numpy array whose ith element in mesh (i.e. the ith element in mesh.groups[group_nr].vertex_indices) corresponds to the fdrake_cell_ordering[i]th 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 ith local vertex of the cth meshmode element.

Warning

Currently, no custom boundary tags are exported along with the mesh. firedrake seems to only allow one marker on each facet, whereas meshmode allows many.

Reference Cells

meshmode.interop.firedrake.reference_cell.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.

Parameters:

• ambient_dim – The spatial dimension

• 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, firedrake uses the standard FIAT simplex and meshmode uses modepy’s unit coordinates.

Returns:

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.

meshmode.interop.firedrake.reference_cell.get_finat_element_unit_nodes(finat_element)

Returns the unit nodes used by the finat element in firedrake’s (equivalently, finat/FIAT’s) reference coordinates

Parameters:

finat_element –

An instance of one of the following finat elements

• finat.fiat_elements.Lagrange

• finat.fiat_elements.DiscontinuousLagrange

• finat.fiat_elements.CrouzeixRaviart

• finat.spectral.GaussLobattoLegendre

• finat.spectral.GaussLegendre

Returns:

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)

Implementation Details

Converting between firedrake and meshmode is in general straightforward. Some language is different:

• In a mesh, a meshmode “element” is a firedrake “cell”

• A Discretization is a firedrake WithGeometry, usually created by calling the function FunctionSpace() and referred to as a “function space”

• In a mesh, any vertices, faces, cells, etc. are firedrake “entities” (see the PETSc documentation on DMPLEX for more info on how topological mesh information is stored in firedrake).

Other than carefully tabulating how and which vertices/faces correspond to other vertices/faces/cells, there are two main difficulties.

1. meshmode has discontinuous polynomial function spaces which may use different unit nodes than firedrake.

2. meshmode requires that all mesh elements be positively oriented, firedrake does not. Meanwhile, when firedrake creates a mesh, it changes the element ordering and the local vertex ordering.

(1.) is easily handled by insisting that the firedrake WithGeometry uses polynomial elements and that the group of the Discretization being converted is a InterpolatoryQuadratureSimplexElementGroup of the same order. Then, on each element, the function space being represented is the same in firedrake and meshmode. We may simply resample to one system or another’s unit nodes.

To handle (2.), once we associate a meshmode element to the correct firedrake cell, we have something like this picture:

(Assume we have already ensured that meshmode and firedrake use the same reference element by mapping firedrake’s reference element onto meshmode’s). If \(f=g\), then we can resample function values from one node set to the other. However, if firedrake has reordered the vertices or if we flipped their order to ensure meshmode has positively-oriented elements, there is some map \(A\) applied to the reference element which implements this permutation of barycentric coordinates. In this case, \(f=g\circ A\). Now, we have a connected diagram:

In short, once we reorder the firedrake nodes so that the mapping from the meshmode and 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 meshmode to firedrake). The information for this whole reordering process is stored in mm2fd_node_mapping, an array which associates each meshmode node to the firedrake node found by tracing the above diagram (i.e. it stores \(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 meshmode and 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 finat and FIAT is used to define what meshmode calls the unit nodes and unit vertices. It is worth noting that firedrake does not require a positive orientation of elements and that its reference triangle is different than specified in modepy.

2. A ~firedrake.mesh.MeshTopology which holds information about connectivity and other topological properties, but nothing about geometry/coordinates etc.

3. A class FunctionSpace created from a finat element and a ~firedrake.mesh.MeshTopology which allows us to define functions mapping the nodes (defined by the finat element) of each element in the ~firedrake.mesh.MeshTopology to some values. Note that the function FunctionSpace() in the firedrake API is used to create objects of class FunctionSpace and WithGeometry (see (6)).

4. A ~firedrake.function.CoordinatelessFunction (in the sense that its domain has no coordinates) which is a function in a FunctionSpace.

5. A ~firedrake.mesh.MeshGeometry created from a FunctionSpace and a ~firedrake.function.CoordinatelessFunction in that FunctionSpace which maps each dof to its geometric coordinates.

6. A WithGeometry which is a FunctionSpace together with a ~firedrake.mesh.MeshGeometry. This is the object returned usually returned to the user by a call to the firedrake function FunctionSpace().

7. A Function is defined on a WithGeometry.

Thus, by the coordinates of a mesh geometry we mean

1. On the hidden back-end: a ~firedrake.function.CoordinatelessFunction f on some function space defined only on the mesh topology.

2. On the front-end: A Function with the values of f but defined on a WithGeometry created from the FunctionSpace f lives in and the ~firedrake.mesh.MeshGeometry f defines.

Basically, it’s this picture (where \(a\to b\) if \(b\) depends on \(a\))

Warning

In general, the FunctionSpace of the coordinates function of a WithGeometry may not be the same FunctionSpace as for functions which live in the WithGeometry. This picture only shows how the class definitions depend on each other.