diff --git a/docs/source/manual.rst b/docs/source/manual.rst index fbdae8db..73635128 100644 --- a/docs/source/manual.rst +++ b/docs/source/manual.rst @@ -10,6 +10,7 @@ Manual manual/groups manual/serialisation manual/triples + manual/mixed manual/spaces diff --git a/docs/source/manual/mixed.rst b/docs/source/manual/mixed.rst new file mode 100644 index 00000000..fb037987 --- /dev/null +++ b/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. diff --git a/fuse/__init__.py b/fuse/__init__.py index 325213b9..a92f315f 100644 --- a/fuse/__init__.py +++ b/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 diff --git a/fuse/mixed.py b/fuse/mixed.py new file mode 100644 index 00000000..f52240e6 --- /dev/null +++ b/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"]) diff --git a/fuse/serialisation.py b/fuse/serialisation.py index 938fec87..ac4bb141 100644 --- a/fuse/serialisation.py +++ b/fuse/serialisation.py @@ -32,6 +32,7 @@ def __init__(self): "Edge": Edge, "Triple": ElementTriple, "VectorTriple": VectorTriple, + "MixedTriple": MixedTriple, "Group": GroupRepresentation, "GroupMember": GroupMemberRep, "PermutationSet": PermutationSetRepresentation, diff --git a/test/test_mixed.py b/test/test_mixed.py new file mode 100644 index 00000000..72f58800 --- /dev/null +++ b/test/test_mixed.py @@ -0,0 +1,175 @@ +import pytest +import numpy as np +import sympy as sp +from firedrake import * +from fuse import * +from FIAT.mixed import MixedElement as FIATMixedElement +from FIAT.quadrature_schemes import create_quadrature +from test_2d_examples_docs import construct_rt +from test_3d_examples_docs import construct_tet_rt +from test_convert_to_fiat import create_cg1, create_cg2_tri, create_dg1_tet +from fuse.mixed import MixedTriple +from fuse.vectorisation import VectorTriple + + +def taylor_hood_like(cell): + """Vector RT velocity block + scalar CG1 pressure block (mixed, block value shape).""" + return MixedTriple(construct_rt(cell), create_cg1(cell)) + + +def test_value_shape_and_num_dofs(): + cell = polygon(3) + rt = construct_rt(cell) + cg1 = create_cg1(cell) + me = MixedTriple(rt, cg1) + + # value shape is the flattened sum of the sub-element value sizes + assert me.get_value_shape() == (int(np.prod(rt.get_value_shape())) + 1,) + assert me.get_value_shape() == (3,) + assert me.num_dofs() == rt.num_dofs() + cg1.num_dofs() + + +def test_requires_common_cell(): + tri = polygon(3) + quad = polygon(4) + with pytest.raises(ValueError): + MixedTriple(create_cg1(tri), create_cg1(quad)) + + +def test_entity_ids_offset_concatenation(): + cell = polygon(3) + rt = construct_rt(cell) + cg1 = create_cg1(cell) + me = MixedTriple(rt, cg1) + + flat = [i for dim in me.entity_ids for ent in me.entity_ids[dim] + for i in me.entity_ids[dim][ent]] + assert sorted(flat) == list(range(me.num_dofs())) + # the second block's ids are all offset above the first block's count + second_block_ids = [i for i in flat if i >= rt.num_dofs()] + assert len(second_block_ids) == cg1.num_dofs() + + +def test_to_ufl_is_mixed_with_per_block_pullback(): + cell = polygon(3) + me = taylor_hood_like(cell) + ue = me.to_ufl() + + assert type(ue).__name__ == "MixedElement" + assert ue.reference_value_shape == (3,) + sub_pullbacks = [type(s.pullback).__name__ for s in ue.sub_elements] + assert sub_pullbacks == ["ContravariantPiola", "IdentityPullback"] + + +def test_to_fiat_block_diagonal_tabulation(): + cell = polygon(3) + rt = construct_rt(cell) + cg1 = create_cg1(cell) + me = MixedTriple(rt, cg1) + + fe = me.to_fiat() + assert isinstance(fe, FIATMixedElement) + assert fe.value_shape() == (3,) + assert fe.space_dimension() == 6 + + ref = cell.to_fiat() + pts = create_quadrature(ref, 2).get_points() + npts = len(pts) + + mixed_tab = fe.tabulate(0, pts)[(0, 0)] + assert mixed_tab.shape == (6, 3, npts) + + rt_tab = rt.to_fiat().tabulate(0, pts)[(0, 0)].reshape(3, 2, npts) + cg_tab = cg1.to_fiat().tabulate(0, pts)[(0, 0)].reshape(3, 1, npts) + + # velocity block occupies dof rows 0:3, components 0:2 + assert np.allclose(mixed_tab[0:3, 0:2, :], rt_tab) + assert np.allclose(mixed_tab[0:3, 2:3, :], 0.0) + # pressure block occupies dof rows 3:6, component 2 + assert np.allclose(mixed_tab[3:6, 2:3, :], cg_tab) + assert np.allclose(mixed_tab[3:6, 0:2, :], 0.0) + + +def _mixed_l2_project(W, gvec, gscal): + sigma, u = TrialFunctions(W) + tau, v = TestFunctions(W) + a = (inner(sigma, tau) + inner(u, v)) * dx + L = (inner(gvec, tau) + inner(gscal, v)) * dx + w = Function(W) + solve(a == L, w, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'}) + return w + + +# The six relative orientations of the two tetrahedra sharing a face. +_TWO_TET_PERMS = [sp.combinatorics.Permutation(p) for p in + ([0, 1, 2, 3], [0, 2, 3, 1], [0, 3, 1, 2], + [0, 1, 3, 2], [0, 3, 2, 1], [0, 2, 1, 3])] + + +@pytest.mark.parametrize("perm", _TWO_TET_PERMS) +def test_two_tet_mixed_orientation(perm): + """A mixed RT x DG element must reproduce a field it can represent exactly + on every relative orientation of the two cells - the decisive check that the + per-block custom orientation is applied correctly through the mixed element.""" + from firedrake.utility_meshes import TwoTetMesh + cell = make_tetrahedron() + me = MixedTriple(construct_tet_rt(cell), create_dg1_tet(cell)) + + mesh = TwoTetMesh(perm=perm, use_fuse=True) + x = SpatialCoordinate(mesh) + gvec = as_vector([1.0, 2.0, 3.0]) # constant vector, reproduced by RT1 + gscal = x[0] # linear scalar, reproduced by DG1 + + W = FunctionSpace(mesh, me.to_ufl()) + sigma_h, u_h = _mixed_l2_project(W, gvec, gscal).subfunctions + + assert errornorm(gvec, sigma_h) < 1e-10 + assert errornorm(gscal, u_h) < 1e-10 + + +def taylor_hood(cell): + """FUSE Taylor-Hood element: vector CG2 velocity x scalar CG1 pressure.""" + return MixedTriple(VectorTriple(create_cg2_tri(cell)), create_cg1(cell)) + + +# Direct solver options robust to the pressure constant nullspace of a Stokes +# saddle-point system (mumps null-pivot detection), so the solve does not depend +# on the fragile default handling of the singular matrix. +_STOKES_PARAMS = {'mat_type': 'aij', 'ksp_type': 'gmres', 'ksp_rtol': 1e-13, + 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps', + 'mat_mumps_icntl_24': 1, 'mat_mumps_icntl_25': 0} + + +def test_taylor_hood_stokes(): + """Solve Stokes with the FUSE Taylor-Hood mixed element on a manufactured + solution that lies exactly in the discrete space, and check it is recovered. + + With u = (y, -x) (divergence free, degree 1) and p = x (degree 1), the forcing + is f = -div(grad(u)) + grad(p) = (1, 0). Taylor-Hood contains this solution + exactly, so the mixed element must reproduce it to solver tolerance.""" + cell = polygon(3) + mesh = UnitSquareMesh(4, 4, use_fuse=True) + x = SpatialCoordinate(mesh) + u_exact = as_vector([x[1], -x[0]]) + p_exact = x[0] + f = as_vector([1.0, 0.0]) + + W = FunctionSpace(mesh, taylor_hood(cell).to_ufl()) + u, p = TrialFunctions(W) + v, q = TestFunctions(W) + a = (inner(grad(u), grad(v)) - p * div(v) - q * div(u)) * dx + L = inner(f, v) * dx + bc = DirichletBC(W.sub(0), u_exact, "on_boundary") + nullspace = MixedVectorSpaceBasis( + W, [W.sub(0), VectorSpaceBasis(constant=True, comm=W.comm)]) + up = Function(W) + solve(a == L, up, bcs=bc, nullspace=nullspace, solver_parameters=_STOKES_PARAMS) + uh, ph = up.subfunctions + + # pressure is determined only up to a constant, so compare with the mean removed + area = assemble(Constant(1.0) * dx(domain=mesh)) + ph0 = ph - assemble(ph * dx) / area + pex0 = p_exact - assemble(p_exact * dx(domain=mesh)) / area + + assert sqrt(assemble(inner(uh - u_exact, uh - u_exact) * dx)) < 1e-10 + assert sqrt(assemble((ph0 - pex0) ** 2 * dx)) < 1e-10