add mixed element infrastructure

docs/source/manual.rst

@@ -10,6 +10,7 @@ Manual
     manual/groups
     manual/serialisation
     manual/triples
+    manual/mixed
     manual/spaces
 
 

docs/source/manual/mixed.rst

@@ -0,0 +1,95 @@
+Mixed elements
+==============
+
+A **mixed element** is the Cartesian product of a sequence of element triples
+defined on a common cell. It is the object used to discretise systems of
+equations in which different fields are approximated in different spaces, for
+example the velocity--pressure pair of the Stokes equations (Taylor--Hood,
+:math:`[P_2]^d \times P_1`) or the flux--potential pair of mixed Poisson
+(:math:`RT_1 \times DG_0`).
+
+Given sub-elements with dual bases
+:math:`\{\mathcal{X}^{(k)}_i\}` and value shapes :math:`s_k`, the mixed element
+has
+
+* dual basis the **disjoint union** :math:`\bigcup_k \{\mathcal{X}^{(k)}_i\}`
+  (the blocks are independent -- there is no coupling between fields at the
+  element level),
+* value shape the flattened **sum** :math:`\left(\sum_k \prod s_k\right)`, with
+  each sub-element occupying its own block of components,
+* the pullback of each block taken from that sub-element (identity, covariant or
+  contravariant Piola), applied per block.
+
+A mixed element is therefore **not** a Ciarlet element: it has no single
+polynomial space and no single nodal basis. In FUSE it is represented by the
+:class:`~fuse.mixed.MixedTriple` class rather than by
+:class:`~fuse.triples.ElementTriple`.
+
+This is distinct from **enrichment** (the direct sum of spaces sharing the *same*
+value shape and mapping, such as the :math:`P_1 \oplus \text{bubble}` velocity of
+the MINI element), which is expressed with addition of triples and yields an
+ordinary :class:`~fuse.triples.ElementTriple`.
+
+Construction
+------------
+
+A mixed element is built from any existing element triples on the same cell.
+Here ``velocity`` is a vector-valued :math:`H(\mathrm{div})` triple (a Raviart--Thomas
+element) and ``pressure`` is a scalar Lagrange triple, both constructed on the same
+cell as shown in :doc:`examples2d`::
+
+    from fuse.mixed import MixedTriple
+
+    mixed = MixedTriple(velocity, pressure)
+
+    mixed.get_value_shape()            # (3,)  = 2 (velocity) + 1 (pressure)
+    mixed.num_dofs()                   # sum of the sub-element dof counts
+
+The element converts to UFL/FInAT and FIAT as a single object::
+
+    ufl_element = mixed.to_ufl()       # finat.ufl.MixedElement
+    fiat_element = mixed.to_fiat()     # FIAT.mixed.MixedElement
+
+In Firedrake the resulting element defines a mixed function space directly::
+
+    W = FunctionSpace(mesh, mixed.to_ufl())
+
+A block may itself be vector valued. The Taylor--Hood element pairs a vector
+Lagrange velocity with a scalar Lagrange pressure, the velocity being obtained by
+vectorising a scalar triple with :class:`~fuse.vectorisation.VectorTriple`::
+
+    from fuse.vectorisation import VectorTriple
+
+    taylor_hood = MixedTriple(VectorTriple(cg2), cg1)   # [P2]^d x P1
+
+Each sub-element keeps its own mapping, so the contravariant Piola pullback of the
+:math:`RT` velocity and the identity pullback of the pressure are applied per
+block. FUSE's custom entity orientation is likewise applied independently to each
+block during assembly, so a mixed element built from FUSE sub-elements behaves
+exactly as those sub-elements do on their own.
+
+Composability
+-------------
+
+A mixed element **can** be combined with:
+
+* ``to_ufl`` / ``to_fiat`` conversion as a unit, and Firedrake assembly,
+  interpolation, projection and solves via ``FunctionSpace(mesh, mixed.to_ufl())``;
+* any FUSE element as a block -- Lagrange, DG, vector-valued, :math:`RT`,
+  Nédélec and BDM elements, and enriched triples formed by adding triples
+  (enabling, e.g., the MINI element);
+* plotting and introspection of the combined degrees of freedom.
+
+A mixed element **cannot** be:
+
+* treated as a Ciarlet element -- it has no single polynomial space, nodal basis
+  or well-defined single degree;
+* immersed onto a facet or passed to a trace (:class:`~fuse.traces.TrH1`,
+  :class:`~fuse.traces.TrHDiv`, ...): the sub-elements have heterogeneous Sobolev
+  spaces and no common trace. Immerse the sub-elements *before* mixing;
+* enriched (added) with another element -- enrichment requires an identical value
+  shape and mapping;
+* used as a factor in a tensor product.
+
+FUSE does not check the well-posedness (inf--sup / LBB condition) of the mixed
+system; that remains the responsibility of the user.

fuse/__init__.py

@@ -5,6 +5,7 @@
 from fuse.traces import TrH1, TrGrad, TrHess, TrHCurl, TrHDiv
 from fuse.tensor_products import tensor_product
 from fuse.vectorisation import VectorTriple
+from fuse.mixed import MixedTriple
 
 from fuse.spaces.element_sobolev_spaces import CellH1, CellL2, CellHDiv, CellHCurl, CellH2
 from fuse.spaces.polynomial_spaces import P0, P1, P2, P3, Q2, PolynomialSpace

fuse/mixed.py

