Merge pull request #13 from firedrakeproject/discontinuous-lagrange-cube

Discontinuous Serendipity element for cuboids (DPC)

FIAT/__init__.py

@@ -15,6 +15,7 @@
 from FIAT.discontinuous_lagrange import DiscontinuousLagrange
 from FIAT.discontinuous_taylor import DiscontinuousTaylor
 from FIAT.discontinuous_raviart_thomas import DiscontinuousRaviartThomas
+from FIAT.discontinuous_pc import DPC
 from FIAT.hermite import CubicHermite
 from FIAT.lagrange import Lagrange
 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre
@@ -55,6 +56,7 @@
                       "FacetBubble": FacetBubble,
                       "Crouzeix-Raviart": CrouzeixRaviart,
                       "Discontinuous Lagrange": DiscontinuousLagrange,
+                      "DPC": DPC,
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
                       "Hermite": CubicHermite,

FIAT/discontinuous_pc.py

@@ -0,0 +1,115 @@
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+
+from FIAT import finite_element, polynomial_set, dual_set, functional
+from FIAT.reference_element import (Point,
+                                    DefaultLine,
+                                    UFCInterval,
+                                    UFCQuadrilateral,
+                                    UFCHexahedron,
+                                    UFCTriangle,
+                                    UFCTetrahedron,
+                                    make_affine_mapping)
+from FIAT.P0 import P0Dual
+import numpy as np
+
+hypercube_simplex_map = {Point(): Point(),
+                         DefaultLine(): DefaultLine(),
+                         UFCInterval(): UFCInterval(),
+                         UFCQuadrilateral(): UFCTriangle(),
+                         UFCHexahedron(): UFCTetrahedron()}
+
+
+class DPC0(finite_element.CiarletElement):
+    def __init__(self, ref_el):
+        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0)
+        dual = P0Dual(ref_el)
+        degree = 0
+        formdegree = ref_el.get_spatial_dimension()  # n-form
+        super(DPC0, self).__init__(poly_set=poly_set,
+                                   dual=dual,
+                                   order=degree,
+                                   ref_el=ref_el,
+                                   formdegree=formdegree)
+
+
+class DPCDualSet(dual_set.DualSet):
+    """The dual basis for DPC elements.  This class works for
+    hypercubes of any dimension.  Nodes are point evaluation at
+    equispaced points.  This is the discontinuous version where
+    all nodes are topologically associated with the cell itself"""
+
+    def __init__(self, ref_el, degree):
+        entity_ids = {}
+        nodes = []
+
+        # Change coordinates here.
+        # Vertices of the simplex corresponding to the reference element.
+        v_simplex = hypercube_simplex_map[ref_el].get_vertices()
+        # Vertices of the reference element.
+        v_hypercube = ref_el.get_vertices()
+        # For the mapping, first two vertices are unchanged in all dimensions.
+        v_ = list(v_hypercube[:2])
+
+        # For dimension 1 upwards,
+        # take the next vertex and map it to the midpoint of the edge/face it belongs to, and shares
+        # with no other points.
+        for d in range(1, ref_el.get_dimension()):
+            v_.append(tuple(np.asarray(v_hypercube[d+1] +
+                            np.average(np.asarray(v_hypercube[2**d:2**(d+1)]), axis=0))))
+        A, b = make_affine_mapping(v_simplex, tuple(v_))  # Make affine mapping to be used later.
+
+        # make nodes by getting points
+        # need to do this dimension-by-dimension, facet-by-facet
+        top = hypercube_simplex_map[ref_el].get_topology()
+        cube_topology = ref_el.get_topology()
+
+        cur = 0
+        for dim in sorted(top):
+            entity_ids[dim] = {}
+            for entity in sorted(top[dim]):
+                pts_cur = hypercube_simplex_map[ref_el].make_points(dim, entity, degree)
+                pts_cur = [tuple(np.matmul(A, np.array(x)) + b) for x in pts_cur]
+                nodes_cur = [functional.PointEvaluation(ref_el, x)
+                             for x in pts_cur]
+                nnodes_cur = len(nodes_cur)
+                nodes += nodes_cur
+                cur += nnodes_cur
+            for entity in sorted(cube_topology[dim]):
+                entity_ids[dim][entity] = []
+
+        entity_ids[dim][0] = list(range(len(nodes)))
+        super(DPCDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids)
+
+
+class HigherOrderDPC(finite_element.CiarletElement):
+    """The DPC finite element.  It is what it is."""
+
+    def __init__(self, ref_el, degree):
+        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], degree)
+        dual = DPCDualSet(ref_el, degree)
+        formdegree = ref_el.get_spatial_dimension()  # n-form
+        super(HigherOrderDPC, self).__init__(poly_set=poly_set,
+                                             dual=dual,
+                                             order=degree,
+                                             ref_el=ref_el,
+                                             formdegree=formdegree)
+
+
+def DPC(ref_el, degree):
+    if degree == 0:
+        return DPC0(ref_el)
+    else:
+        return HigherOrderDPC(ref_el, degree)

