I use RT element to solve problem A and denote the finite element solution by sigma, sigma is a matrix-valued finite element function. I need to use the divergence of sigma to be the source term of another problem, the trial function space is vector-valued functions. But it doesn't work. The code is
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(8,8)
RT = VectorElement("RT",mesh.ufl_cell(), 1,2)
Q = FunctionSpace(mesh, RT)
V = VectorFunctionSpace(mesh, 'CG',1)
f = Expression((('x[0]*x[1]','x[0]'),('x[0]','x[0]*x[1]')),degree=2)
sigma = TrialFunction(Q)
tau = TestFunction(Q)
a = inner(sigma, tau)*dx
L = inner(f,tau)*dx
sigma = Function(Q)
solve(a == L, sigma)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, Constant((0.0, 0.0)), boundary)
u = TrialFunction(V)
v = TestFunction(V)
au = inner(grad(u), grad(v))*dx
Lu = inner(sigma,grad(v))*dx
u = Function(V)
solve(au == Lu, u, bc)
