Thank you, that was very useful! I followed your suggestion to use a mixed form, and it works like a charm.
Here is what I have done: I want to solve \nabla^2 \nabla^2 u = f in \Omega = [0,1] \times [0,1] with BCs u = g and \nabla^2 u = h on \partial \Omega.
I introduce the auxiliary fields v, w defined by \nabla^2 u = v, \nabla^2 v = w. Notice that I introduced w only to extract and plot \nabla^2 \nabla^2 u from the solution output, I could not just take the Laplacian of v because this would have introduced blow-ups according to your comment.
I obtain the mixed problem for u, v, w
\nabla^2 v = f,
\nabla^2 u = v,
\nabla^2 v = w
with BCs u = g, v = h on \partial \Omega. I solve this with the following code by building an example where I know the solution, u = \cos(x+y) \sin(x-y), f = 8 [\sin(2x) - \sin(2y)], and similarly for g, h.
from fenics import *
import argparse
from mshr import *
import ufl as ufl
from dolfin import *
L = 1.0
h = 1.0
function_space_degree = 4
i, j, k, l = ufl.indices( 4 )
parser = argparse.ArgumentParser()
parser.add_argument( "input_directory" )
parser.add_argument( "output_directory" )
args = parser.parse_args()
xdmffile_u = XDMFFile( (args.output_directory) + "/u.xdmf" )
xdmffile_v = XDMFFile( (args.output_directory) + "/v.xdmf" )
xdmffile_w = XDMFFile( (args.output_directory) + "/w.xdmf" )
xdmffile_check = XDMFFile( (args.output_directory) + "/check.xdmf" )
xdmffile_check.parameters.update( {"functions_share_mesh": True, "rewrite_function_mesh": False} )
# read the mesh
mesh = Mesh()
xdmf = XDMFFile( mesh.mpi_comm(), (args.input_directory) + "/triangle_mesh.xdmf" )
xdmf.read( mesh )
# radius of the smallest cell in the mesh
r_mesh = mesh.hmin()
print( f"Mesh radius = {r_mesh}" )
# read the triangles
mvc = MeshValueCollection( "size_t", mesh, mesh.topology().dim() )
with XDMFFile( (args.input_directory) + "/triangle_mesh.xdmf" ) as infile:
infile.read( mvc, "name_to_read" )
cf = cpp.mesh.MeshFunctionSizet( mesh, mvc )
xdmf.close()
# read the lines
mvc = MeshValueCollection( "size_t", mesh, mesh.topology().dim() - 1 )
with XDMFFile( (args.input_directory) + "/line_mesh.xdmf" ) as infile:
infile.read( mvc, "name_to_read" )
sf = cpp.mesh.MeshFunctionSizet( mesh, mvc )
xdmf.close()
boundary = 'on_boundary'
boundary_lr = 'near(x[0], 0) || near(x[0], 1.0)'
boundary_tb = 'near(x[1], 0) || near(x[1], 1.0)'
dx = Measure( "dx", domain=mesh, subdomain_data=cf, subdomain_id=1 )
ds_l = Measure( "ds", domain=mesh, subdomain_data=sf, subdomain_id=2 )
ds_r = Measure( "ds", domain=mesh, subdomain_data=sf, subdomain_id=3 )
ds_t = Measure( "ds", domain=mesh, subdomain_data=sf, subdomain_id=4 )
ds_b = Measure( "ds", domain=mesh, subdomain_data=sf, subdomain_id=5 )
ds = ds_l + ds_r + ds_t + ds_b
n = FacetNormal( mesh )
P_u = FiniteElement( 'P', triangle, function_space_degree )
P_v = FiniteElement( 'P', triangle, function_space_degree )
P_w = FiniteElement( 'P', triangle, function_space_degree )
element = MixedElement( [P_u, P_v, P_w] )
Q = FunctionSpace( mesh, element )
Q_u = Q.sub( 0 ).collapse()
Q_v = Q.sub( 1 ).collapse()
Q_w = Q.sub( 2 ).collapse()
Q_grad_v = VectorFunctionSpace( mesh, 'P', function_space_degree )
class u_exact_expression( UserExpression ):
def eval(self, values, x):
values[0] = cos(x[0]+x[1]) * sin(x[0]-x[1])
def value_shape(self):
return (1,)
class v_exact_expression( UserExpression ):
def eval(self, values, x):
values[0] = - 4 * cos(x[0])*sin(x[0]) + 4 * cos(x[1])*sin(x[1])
def value_shape(self):
return (1,)
class w_exact_expression( UserExpression ):
def eval(self, values, x):
values[0] = 8 * (sin(2*x[0]) - sin(2*x[1]))
def value_shape(self):
return (1,)
psi = Function( Q )
nu_u, nu_v, nu_w = TestFunctions( Q )
grad_v = Function( Q_grad_v )
u_output = Function( Q_u )
v_output = Function( Q_v )
w_output = Function( Q_w )
u_exact = Function( Q_u )
v_exact = Function( Q_v )
w_exact = Function( Q_w )
f = Function( Q_w )
J_uvw = TrialFunction( Q )
u, v, w = split( psi )
u_exact.interpolate( u_exact_expression( element=Q_u.ufl_element() ) )
v_exact.interpolate( v_exact_expression( element=Q_v.ufl_element() ) )
w_exact.interpolate( w_exact_expression( element=Q_w.ufl_element() ) )
f.interpolate( w_exact_expression( element=Q_w.ufl_element() ) )
u_profile = Expression( 'cos(x[0]+x[1]) * sin(x[0]-x[1])', L=L, h=h, element=Q.sub( 0 ).ufl_element() )
v_profile = Expression( '- 4 * cos(x[0])*sin(x[0]) + 4 * cos(x[1])*sin(x[1])', L=L, h=h, element=Q.sub( 1 ).ufl_element() )
bc_u = DirichletBC( Q.sub( 0 ), u_profile, boundary )
bc_v = DirichletBC( Q.sub( 1 ), v_profile, boundary )
F_v = ((v.dx( i )) * (nu_u.dx( i )) + f * nu_u) * dx \
- n[i] * (v.dx( i )) * nu_u * ds
F_u = ((u.dx( i )) * (nu_v.dx( i )) + v * nu_v) * dx \
- n[i] * (u.dx( i )) * nu_v * ds
F_w = ((v.dx( i )) * (nu_w.dx( i )) + w * nu_w) * dx \
- n[i] * (v.dx( i )) * nu_w * ds
F = F_u + F_v + F_w
bcs = [bc_u, bc_v]
J = derivative( F, psi, J_uvw )
problem = NonlinearVariationalProblem( F, psi, bcs, J )
solver = NonlinearVariationalSolver( problem )
solver.solve()
u_output, v_output, w_output = psi.split( deepcopy=True )
xdmffile_u.write( u_output, 0 )
xdmffile_v.write( v_output, 0 )
xdmffile_w.write( w_output, 0 )
print( "BCs check: " )
print( f"\t<<(u-u_exact)^2>>_\partial Omega = {assemble( (u_output - u_exact) ** 2 * ds ) / assemble( Constant( 1.0 ) * ds )}" )
print( f"\t<<(v-v_exact)^2>>_\partial Omega = {assemble( (v_output - v_exact) ** 2 * ds ) / assemble( Constant( 1.0 ) * ds )}" )
print( f"\t<<(w-w_exact)^2>>_\partial Omega = {assemble( (w_output - w_exact) ** 2 * ds ) / assemble( Constant( 1.0 ) * ds )}" )
print( "Comparison with exact solution: " )
print( f"\t<<(u - u_exact)^2>> = {sqrt( assemble( ((u_output - u_exact) ** 2) * dx ) / assemble( Constant( 1.0 ) * dx ) )}" )
print( f"\t<<(v - v_exact)^2>> = {sqrt( assemble( ((v_output - v_exact) ** 2) * dx ) / assemble( Constant( 1.0 ) * dx ) )}" )
print("Check that the PDE is satisfied: ")
print( f"\t<<(w-f)^2>>_\partial Omega = {assemble( (w_output - f) ** 2 * ds ) / assemble( Constant( 1.0 ) * ds )}" )
which runs ok
Mesh radius = 0.017499012335440677
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.454e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.921e-12 (tol = 1.000e-10) r (rel) = 1.321e-14 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
BCs check:
<<(u-u_exact)^2>>_\partial Omega = 1.4444230296665588e-34
<<(v-v_exact)^2>>_\partial Omega = 2.253436960751194e-33
<<(w-w_exact)^2>>_\partial Omega = 1.1049231523976149e-11
Comparison with exact solution:
<<(u - u_exact)^2>> = 4.221561446442978e-13
<<(v - v_exact)^2>> = 1.6254659235591433e-12
Check that the PDE is satisfied:
<<(w-f)^2>>_\partial Omega = 1.1049231523976149e-11
The last line shows that \nabla^2 \nabla^2 u is very close to f.
Is this the solution that you meant ?