FIAT/finite_element.py

@@ -120,8 +120,8 @@ class CiarletElement(FiniteElement):
     basis generated from polynomials encoded in a `PolynomialSet`.
     """
 
-    def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine"):
-        ref_el = poly_set.get_reference_element()
+    def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref_el=None):
+        ref_el = ref_el or poly_set.get_reference_element()
         super(CiarletElement, self).__init__(ref_el, dual, order, formdegree, mapping)
 
         # build generalized Vandermonde matrix

FIAT/quadrature.py

@@ -252,6 +252,14 @@ def make_quadrature(ref_el, m):
         return CollapsedQuadratureTriangleRule(ref_el, m)
     elif ref_el.get_shape() == reference_element.TETRAHEDRON:
         return CollapsedQuadratureTetrahedronRule(ref_el, m)
+    elif ref_el.get_shape() == reference_element.QUADRILATERAL:
+        line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m)
+        return make_tensor_product_quadrature(line_rule, line_rule)
+    elif ref_el.get_shape() == reference_element.HEXAHEDRON:
+        line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m)
+        return make_tensor_product_quadrature(line_rule, line_rule, line_rule)
+    else:
+        raise ValueError("Unable to make quadrature for cell: %s" % ref_el)
 
 
 def make_tensor_product_quadrature(*quad_rules):

test/unit/test_discontinuous_pc.py

@@ -0,0 +1,51 @@
+# Copyright (C) 2016 Imperial College London and others
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+#
+# Authors:
+#
+# David Ham
+
+import pytest
+import numpy as np
+
+
+@pytest.mark.parametrize("dim, degree", [(dim, degree)
+                                         for dim in range(1, 4)
+                                         for degree in range(4)])
+def test_basis_values(dim, degree):
+    """Ensure that integrating a simple monomial produces the expected results."""
+    from FIAT import ufc_cell, make_quadrature
+    from FIAT.discontinuous_pc import DPC
+
+    cell = np.array([None, 'interval', 'quadrilateral', 'hexahedron'])
+    s = ufc_cell(cell[dim])
+    q = make_quadrature(s, degree + 1)
+
+    fe = DPC(s, degree)
+    tab = fe.tabulate(0, q.pts)[(0,) * dim]
+
+    for test_degree in range(degree + 1):
+        coefs = [n(lambda x: x[0]**test_degree) for n in fe.dual.nodes]
+        integral = np.float(np.dot(coefs, np.dot(tab, q.wts)))
+        reference = np.dot([x[0]**test_degree
+                            for x in q.pts], q.wts)
+        assert np.isclose(integral, reference, rtol=1e-14)
+
+
+if __name__ == '__main__':
+    import os
+    pytest.main(os.path.abspath(__file__))