Project gradient on BoundaryMesh

variational formulation
1

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)


gradient = project(grad(u), VectorFunctionSpace(boundarymesh, "CG", 1)) #doesn't work

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

Thank you in advance for your help !

Rem

2

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[0])
def PDE_res(u,v):
    return inner(grad(u),grad(v))*dx - f*v*dx

# 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