Merge branch 'lzlarryli/regge-element-changes'

FIAT/reference_element.py

@@ -16,6 +16,7 @@
 # along with FIAT. If not, see <http://www.gnu.org/licenses/>.
 #
 # Modified by David A. Ham (david.ham@imperial.ac.uk), 2014
+# Modified by Lizao Li (lzlarryli@gmail.com), 2016
 """
 Abstract class and particular implementations of finite element
 reference simplex geometry/topology.
@@ -305,6 +306,20 @@ def compute_face_tangents(self, face_i):
                                 self.get_vertices_of_subcomplex(t[2][face_i])))
         return (v1 - v0, v2 - v0)
 
+    def compute_face_edge_tangents(self, dim, entity_id):
+        """Computes all the edge tangents of any k-face with k>=1.
+        The result is a array of binom(dim+1,2) vectors.
+        This agrees with `compute_edge_tangent` when dim=1.
+        """
+        vert_ids = self.get_topology()[dim][entity_id]
+        vert_coords = [numpy.array(x)
+                       for x in self.get_vertices_of_subcomplex(vert_ids)]
+        edge_ts = []
+        for source in range(dim):
+            for dest in range(source + 1, dim + 1):
+                edge_ts.append(vert_coords[dest] - vert_coords[source])
+        return edge_ts
+
     def make_lattice(self, n, interior=0):
         """Constructs a lattice of points on the simplex.  For
         example, the 1:st order lattice will be just the vertices.

FIAT/regge.py

@@ -1,4 +1,7 @@
-# Copyright (C) 2015-2017 Lizao Li
+# -*- coding: utf-8 -*-
+"""Implementation of the generalized Regge finite elements."""
+
+# Copyright (C) 2015-2018 Lizao Li
 #
 # This file is part of FIAT.
 #
@@ -14,144 +17,92 @@
 #
 # 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 __future__ import absolute_import, print_function, division
 
-import numpy
-
 from FIAT.finite_element import FiniteElement
 from FIAT.dual_set import DualSet
 from FIAT.polynomial_set import ONSymTensorPolynomialSet
 from FIAT.functional import PointwiseInnerProductEvaluation as InnerProduct
-from FIAT.functional import index_iterator
 
 
 class ReggeDual(DualSet):
-
+    """Degrees of freedom for generalized Regge finite elements."""
     def __init__(self, cell, degree):
-        (dofs, ids) = self.generate_degrees_of_freedom(cell, degree)
-        super(ReggeDual, self).__init__(dofs, cell, ids)
-
-    def generate_degrees_of_freedom(self, cell, degree):
-        """
-        Suppose f is a k-face of the reference n-cell. Let t1,...,tk be a
-        basis for the tangent space of f as n-vectors. Given a symmetric
-        2-tensor field u on Rn. One set of dofs for Regge(r) on f is
-        the moment of each of the (k+1)k/2 scalar functions
-          [u(t1,t1),u(t1,t2),...,u(t1,tk),
-           u(t2,t2), u(t2,t3),...,..., u(tk,tk)]
-        aginst scalar polynomials of degrees (r-k+1). Here this is
-        implemented as pointwise evaluations of those scalar functions.
-
-        Below is an implementation for dimension 2--3. In dimension 1,
-        Regge(r)=DG(r). It is awkward in the current FEniCS interface to
-        implement the element uniformly for all dimensions due to its edge,
-        facet=trig, cell style.
-        """
-
-        dofs = []
-        ids = {}
-
-        d = cell.get_spatial_dimension()
-        if (d < 2) or (d > 3):
-            raise ValueError("Regge elements only implemented for dimension 2--3.")
-
-        # No vertex dof
-        ids[0] = dict(list(zip(list(range(d + 1)), ([] for i in range(d + 1)))))
+        dim = cell.get_spatial_dimension()
+        if (dim < 2) or (dim > 3):
+            raise ValueError("Generalized Regge elements are implemented only "
+                             "for dimension 2--3. For 1D, it is just DG(r).")
+
+        # construct the degrees of freedoms
+        dofs = []               # list of functionals
+        # dof_ids[i][j] contains the indices of dofs that are associated with
+        # entity j in dim i
+        dof_ids = {}
+
+        # no vertex dof
+        dof_ids[0] = {i: [] for i in range(dim + 1)}
         # edge dofs
