How To Matrix of TestFunctions

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)

# Output results:
File("u.pvd") << uh

Thank you! This seems to work just fine!