I have a PDE of form
div(D(u) * E) = g
where D(u) is a 3x3 matrix function of u(x), and E is a vector in R^3. I need to solve for u(x) but I’m not sure how to implement a matrix of TrialFunctions in my bilinear/linear forms. Does Fenics have this capability, and, if so, how would I go about doing this? I can’t find anything online or in resources.
You can use the as_matrix
function to build up arbitrary matrices in UFL, as demonstrated below:
from dolfin import *
# Set up mesh, function space, and BC:
N = 8
mesh = UnitCubeMesh(N,N,N)
V = FunctionSpace(mesh,"CG",1)
bc = DirichletBC(V,Constant(0),"on_boundary")
# Example problem data:
def D(u):
# Gradient of `u` along diagonal of a matrix:
return -as_matrix([[u.dx(0), 0 , 0 ],
[0 , u.dx(1), 0 ],
[0 , 0 , u.dx(2)]])
E = Constant((1,1,1))
g = Constant(1)
# Set up weak form and solve:
u = TrialFunction(V)
v = TestFunction(V)
# Note: With these definitions, `D(u)*E` is
# equivalent to `-grad(u)`, so we're just solving the
# Poisson equation with source term `g`.
a = -dot(D(u)*E,grad(v))*dx # (Integrated `div` by parts)
L = g*v*dx
uh = Function(V)
solve(a==L,uh,bc)
# Output results:
File("u.pvd") << uh
3 Likes
Thank you! This seems to work just fine!