Div-decoupling dofs

FIAT/functional.py

@@ -264,11 +264,72 @@ def __init__(self, ref_el, facet_no, pt):
 
         self.tau = tau
         self.alphas = alphas
-        dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]}
+        dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(len(alphas))]}
 
         super().__init__(ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv")
 
 
+class PointTangentialDerivative(Functional):
+    """Represents d/dt at a point on an edge."""
+    def __init__(self, ref_el, edge_no, pt, comp=(), shp=()):
+        t = ref_el.compute_edge_tangent(edge_no)
+        sd = ref_el.get_spatial_dimension()
+
+        alphas = []
+        for i in range(sd):
+            alpha = [0] * sd
+            alpha[i] = 1
+            alphas.append(tuple(alpha))
+        dpt_dict = {tuple(pt): [(t[i], alphas[i], comp) for i in range(sd)]}
+
+        super().__init__(ref_el, shp, {}, dpt_dict, "PointTangentialDeriv")
+
+
+class PointTangentialSecondDerivative(Functional):
+    """Represents d^/dt^2 at a point on an edge."""
+    def __init__(self, ref_el, edge_no, pt, comp=(), shp=()):
+        t = ref_el.compute_edge_tangent(edge_no)
+        sd = ref_el.get_spatial_dimension()
+        tau = numpy.zeros((sd*(sd+1)//2,))
+
+        alphas = []
+        cur = 0
+        for i in range(sd):
+            for j in range(i, sd):
+                alpha = [0] * sd
+                alpha[i] += 1
+                alpha[j] += 1
+                alphas.append(tuple(alpha))
+                tau[cur] = t[i]*t[j] * (1 + (i != j))
+                cur += 1
+
+        dpt_dict = {tuple(pt): [(tau[i], alphas[i], comp) for i in range(len(alphas))]}
+
+        super().__init__(ref_el, shp, {}, dpt_dict, "PointTangentialSecondDeriv")
+
+
+class PointSecondDerivative(Functional):
+    """Represents d/ds1 d/ds2 at a point on an edge."""
+    def __init__(self, ref_el, s1, s2, pt, comp=(), shp=()):
+        sd = ref_el.get_spatial_dimension()
+        tau = numpy.zeros((sd*(sd+1)//2,))
+
+        alphas = []
+        cur = 0
+        for i in range(sd):
+            for j in range(i, sd):
+                alpha = [0] * sd
+                alpha[i] += 1
+                alpha[j] += 1
+                alphas.append(tuple(alpha))
+                tau[cur] = s1[i]*s2[j] / (1 + (i == j))
+                cur += 1
+
+        dpt_dict = {tuple(pt): [(tau[i], alphas[i], comp) for i in range(len(alphas))]}
+
+        super().__init__(ref_el, shp, {}, dpt_dict, "PointSecondDeriv")
+
+
 class PointDivergence(Functional):
     """Class representing point divergence of vector
     functions at a particular point x."""
@@ -313,6 +374,26 @@ def __call__(self, fn):
         return result
 
 
+class IntegralMomentOfDerivative(Functional):
+    """Functional giving directional derivative integrated against some function on a facet."""
+
+    def __init__(self, ref_el, s, Q, f_at_qpts, comp=(), shp=()):
+        self.f_at_qpts = f_at_qpts
+        self.Q = Q
+
+        sd = ref_el.get_spatial_dimension()
+
+        points = Q.get_points()
+        weights = numpy.multiply(f_at_qpts, Q.get_weights())
+
+        alphas = tuple(map(tuple, numpy.eye(sd, dtype=int)))
+        dpt_dict = {tuple(pt): [(wt*s[i], alphas[i], comp) for i in range(sd)]
+                    for pt, wt in zip(points, weights)}
+
+        super().__init__(ref_el, shp,
+                         {}, dpt_dict, "IntegralMomentOfDerivative")
+
+
 class IntegralMomentOfNormalDerivative(Functional):
     """Functional giving normal derivative integrated against some function on a facet."""
 

FIAT/stokes.py

@@ -8,8 +8,15 @@
 import scipy
 
 from FIAT import finite_element, dual_set
-from FIAT.functional import ComponentPointEvaluation, FrobeniusIntegralMoment
+from FIAT.functional import (ComponentPointEvaluation,
+                             PointTangentialDerivative,
+                             PointTangentialSecondDerivative,
+                             PointSecondDerivative,
+                             IntegralMomentOfDerivative,
+                             IntegralMoment,
+                             FrobeniusIntegralMoment)
 from FIAT.polynomial_set import make_bubbles, ONPolynomialSet
+from FIAT.expansions import polynomial_dimension
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.reference_element import symmetric_simplex
@@ -25,17 +32,10 @@ def eps(table_u):
     return numpy.transpose(eps_u, (1, 0, 2)).reshape(num_members, d, d, -1)
 
 
-def grad(table_u):
-    grad_u = {alpha.index(1): table_u[alpha] for alpha in table_u if sum(alpha) == 1}
-    grad_u = [grad_u[k] for k in sorted(grad_u)]
-    return numpy.transpose(grad_u, (1, 0, 2, 3))
-
-
 def inner(v, u, Qwts):
     return numpy.tensordot(numpy.multiply(v, Qwts), u, axes=(range(1, v.ndim), range(1, u.ndim)))
 
 
-
 def map_duals(ref_el, dim, entity, mapping, Q_ref, Phis):
     Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
     if mapping == "normal":
@@ -60,7 +60,8 @@ def __init__(self, ref_el, degree):
         entity_ids = {}
         top = ref_el.get_topology()
         sd = ref_el.get_spatial_dimension()
-        self.shp = (sd,)
+        self.sd = sd
+        shp = (sd,)
 
         mapping = None
         for dim in sorted(top):
@@ -69,73 +70,99 @@ def __init__(self, ref_el, degree):
                 for entity in sorted(top[dim]):
                     cur = len(nodes)
                     pts = ref_el.make_points(dim, entity, degree)
-                    nodes.extend(ComponentPointEvaluation(ref_el, (i,), self.shp, pt)
+                    nodes.extend(ComponentPointEvaluation(ref_el, (i,), shp, pt)
                                  for pt in pts for i in range(sd))
                     entity_ids[dim][entity] = list(range(cur, len(nodes)))
             else:
-                Q_ref, Phis = self._reference_duals(dim, degree)
+                # degree of bubbles
+                moment_degree = degree
+                if dim == 1:
+                    moment_degree -= 2*sd
+                elif dim == 2 and sd == 3:
+                    moment_degree -= 6
+                Q_ref, Phis = self._reference_duals(dim, degree, moment_degree)
+
                 for entity in sorted(top[dim]):
                     cur = len(nodes)
+                    if dim == 1:
+                        verts = ref_el.get_vertices_of_subcomplex(top[dim][entity])
+                        ells = [PointTangentialDerivative]
+                        if sd == 3:
+                            ells.append(PointTangentialSecondDerivative)
+                        nodes.extend(ell(ref_el, entity, pt, comp=(i,), shp=shp)
+                                     for pt in verts for ell in ells for i in range(sd))
+
+                    elif dim == 2 and sd == 3:
+                        # Face vertex dofs
+                        verts = numpy.array(ref_el.get_vertices_of_subcomplex(top[dim][entity]))
+                        for i in range(len(verts)):
+                            tangents = [verts[j] - verts[i] for j in range(len(verts)) if j != i]
+                            nodes.extend(PointSecondDerivative(ref_el, *tangents, verts[i], comp=(k,), shp=shp)
+                                         for k in range(sd))
+
+                        # Face edge dofs
+                        face = set(top[dim][entity])
+                        edges = [e for e in top[1] if set(top[1][e]) < face]
+
+                        phi_degree = degree - 5
+                        ref_edge = ref_el.construct_subelement(1)
+                        Q_edge = create_quadrature(ref_edge, degree + phi_degree)
+                        Phis_edge = ONPolynomialSet(ref_edge, phi_degree).tabulate(Q_edge.get_points())[(0,)]
+                        n = ref_el.compute_scaled_normal(entity)
+                        for e in edges:
+                            t = ref_el.compute_edge_tangent(e)
+                            s = numpy.cross(n, t)
+                            Q, phis = map_duals(ref_el, 1, e, mapping, Q_edge, Phis_edge)
+                            nodes.extend(IntegralMomentOfDerivative(ref_el, s, Q, phi, comp=(k,), shp=shp)
+                                         for phi in phis for k in range(sd))
+
+                    # Rest of the facet moments
                     Q, phis = map_duals(ref_el, dim, entity, mapping, Q_ref, Phis)
-                    nodes.extend(FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis)
+                    if dim == sd:
+                        nodes.extend(FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis)
+                    else:
+                        nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(k,), shp=shp)
+                                     for phi in phis for k in range(sd))
                     entity_ids[dim][entity] = list(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
-    def _reference_duals(self, dim, degree):
+    def _reference_duals(self, dim, degree, moment_degree):
         facet = symmetric_simplex(dim)
         Q = create_quadrature(facet, 2 * degree)
         Qpts, Qwts = Q.get_points(), Q.get_weights()
 
-        P = ONPolynomialSet(facet, degree, self.shp, scale="orthonormal")
-        P_at_qpts = P.tabulate(Qpts, 1)
-        trial = P_at_qpts[(0,) * dim]
-
-        if dim == self.shp[0]:
-            dtrial = eps(P_at_qpts)
-            moments = self._interior_moments
-        else:
-            dtrial = grad(P_at_qpts)
-            moments = self._facet_moments
+        shp = (dim,) if dim == self.sd else tuple()
+        V = ONPolynomialSet(facet, degree, shp, scale="orthonormal")
 
-        K = moments(facet, degree, Qpts, Qwts, dtrial)
-        duals = numpy.tensordot(K, P_at_qpts[(0,) * dim], axes=(1, 0))
+        moments = self._interior_moments if dim == self.sd else self._facet_moments
+        duals = moments(facet, moment_degree, Qpts, Qwts, V)
         return Q, duals
 
-    def _facet_moments(self, facet, degree, Qpts, Qwts, trial):
+    def _facet_moments(self, facet, moment_degree, Qpts, Qwts, V):
         """Integrate trial expressions against an orthonormal basis for
            the exterior derivative of bubbles.
         """
         dim = facet.get_spatial_dimension()
-        # Get bubbles
-        B = make_bubbles(facet, degree, shape=self.shp)
 
-        # Tabulate the exterior derivate
-        B_at_qpts = B.tabulate(Qpts, 1)
-        dtest = grad(B_at_qpts)
-        test = B_at_qpts[(0,) * dim]
-        expr = dtest
-        if len(dtest) > 0:
-            # Build an orthonormal basis, remove nullspace
-            B = inner(test, test, Qwts)
-            A = inner(dtest, dtest, Qwts)
-            sig, S = scipy.linalg.eigh(A, B)
-            tol = sig[-1] * 1E-12
-            nullspace_dim = len([s for s in sig if abs(s) <= tol])
-            S = S[:, nullspace_dim:]
-            S *= numpy.sqrt(1 / sig[None, nullspace_dim:])
-            # Apply change of basis
-            expr = numpy.tensordot(S, expr, axes=(0, 0))
-        return inner(expr, trial, Qwts)
-
-
-    def _interior_moments(self, facet, degree, Qpts, Qwts, eps_trial):
+        V_at_qpts = V.tabulate(Qpts)
+        trial = V_at_qpts[(0,) * dim]
+        test = trial[:polynomial_dimension(facet, moment_degree)]
+        K = inner(test, trial, Qwts)
+        return numpy.tensordot(K, trial, axes=(1, 0))
+
+    def _interior_moments(self, facet, moment_degree, Qpts, Qwts, V):
         """Integrate trial expressions against an orthonormal basis for
            the exterior derivative of bubbles.
         """
         dim = facet.get_spatial_dimension()
+
+        V_at_qpts = V.tabulate(Qpts, 1)
+        trial = V_at_qpts[(0,) * dim]
+        eps_trial = eps(V_at_qpts)
+
         # Get bubbles
-        B = make_bubbles(facet, degree, shape=self.shp)
+        B = make_bubbles(facet, V.degree, shape=(dim,))
 
         # Tabulate the exterior derivate
         B_at_qpts = B.tabulate(Qpts, 1)
@@ -158,8 +185,8 @@ def _interior_moments(self, facet, degree, Qpts, Qwts, eps_trial):
         # Apply change of basis
         expr = [inner(numpy.tensordot(S1, div_test, axes=(0, 0)), div_trial, Qwts),
                 inner(numpy.tensordot(S2, eps_test, axes=(0, 0)), eps_trial, Qwts)]
-        expr = numpy.concatenate(expr, axis=0)
-        return expr
+        K = numpy.concatenate(expr, axis=0)
+        return numpy.tensordot(K, trial, axes=(1, 0))
 
 
 class Stokes(finite_element.CiarletElement):
@@ -170,50 +197,43 @@ def __init__(self, ref_el, degree):
             raise ValueError(f"{type(self).__name__} elements only valid for k >= {2*sd}")
         poly_set = ONPolynomialSet(ref_el, degree, shape=(sd,), variant="bubble")
         dual = StokesDual(ref_el, degree)
-
-        formdegree = 0  # 0-form
-        super().__init__(poly_set, dual, degree, formdegree)
+        formdegree = sd-1  # (n-1)-form
+        super().__init__(poly_set, dual, degree, formdegree,
+                         mapping="contravariant piola")
 
 
 if __name__ == "__main__":
     import matplotlib.pyplot as plt
-    from FIAT import ufc_simplex
 
     dim = 2
     degree = 10
     ref_el = symmetric_simplex(dim)
-    # ref_el = ufc_simplex(dim)
     Q = create_quadrature(ref_el, 2 * degree)
     Qpts, Qwts = Q.get_points(), Q.get_weights()
 
-    space_dict = {"H1": (Stokes, eps)}
-    spaces = list(space_dict.keys())
-
-    fig, axes = plt.subplots(ncols=2, nrows=len(spaces), figsize=(12, 6*len(spaces)))
+    fig, axes = plt.subplots(ncols=2, nrows=1, figsize=(12, 6))
     axes = axes.flat
 
-    for space in spaces:
-        element, d = space_dict[space]
-        fe = element(ref_el, degree)
-        phi_at_qpts = fe.tabulate(1, Qpts)
-
-        Veps = d(phi_at_qpts)
-        Vdiv = div(phi_at_qpts)
-        Aeps = inner(Veps, Veps, Qwts)
-        Adiv = inner(Vdiv, Vdiv, Qwts)
-
-        title = f"{type(fe).__name__}({degree})"
-        names = ("eps", "div")
-        mats = (Aeps, Adiv)
-        for name, A in zip(names, mats):
-            A[abs(A) < 1E-10] = 0.0
-            scipy_mat = scipy.sparse.csr_matrix(A)
-            nnz = scipy_mat.count_nonzero()
-            ms = 0
-            ax = next(axes)
-            ax.spy(A, markersize=ms)
-            if ms == 0:
-                ax.pcolor(numpy.log(abs(A)))
-            ax.set_title(f"{title} {name} nnz {nnz}")
+    fe = Stokes(ref_el, degree)
+    phi_at_qpts = fe.tabulate(1, Qpts)
+
+    Veps = eps(phi_at_qpts)
+    Vdiv = numpy.trace(Veps, axis1=1, axis2=2)
+    Aeps = inner(Veps, Veps, Qwts)
+    Adiv = inner(Vdiv, Vdiv, Qwts)
+
+    title = f"{type(fe).__name__}({degree})"
+    names = ("eps", "div")
+    mats = (Aeps, Adiv)
+    for name, A in zip(names, mats):
+        A[abs(A) < 1E-10] = 0.0
+        scipy_mat = scipy.sparse.csr_matrix(A)
+        nnz = scipy_mat.count_nonzero()
+        ms = 0
+        ax = next(axes)
+        ax.spy(A, markersize=ms)
+        if ms == 0:
+            ax.pcolor(numpy.log(abs(A)))
+        ax.set_title(f"{title} {name} nnz {nnz}")
 
     plt.show()