@@ -0,0 +1,138 @@
+import numpy as np
+import finat.ufl
+from FIAT.mixed import MixedElement as FIATMixedElement
+from fuse.triples import ElementTriple
+from fuse.dof import DeltaPairing, L2Pairing, FuseFunction, PointKernel
+from fuse.traces import Trace
+
+
+class MixedTriple():
+    """
+    A mixed element: the Cartesian product of a sequence of element triples
+    defined on a common cell.
+
+    Unlike an :class:`ElementTriple`, a mixed element is not itself a Ciarlet
+    element - it has no single polynomial space or nodal basis, and so does not
+    subclass :class:`ElementTriple`. Its degrees of freedom are the disjoint union
+    of the sub-element functionals, its value shape is the sum of the sub-element
+    value sizes, and each sub-element retains its own mapping. The name follows the
+    ``VectorTriple`` / ``TensorProductTriple`` family and avoids clashing with
+    ``ufl.MixedElement``.
+
+    :param sub_elements: The element triples to combine. May be passed either as
+        separate arguments or as a single iterable.
+    """
+
+    def __init__(self, *sub_elements):
+        if len(sub_elements) == 1 and not isinstance(sub_elements[0], ElementTriple):
+            sub_elements = tuple(sub_elements[0])
+        if len(sub_elements) < 2:
+            raise ValueError("A mixed element requires at least two sub-elements")
+        for e in sub_elements:
+            if not isinstance(e, ElementTriple):
+                raise ValueError("Mixed element sub-elements must be ElementTriples")
+
+        def cell_signature(cell):
+            return (cell.get_spatial_dimension(), len(cell.vertices()))
+
+        ref = cell_signature(sub_elements[0].cell)
+        if not all(cell_signature(e.cell) == ref for e in sub_elements):
+            raise ValueError("Mixed element sub-elements must be defined on the same cell")
+
+        self._sub_elements = tuple(sub_elements)
+        self.cell = sub_elements[0].cell
+
+    @property
+    def sub_elements(self):
+        return list(self._sub_elements)
+
+    def __repr__(self):
+        return "MixedTriple(%s)" % ", ".join(repr(e) for e in self._sub_elements)
+
+    def _sub_value_size(self, e):
+        return int(np.prod(e.get_value_shape(), dtype=int))
+
+    def get_value_shape(self):
+        return (sum(self._sub_value_size(e) for e in self._sub_elements),)
+
+    def num_dofs(self):
+        return sum(e.num_dofs() for e in self._sub_elements)
+
+    @property
+    def entity_ids(self):
+        """Entity to dof-id map, offset-concatenated across the sub-elements."""
+        combined = None
+        offset = 0
+        for e in self._sub_elements:
+            e.to_ufl()
+            if combined is None:
+                combined = {dim: {ent: [] for ent in e.entity_ids[dim]}
+                            for dim in e.entity_ids}
+            for dim in e.entity_ids:
+                for ent in e.entity_ids[dim]:
+                    combined[dim][ent] += [i + offset for i in e.entity_ids[dim][ent]]
+            offset += e.num_dofs()
+        return combined
+
+    def to_ufl(self):
+        return finat.ufl.MixedElement(*[e.to_ufl() for e in self._sub_elements])
+
+    def to_fiat(self):
+        return FIATMixedElement([e.to_fiat() for e in self._sub_elements])
+
+    def generate(self):
+        dofs = []
+        for e in self._sub_elements:
+            dofs.extend(e.generate())
+        return dofs
+
+    def plot(self, filename="temp.png"):
+        import matplotlib.pyplot as plt
+        if self.cell.dimension == 0:
+            raise ValueError("Dimension 0 cells cannot be plotted")
+        if self.cell.dimension > 3:
+            raise ValueError("Plotting not supported in this dimension")
+
+        identity = FuseFunction(lambda *x: x)
+        fig = plt.figure()
+        if self.cell.dimension < 3:
+            ax = plt.gca()
+            self.cell.plot(show=False, plain=True, ax=ax)
+        else:
+            ax = fig.add_subplot(projection='3d')
+            self.cell.plot3d(show=False, ax=ax)
+
+        for block, e in enumerate(self._sub_elements):
+            for dof in e.generate():
+                center, color = e.get_dof_info(dof)
+                if center is None:
+                    center = [0, 0, 0]
+                if isinstance(dof.pairing, DeltaPairing) and isinstance(dof.kernel, PointKernel):
+                    coord = dof.eval(identity, pullback=False)
+                elif isinstance(dof.pairing, L2Pairing):
+                    coord = center
+                else:
+                    coord = center
+                if len(coord) == 1:
+                    coord = (coord[0], 0)
+                if isinstance(dof.target_space, Trace):
+                    dof.target_space.plot(ax, coord, dof.cell_defined_on, color=color)
+                else:
+                    ax.scatter(*coord, color=color)
+                ax.text(*coord, "%d.%d" % (block, dof.id))
+        plt.axis('off')
+        ax.get_xaxis().set_visible(False)
+        ax.get_yaxis().set_visible(False)
+        if filename:
+            fig.savefig(filename)
+        else:
+            plt.show()
+
+    def _to_dict(self):
+        return {"sub_elements": list(self._sub_elements)}
+
+    def dict_id(self):
+        return "MixedTriple"
+
+    def _from_dict(o_dict):
+        return MixedTriple(o_dict["sub_elements"])

fuse/serialisation.py

@@ -32,6 +32,7 @@ def __init__(self):
                           "Edge": Edge,
                           "Triple": ElementTriple,
                           "VectorTriple": VectorTriple,
+                          "MixedTriple": MixedTriple,
                           "Group": GroupRepresentation,
                           "GroupMember": GroupMemberRep,
                           "PermutationSet": PermutationSetRepresentation,