Hi all,

I have a surface problem that takes into account the gradient of a solution from a volume problem. I have then 2 meshes, one being the BoundaryMesh of the other. I’d like to compute the gradient of solution u and then use it in the surface formulation F .

from fenics import *

mesh = UnitCubeMesh(15, 15, 15)
boundarymesh = BoundaryMesh(mesh, 'exterior')
n = FacetNormal(mesh)

V = FunctionSpace(mesh, "CG", 1)
u = Function(V)

W = FunctionSpace(boundarymesh, "CG", 1)
w = TrialFunction(W)
w_test = TestFunction(W)

F = w*w_test*dx + dot(gradient, n)*w_test*dx #dummy formulation that requires the normal gradient


Rem

At the risk of reading too much into the question, I will first assume that you specifically need \nabla u\cdot\mathbf{n} on the boundary, rather than the full vector \nabla u (which avoids the extra complication of how to get the normal \mathbf{n} on boundarymesh). If you have \nabla\cdot\mathbf{n} as a nodal Function on mesh, then you can interpolate it in the space W on boundarymesh relatively easily. The question is then “What is the best way to obtain a nodal function that represents \nabla u\cdot\mathbf{n}?” If I assume further that \nabla u\cdot\mathbf{n} is the flux of a PDE that was solved on mesh, then the following is a good way to extract a nodal flux that is consistent with the interior solution (cf. Exercise 8 in Section 2.12 of Tom Hughes’s FEM book):

from fenics import *

mesh = UnitCubeMesh(15, 15, 15)
boundarymesh = BoundaryMesh(mesh, 'exterior')
n = FacetNormal(mesh)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
W = FunctionSpace(boundarymesh, "CG", 1)

# Replace with the actual residual of the PDE used to compute u; in this case
# the flux for this problem will be the normal gradient, as in the MWE.
x = SpatialCoordinate(mesh)
f = sin(pi*x)
def PDE_res(u,v):

# Dirichlet BC on boundary where you want to extract flux:
bcStr = ""
for i in range(0,3):
bcStr += "near(x["+str(i)+"],0.0)||near(x["+str(i)+"],1.0)"
if(i<2):
bcStr += "||"
bcStr = "("+bcStr+")"
uBC = DirichletBC(V,Constant(0.0),bcStr)

# Calculate something for u:
v = TestFunction(V)
res = PDE_res(TrialFunction(V),v)
solve(lhs(res)==rhs(res),u,uBC)

# Calculate a nodal flux consistent with the conservation law, i.e., find a
# nodal field h such that, if it were applied as Neumann data, it would have
# produced the solution u, which is considered given in this problem.
resFlux = PDE_res(u,v) - TrialFunction(V)*v*ds
h = Function(V)
# (Need to pass special keep_diagonal argument to LHS assembly, because
# the TrialFunction is only assembled on the boundary, and there will
# otherwise be a zero block in the sparse structure when we go to apply BCs
# on the interior, leading to errors.)
ABdry = assemble(lhs(resFlux),keep_diagonal=True)
bBdry = assemble(rhs(resFlux))
# The flux is defined only on boundary nodes, so set everything else to zero
# to obtain a well-posed linear algebra problem.
hBC = DirichletBC(V,Constant(0.0),"!"+bcStr)
hBC.apply(ABdry,bBdry)
solve(ABdry,h.vector(),bBdry)

# Interpolate the flux h defined on boundary nodes of mesh onto the space
# W defined on boundarymesh.  (This step is exact, because W is the
# trace space of V.)
h.set_allow_extrapolation(True)
hbdry = interpolate(h, W)

# Test:  Net flux should balance net source term from corresponding PDE
# (a.k.a. the "compatibility condition"):
print("Flux: "+str(assemble(hbdry*dx)))
print("Source term: "+str(assemble(f*dx(domain=mesh))))

# (This also verifies that flux hbdry interpolated on boundary mesh can be
# used in variational forms.)

2 Likes