I had exactly the same problem before. Here are some my solutions.
In 2D, UFL has enough expressiveness to achieve it. The following is due to Douglas Arnold:
# Canonical interpolant onto N1e1
def Pih(u0, mesh):
V = FunctionSpace(mesh, 'N1curl', r)
u = TrialFunction(V)
if r == 1:
W = FunctionSpace(mesh, 'N1curl', 1)
v0 = TestFunction(W)
else:
V0 = FunctionSpace(mesh, 'N1curl', r, restriction='facet')
V1 = VectorFunctionSpace(mesh, 'DG', r-2)
W = MixedFunctionSpace([V0, V1])
v0, v1 = TestFunctions(W)
# Define the tangent to element facets (edges)
n = FacetNormal(mesh)
t = as_vector((-n[1], n[0]))
# Define the bilinear form
#
# A(u, (v0, v1)) = int_Gamma u.t v0.t ds + int u.v1 dx
#
# where the first integral is over all the edges of the mesh, and the second
# over the domain. v0 is in the N2curl space with all non-edge DOFs set to zero
# so v0 -> v0.t is an isomorphism onto DG_{r-1} on the edges. v1 is a vector DG
A = (inner(u, t) * inner(v0, t))('+') * dS + inner(u, t) * inner(v0, t) * ds
# Define the right hand side, L(v0, v1) = A(u0, (v0, v1))
L = (inner(u0, t) * inner(v0, t))('+') * dS + inner(u0, t) * inner(v0, t) * ds
if r > 1:
A = A + inner(u, v1) * dx
L = L + inner(u0, v1) * dx
# Solve the resulting linear system for u = Pih(u0)
u = Function(V)
solve(A == L, u)
return u
In 3D, the above does not work because UFL does not have the expression for edge integrals in 3D. I wrote the following code computating the canonical interpolant by hand (does not work in parallel, only for the lowest order):
from dolfin import *
from FIAT import *
import numpy
def tangential_edge_integral(mesh, f, quad_degree):
"""
This function uses quadrature to compute the average of
of the dot product f.t on each edge of the given mesh.
Here f is a given vector valued function and t is the vector
pointing from the first to the second vertex of the edge.
The returned value is a numpy array edge_int where
edge_int[edge.index()] = avg of f.t on that edge
"""
# Reference interval
L = reference_element.DefaultLine()
# Coordinates for the reference interval
x = numpy.array(L.get_vertices())
# Make quadrature rule on the line
Q = make_quadrature(L, quad_degree)
# Quadrature points
qp = Q.get_points()
# Quadrature weights
qw = Q.get_weights()
# Get mesh vertex coordinates
coor = mesh.coordinates()
# Create an array to save the result
edge_int = numpy.zeros(mesh.num_edges())
# Integrate over all the edges
for e in edges(mesh):
# Extract end points of the edge
end_pts = e.entities(0)
y = numpy.vstack([coor[end_pts[0]], coor[end_pts[1]]])
# Compute the tangent. Following the FEniCS convention this
# is the vector pointing from the first to the end point
# of the edge.
t = y[1] - y[0]
# Make the affine map from the reference line
(A, b) = reference_element.make_affine_mapping(x, y)
# Quadrature points on the physical edge
phys_qp = map(lambda x: A.dot(x)+b, qp)
# Evaluate f.t at quadrature points
vals = map(lambda x: t.dot(f(x)), phys_qp)
# Approximate average of f.t on the edge using the quadrature rule
edge_int[e.index()] = qw.dot(vals)/2.0
return edge_int
def build_edge2dof_map(V):
"""
This function takes a N1Curl(1) space and return an integer valued array edge2dof.
This array has the number of edges as its length. In particular
edge2dof[i] = j
means that dof #i, that is u.vector()[i], is associated to edge #j.
"""
# Extract the cell to edge map (given an cell index, it returns the indices of its edges)
cell2edges = V.mesh().topology()(3, 1)
# Extract the cell dofmap (given a cell index, it returns the dof numbers)
cell2dofs = V.dofmap().cell_dofs
# Array to save the result
edge2dof = numpy.zeros(mesh.num_edges(), dtype="int")
# Iterate over cells, associating the edges to the dofs for that cell
for c in range(mesh.num_cells()):
# get the global edge numbers for this cell
c_edges = cell2edges(c)
# get the global dof numbers for this cell
c_dofs = cell2dofs(c)
# associate the edge numbers to the corresponding dof numbers
edge2dof[c_dofs] = c_edges
# This algorithm might not look fast as it does quite some redundant work. In actual
# runs, for most meshes, this is not the most time consuming step and does not take
# more than a milisecond.
return edge2dof
def n1curl_1_canonical_projection(f, mesh, quad_degree):
# Initialize the mesh
mesh.init()
# Compute the average edge integrals
edge_int = tangential_edge_integral(mesh, f, quad_degree)
# Create the return value
V = FunctionSpace(mesh, "N1curl", 1)
u = Function(V)
# Permute the average edge integrals to match the order of the dofs.
u.vector()[:] = edge_int.take(build_edge2dof_map(V))
return u
f = Expression(("1.0 - x[1] + 2*x[2]", "3.0 + x[0] + 3*x[2]", "2.0 - 2*x[0] - 3*x[1]"))
domain = Box(0., 0., 0., 1., 1., 1.)
n = 8; mesh = Mesh(domain, n)
quad_degree = 2
V = FunctionSpace(mesh, "N1curl", 1)
u = n1curl_1_canonical_projection(f, mesh, quad_degree)
v = interpolate(f, V)
print(assemble(dot(u-v,u-v)*dx))
It should be easy to use the above code to assemble the matrix you want by hand (not efficient though).