-        (_dofs, _ids) = self._generate_edge_dofs(cell, degree, 0)
+        (_dofs, _dof_ids) = self._generate_dofs(cell, 1, degree, 0)
         dofs.extend(_dofs)
-        ids[1] = _ids
+        dof_ids[1] = _dof_ids
         # facet dofs for 3D
-        if d == 3:
-            (_dofs, _ids) = self._generate_facet_dofs(cell, degree, len(dofs))
+        if dim == 3:
+            (_dofs, _dof_ids) = self._generate_dofs(cell, 2, degree, len(dofs))
             dofs.extend(_dofs)
-            ids[2] = _ids
-        # Cell dofs
-        (_dofs, _ids) = self._generate_cell_dofs(cell, degree, len(dofs))
+            dof_ids[2] = _dof_ids
+        # cell dofs
+        (_dofs, _dof_ids) = self._generate_dofs(cell, dim, degree, len(dofs))
         dofs.extend(_dofs)
-        ids[d] = _ids
-        return (dofs, ids)
+        dof_ids[dim] = _dof_ids
 
-    def _generate_edge_dofs(self, cell, degree, offset):
-        """Generate dofs on edges."""
-        dofs = []
-        ids = {}
-        for s in range(len(cell.get_topology()[1])):
-            # Points to evaluate the inner product
-            pts = cell.make_points(1, s, degree + 2)
-            # Evalute squared length of the tagent vector along an edge
-            t = cell.compute_edge_tangent(s)
-            # Fill dofs
-            dofs += [InnerProduct(cell, t, t, p) for p in pts]
-            # Fill ids
-            i = len(pts) * s
-            ids[s] = list(range(offset + i, offset + i + len(pts)))
-        return (dofs, ids)
+        super(ReggeDual, self).__init__(dofs, cell, dof_ids)
 
-    def _generate_facet_dofs(self, cell, degree, offset):
-        """Generate dofs on facets in 3D."""
-        # Return empty if there is no such dofs
-        dofs = []
-        ids = dict(list(zip(list(range(4)), ([] for i in range(4)))))
-        if degree == 0:
-            return (dofs, ids)
-        # Compute dofs
-        for s in range(len(cell.get_topology()[2])):
-            # Points to evaluate the inner product
-            pts = cell.make_points(2, s, degree + 2)
-            # Let t1 and t2 be the two tangent vectors along a triangle
-            # we evaluate u(t1,t1), u(t1,t2), u(t2,t2) at each point.
-            (t1, t2) = cell.compute_face_tangents(s)
-            # Fill dofs
-            for p in pts:
-                dofs += [InnerProduct(cell, t1, t1, p),
-                         InnerProduct(cell, t1, t2, p),
-                         InnerProduct(cell, t2, t2, p)]
-            # Fill ids
-            i = len(pts) * s * 3
-            ids[s] = list(range(offset + i, offset + i + len(pts) * 3))
-        return (dofs, ids)
+    @staticmethod
+    def _generate_dofs(cell, entity_dim, degree, offset):
+        """generate degrees of freedom for enetities of dimension entity_dim
+
+        Input: all obvious except
+           offset  -- the current first available dof id.
 
-    def _generate_cell_dofs(self, cell, degree, offset):
-        """Generate dofs for cells."""
-        # Return empty if there is no such dofs
+        Output:
+           dofs    -- an array of dofs associated to entities in that dim
+           dof_ids -- a dict mapping entity_id to the range of indices of dofs
+                      associated to it.
+
+        On a k-face for degree r, the dofs are given by the value of
+           t^T u t
+        evaluated at points enough to control P(r-k+1) for all the edge
+        tangents of the face.
+        `cell.make_points(entity_dim, entity_id, degree + 2)` happens to
+        generate exactly those points needed.
+        """
         dofs = []
