Compute gradient of scalar field on BoundaryMesh

variational formulation
2

For domains with corners, I would recommend extracting the flux from each smooth section of the boundary separately. When extracting flux from only a portion of the boundary, there is an extra “trick” that needs to be added on top of what I wrote up in the earlier post that you linked to.

Take a look at the following example, which I prepared some time ago for a research collaborator of mine; it allows for 3/2 rate of convergence in L^2 on the boundary, which is, I think, about as good as you can hope to get with linear elements:

from dolfin import *
Nel = 100
mesh = UnitSquareMesh(Nel,Nel)
n = FacetNormal(mesh)
V = FunctionSpace(mesh,"Lagrange",1)

# The full boundary, on which we apply a Dirichlet BC:
class Bdry(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary

# Interior of the domain:
class Interior(SubDomain):
    def inside(self,x,on_boundary):
        return (not (near(x[0],0) or near(x[0],1)
                     or near(x[1],0) or near(x[1],1)))
    
# The boundary on the right side of the domain, on which we want to
# extract the flux:
class BdryOfInterest(SubDomain):
    def inside(self,x,on_boundary):
        return (x[0] > 1.0 - DOLFIN_EPS) and on_boundary
    
# Everything BUT the right side of the domain from which we want to
# extract the boundary flux; this must include all nodes not used to
# approximate flux on the boundary of interest, including those in the
# interior of the domain (so `on_boundary` should not be used in the
# return value).
class AntiBdry(SubDomain):
    def inside(self,x,on_boundary):
        return x[0] < 1.0 - DOLFIN_EPS

# Mark parts of the boundary to which integrals will need to be restricted.
# (Parts not explicitly marked are flagged with zero.)
FLUX_BDRY = 1
COMPLEMENT_FLUX_BDRY = 2
boundaryMarkers = MeshFunction("size_t",mesh,mesh.topology().dim()-1,
                               COMPLEMENT_FLUX_BDRY)
BdryOfInterest().mark(boundaryMarkers,FLUX_BDRY)
ds_marked = ds(subdomain_data=boundaryMarkers)

# Desired exact solution:
u_exact = Expression("x[0]*x[0] + x[1]*x[1]", degree=2, domain=mesh)

# Boundary conditions to apply to solution, `u_h`:
BCs = [DirichletBC(V,u_exact,Bdry()),]

# Boundary conditions to apply when solving for flux on whole boundary:
interiorBCs = [DirichletBC(V,Constant(0.0),Interior()),]

# Boundary conditions to apply to the flux solution when we are only
# interested in flux on the right side of the domain:
antiBCs = [DirichletBC(V,Constant(0.0),AntiBdry()),]

# Residual of the Poisson equation
f = -div(grad(u_exact))
def R(u,v):
    return (inner(grad(u),grad(v)) - f*v)*dx

u = TrialFunction(V)
v = TestFunction(V)
F = R(u,v)
a = lhs(F)
L = rhs(F)

# Solve and put solution in `u_h`
u_h = Function(V)
solve(a==L,u_h,BCs)

# Solve in the trace space for what Neumann data, i.e.,
#
#  \mathbf{q}\cdot\mathbf{n} = \nabla u\cdot\mathbf{n}
#
# would have produced the solution from the Dirichlet problem; this is the
# notion of flux that satisfies the underlying conservation law.
qn = TrialFunction(V)

################################################################################

# The trick:  Since we want to use the corner nodes to approximate the
# flux on our boundary of interest, test functions will end up being
# nonzero on an $O(h)$ part of the complement of the boundary of interest.
# Thus we need to integrate a consistency term on that part of the boundary.

n = FacetNormal(mesh)
consistencyTerm = inner(grad(u_h),n)*v*ds_marked(COMPLEMENT_FLUX_BDRY)
FBdry = qn*v*ds_marked(FLUX_BDRY) - R(u_h,v) + consistencyTerm

FBdry_inconsistent = qn*v*ds_marked(FLUX_BDRY) - R(u_h,v)

# Applying the consistent flux extraction on the full boundary
# and then restricting the result is also sub-optimal; this flux extraction
# technique doesn't appear to play nicely with corners.
FBdry_full = qn*v*ds - R(u_h,v)

################################################################################

# Get $\mathbf{q}\cdot\mathbf{n}$ on the boundary of interest with and
# without the consistency term:
def solveFor_qn_h(FBdry, BCs):
    aBdry = lhs(FBdry)
    LBdry = rhs(FBdry)
    ABdry = assemble(aBdry,keep_diagonal=True)
    bBdry = assemble(LBdry)
    [BC.apply(ABdry,bBdry) for BC in BCs]
    qn_h = Function(V)
    solve(ABdry,qn_h.vector(),bBdry)
    return qn_h

qn_h = solveFor_qn_h(FBdry, antiBCs)
qn_h_inconsistent = solveFor_qn_h(FBdry_inconsistent, antiBCs)
qn_h_full = solveFor_qn_h(FBdry_full, interiorBCs)

# Compare fluxes with the exact solution:
import math
def fluxErr(qn):
    err = inner(grad(u_exact),n) - qn
    return math.sqrt(assemble(err*err*ds_marked(FLUX_BDRY)))

# Converges at 3/2 order in $L^2$, as expected for smooth solutions, according
# to this paper:
#
#  https://link.springer.com/article/10.1007/BF01385871
#
print("Error in consistent flux on boundary of interest: "
      +str(fluxErr(qn_h)))

# Converges at only 1/2 order:
print("Error in restricted consistent flux from full boundary: "
      +str(fluxErr(qn_h_full)))

# Converges at only 1/2 order:
print("Error in inconsistent flux on boundary of interst: "
      +str(fluxErr(qn_h_inconsistent)))

# Converges at first order:
print("Error in direct flux: "+str(fluxErr(inner(grad(u_h),n))))

# Slow for large meshes; also first-order:
#print("Error in projection to linears: "
#      +str(fluxErr(inner(project(grad(u_h),
#                                 VectorFunctionSpace(mesh,"CG",1)),n))))
3 Likes