H(curl) constraints for BDDC

Author: pbrubeck

firedrake/preconditioners/bddc.py

@@ -104,7 +104,8 @@ def initialize(self, pc):
                 div_args = (divergence,)
                 div_kwargs = dict()
             bddcpc.setBDDCDivergenceMat(*div_args, **div_kwargs)
-        elif tdim >= 2 and V.finat_element.formdegree == 1:
+
+        elif tdim >= 3 and V.finat_element.formdegree == 1:
             get_gradient = appctx.get("get_discrete_gradient", get_discrete_gradient)
             gradient = get_gradient(V)
             try:
@@ -130,12 +131,19 @@ def applyTranspose(self, pc, x, y):
         self.pc.applyTranspose(x, y)
 
 
-def get_vertex_dofs(V):
-    W = FunctionSpace(V.mesh(), restrict(V.ufl_element(), "vertex"))
+def get_dof_mapping(V, W):
     indices = get_restriction_indices(V, W)
     indices = V.dof_dset.lgmap.apply(indices)
-    vertex_dofs = PETSc.IS().createGeneral(indices, comm=V.comm)
-    return vertex_dofs
+    return PETSc.IS().createGeneral(indices, comm=V.comm)
+
+
+def get_restricted_dofs(V, domain):
+    W = FunctionSpace(V.mesh(), restrict(V.ufl_element(), domain))
+    return get_dof_mapping(V, W)
+
+
+def get_low_order_dofs(V):
+    return get_dof_mapping(V, V.reconstruct(degree=1))
 
 
 def get_divergence_mat(V, mat_type="aij", allow_repeated=False):
@@ -153,13 +161,41 @@ def get_divergence_mat(V, mat_type="aij", allow_repeated=False):
 
 
 def get_discrete_gradient(V):
-    degree = max(as_tuple(V.ufl_element().degree()))
-    variant = V.ufl_element().variant()
+    from firedrake import VectorSpaceBasis
+
     Q = FunctionSpace(V.mesh(), curl_to_grad(V.ufl_element()))
+    degree = max(as_tuple(Q.ufl_element().degree()))
+    variant = Q.ufl_element().variant()
     gradient = tabulate_exterior_derivative(Q, V)
-    if variant == 'fdm':
-        corners = get_vertex_dofs(Q)
-        gradient.compose('_elements_corners', corners)
+
+    if variant == "fdm":
+        vdofs = get_restricted_dofs(Q, "vertex")
+        gradient.compose('_elements_corners', vdofs)
+
+        # Constant moments along edges
+        wdofs = get_low_order_dofs(V)
+        v0 = Function(V)
+        with v0.dat.vec as y:
+            y.setValues(wdofs, numpy.ones((wdofs.getSizes()[0],), dtype=PETSc.RealType))
+
+        # Linear moments along edges
+        # Constant moments of gradients of H1 edge bubbles
+        p = Function(Q).interpolate(1)
+        q = Function(Q)
+
+        edofs = get_restricted_dofs(Q, "edge")
+        edofs = edofs.difference(vdofs)
+        with q.dat.vec as qvec, p.dat.vec as pvec:
+            qvec.setValues(edofs, pvec.getValues(edofs))
+
+        v1 = Function(V)
+        with v1.dat.vec as y, q.dat.vec as x:
+            gradient.mult(x, y)
+            y.chop(1E-14)
+
+        nsp = VectorSpaceBasis([v0, v1])
+        gradient.compose("_constraints", nsp.nullspace())
+
     grad_args = (gradient,)
     grad_kwargs = {'order': degree}
     return grad_args, grad_kwargs