Hi all,
Given a FE-solution u, I want to solve another PDE by using u in the boundary condition for this PDE:
I can't handle the boundary condition, as I have to compute 
and 
.
The first problem is the computation of :
Although there is no error message, I'm not sure, if using div(grad(u)) in the boundary condition is right, as I'm using FunctionSpace(mesh, 'Lagrange', 1) ?
The second problem is the computation of :
Using dot(nabla_grad(u),n) as boundary condition provides an error: ufl.log.UFLException: Integral of type cell cannot contain a FacetNormal.
I found a solution by manually calculating the FacetNormal for 3D-meshes here, however, my problem is defined on a 2D-mesh. Is there a simple implementation in Fenics for this, or have I to compute the FacetNormal manually (maybe similar to the solution in the link)?
Here is the code providing the error message:
from dolfin import *
mesh = UnitSquareMesh(2,2)
V = FunctionSpace(mesh, 'Lagrange', 1)
# Let u be a FE-solution
def u0_boundary(x, on_boundary):
return on_boundary
bc_u = DirichletBC(V, Constant(3.0) , u0_boundary)
u = TrialFunction(V)
v = TestFunction(V)
a = -inner(nabla_grad(u), nabla_grad(v))*dx
L = Constant(-6.0)*v*dx
u = Function(V)
solve(a == L, u, bc_u)
# define boundary conditions for new problem
n=FacetNormal(mesh)
bc=div(grad(u)) - (dot(nabla_grad(u),n)-Constant(1))
bc_w= DirichletBC(V, bc , u0_boundary)
# solve another PDE with special boundary condition
w = TrialFunction(V)
v1 = TestFunction(V)
a1 = inner(nabla_grad(w), nabla_grad(v1))*dx
L1 = Constant(-6.0)*v1*dx
w = Function(V)
solve(a1 == L1, w, bc_w)
