Hello,
Following up on this previous post, I am trying to solve a fourth-order PDE. I constructed this minimal example in a unit square \Omega for you:
\nabla^2 [\partial_i (z \partial_i z)] = f
with BCs z = z_0 and n^i \partial_i z = 0 on \partial \Omega. I considered a case where I know the exact solution, z = (x^4+y^4)/48, and constructed f and the BCs accordingly.
How can I compute \nabla^2 [\partial_i (z \partial_i z)] and verify that it is equal to f for all mesh points in \Omega?
In the minimal working example below, I solve the problem by introducing the auxiliary fields \omega_i \equiv \partial_i z , \mu = \partial_i (z \omega_i), \rho_i \equiv \partial_i \mu and \tau \equiv \partial_i \mu_i. I should have \tau = f in \Omega, but there are some odd deviations at the boundary, here is \tau - f:
Here is the code:
from fenics import *
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 )
xdmffile_check = XDMFFile( "solution/check.xdmf" )
xdmffile_check.parameters.update( {"functions_share_mesh": True, "rewrite_function_mesh": False} )
mesh = Mesh()
mesh = RectangleMesh( Point( 0, 0 ), Point( L, h ), 20, 20 )
r_mesh = mesh.hmin()
mvc = MeshValueCollection( "size_t", mesh, mesh.topology().dim() )
cf = dolfin.cpp.mesh.MeshFunctionSizet( mesh, mvc )
mvc = MeshValueCollection( "size_t", mesh, mesh.topology().dim() - 1 )
sf = dolfin.cpp.mesh.MeshFunctionSizet( mesh, mvc )
boundary = 'on_boundary'
dx = Measure( "dx", domain=mesh, subdomain_data=cf )
ds = Measure( "ds", domain=mesh, subdomain_data=sf )
n = FacetNormal( mesh )
P_z = FiniteElement( 'P', triangle, function_space_degree )
P_omega = VectorElement( 'P', triangle, function_space_degree )
P_mu = FiniteElement( 'P', triangle, function_space_degree )
P_rho = VectorElement( 'P', triangle, function_space_degree )
P_tau = FiniteElement( 'P', triangle, function_space_degree )
element = MixedElement( [P_z, P_omega, P_mu, P_rho, P_tau] )
Q = FunctionSpace( mesh, element )
Q_z = Q.sub( 0 ).collapse()
Q_omega = Q.sub( 1 ).collapse()
Q_mu = Q.sub( 2 ).collapse()
Q_rho = Q.sub( 3 ).collapse()
Q_tau = Q.sub( 4 ).collapse()
assigner = FunctionAssigner( Q, [Q_z, Q_omega, Q_mu, Q_rho, Q_tau] )
class z_exact_expression( UserExpression ):
def eval(self, values, x):
values[0] = (x[0] ** 4 + x[1] ** 4) / 48.0
def value_shape(self):
return (1,)
class omega_exact_expression( UserExpression ):
def eval(self, values, x):
values[0] = (x[0] ** 3) / 12.0
values[1] = (x[1] ** 3) / 12.0
def value_shape(self):
return (2,)
class mu_exact_expression( UserExpression ):
def eval(self, values, x):
values[0] = (7 * x[0] ** 6 + 3 * x[0] ** 4 * x[1] ** 2 + 3 * x[0] ** 2 * x[1] ** 4 + 7 * x[1] ** 6) / 576.0
def value_shape(self):
return (1,)
class rho_exact_expression( UserExpression ):
def eval(self, values, x):
values[0] = x[0] * (7 * x[0] ** 4 + 2 * x[0] ** 2 * x[1] ** 2 + x[1] ** 4) / 96.0
values[1] = x[1] * (x[0] ** 4 + 2 * x[0] ** 2 * x[1] ** 2 + 7 * x[1] ** 4) / 96.0
def value_shape(self):
return (2,)
class f_exact_expression( UserExpression ):
def eval(self, values, x):
# values[0] = -16 * (np.cos( 4 * x[0] ) + np.cos( 4 * x[1] ) + np.sin( 2 * x[0] ) * np.sin( 2 * x[1] ))
values[0] = 1 / 8.0 * (3 * x[0] ** 4 + x[0] ** 2 * x[1] ** 2 + 3 * x[1] ** 4)
def value_shape(self):
return (1,)
# Define variational problem
psi = Function( Q )
nu_z, nu_omega, nu_mu, nu_rho, nu_tau = TestFunctions( Q )
z_exact = Function( Q_z )
omega_exact = Function( Q_omega )
mu_exact = Function( Q_mu )
rho_exact = Function( Q_rho )
tau_exact = Function( Q_tau )
f = Function( Q_z )
J_Q = TrialFunction( Q )
z, omega, mu, rho, tau = split( psi )
z_exact.interpolate( z_exact_expression( element=Q_z.ufl_element() ) )
omega_exact.interpolate( omega_exact_expression( element=Q_omega.ufl_element() ) )
mu_exact.interpolate( mu_exact_expression( element=Q_mu.ufl_element() ) )
rho_exact.interpolate( rho_exact_expression( element=Q_rho.ufl_element() ) )
tau_exact.interpolate( f_exact_expression( element=Q_tau.ufl_element() ) )
f.interpolate( f_exact_expression( element=Q_z.ufl_element() ) )
z_profile = Expression( '(pow(x[0], 4) + pow(x[1], 4)) / 48.0', element=Q.sub( 0 ).ufl_element() )
bc_z = DirichletBC( Q.sub( 0 ), z_profile, boundary )
assigner.assign( psi, [f, omega_exact, mu_exact, rho_exact, tau_exact] )
F_z = ((mu.dx( j )) * (nu_z.dx( j )) + f * nu_z) * dx
F_omega = (z * ((nu_omega[i]).dx( i )) + omega[i] * nu_omega[i]) * dx \
- n[i] * z * nu_omega[i] * ds
F_mu = (z * omega[i] * (nu_mu.dx( i )) + mu * nu_mu) * dx \
- n[i] * z * omega_exact[i] * nu_mu * ds
F_rho = (mu * ((nu_rho[i]).dx( i )) + rho[i] * nu_rho[i]) * dx \
- n[i] * mu * nu_rho[i] * ds
F_tau = (tau * nu_tau + rho[i] * (nu_tau.dx( i ))) * dx \
- n[i] * rho[i] * nu_tau * ds
F = F_omega + F_z + F_mu + F_rho + F_tau
bcs = [bc_z]
J = derivative( F, psi, J_Q )
problem = NonlinearVariationalProblem( F, psi, bcs, J )
solver = NonlinearVariationalSolver( problem )
solver.solve()
z_out, omega_out, mu_out, rho_out, tau_out = psi.split( deepcopy=True )
xdmffile_check.write( project( tau_out - f, Q_z ), 0 )
Please run it with python3 solve.py.
Thank you! 