-        d = cell.get_spatial_dimension()
-        if (d == 2 and degree == 0) or (d == 3 and degree <= 1):
-            return ([], {0: []})
-        # Compute dofs. There is only one cell. So no need to loop here~
-        # Points to evaluate the inner product
-        pts = cell.make_points(d, 0, degree + 2)
-        # Let {e1,..,ek} be the Euclidean basis. We evaluate inner products
-        #   u(e1,e1), u(e1,e2), u(e1,e3), u(e2,e2), u(e2,e3),..., u(ek,ek)
-        # at each point.
-        e = numpy.eye(d)
-        # Fill dofs
-        for p in pts:
-            dofs += [InnerProduct(cell, e[i], e[j], p)
-                     for [i, j] in index_iterator((d, d)) if i <= j]
-        # Fill ids
-        ids = {0: list(range(offset, offset + len(pts) * d * (d + 1) // 2))}
-        return (dofs, ids)
+        dof_ids = {}
+        num_entities = len(cell.get_topology()[entity_dim])
+        for entity_id in range(num_entities):
+            pts = cell.make_points(entity_dim, entity_id, degree + 2)
+            tangents = cell.compute_face_edge_tangents(entity_dim, entity_id)
+            dofs += [InnerProduct(cell, t, t, pt)
+                     for pt in pts
+                     for t in tangents]
+            num_new_dofs = len(pts) * len(tangents)
+            dof_ids[entity_id] = list(range(offset, offset + num_new_dofs))
+            offset += num_new_dofs
+        return (dofs, dof_ids)
 
 
 class Regge(FiniteElement):
-    """The Regge elements on triangles and tetrahedra: the polynomial
-    space described by the full polynomials of degree k with degrees
-    of freedom to ensure its pullback as a metric to each interior
-    facet and edge is single-valued.
-
+    """The generalized Regge elements for symmetric-matrix-valued functions.
+       REG(r) in dimension n is the space of polynomial symmetric-matrix-valued
+       functions of degree r or less with tangential-tangential continuity.
     """
-
     def __init__(self, cell, degree):
-        # Check degree
         assert degree >= 0, "Regge start at degree 0!"
-        # Construct polynomial basis for d-vector fields
+        # shape functions
         Ps = ONSymTensorPolynomialSet(cell, degree)
-        # Construct dual space
+        # degrees of freedom
         Ls = ReggeDual(cell, degree)
-        # Set mapping
-        mapping = "pullback as metric"
-        # Call init of super-class
+        # mapping under affine transformation
+        mapping = "double covariant piola"
+
         super(Regge, self).__init__(Ps, Ls, degree, mapping=mapping)

doc/sphinx/source/releases/next.rst

@@ -11,6 +11,17 @@ Summary of changes
           be published (and renamed) to list the most important
           changes in the new release.
 
+- More elegant edge-based degrees of freedom are used for generalized Regge
+  finite elements.  This is a internal change and is not visible to other parts
+  of FEniCS.
+- The name of the mapping for generalized Regge finite element is changed to
+  "double covariant piola" from "pullback as metric". Geometrically, this
+  mapping is just the pullback of covariant 2-tensor fields in terms of proxy
+  matrix-fields. Because the mapping for 1-forms in FEniCS is currently named
+  "covariant piola", this mapping for symmetric tensor product of 1-forms is
+  thus called "double covariant piola". This change causes multiple internal
+  changes downstream in UFL and FFC. But this change should not be visible to
+  the end-user.
 
 Detailed changes
 ================

test/regression/fiat-reference-data-id

@@ -1 +1 @@
-2edbbad2000f1a15c775824fa90193c5d411ad08
+62da32c6c3ecde7d73436e4829f8832969c77519

test/regression/test_regression.py

@@ -266,7 +266,8 @@ def _perform_test(family, dim, degree, reference_table):
             assert eval(dtuple) in table
             assert table[eval(dtuple)].shape == reference_table[dtuple].shape
             diff = table[eval(dtuple)] - reference_table[dtuple]
-            assert (abs(diff) < tolerance).all()
+            assert (abs(diff) < tolerance).all(), \
+                "quadrature case %s %s %s failed!" % (family, dim, degree)
 
     filename = os.path.join(ref_path, "reference.json")