Trying to resolve an error with expressions of different shapes

2

It’s difficult to decipher how you derived the Jacobian but rhoJ doesn’t look dimensionally correct. Consider this post.

With the following code:

from fenics import *
mesh  = UnitSquareMesh(20, 20)
epsi = Constant(20)

V = VectorFunctionSpace(mesh, 'P', 1)
Q = FunctionSpace(mesh, 'Lagrange', 1)

E = Function(V)

rho_2 = TrialFunction(Q)
u = TestFunction(Q)
rho_2 = Function(Q)

b = ((2*rho_2 )/ epsi)*dx
b_compiled = fem.form.Form(b)
print(b_compiled.rank())

a = dot(E,grad(grad(rho_2)))*dx
a_compiled = fem.form.Form(a)
print(a_compiled.rank())

The first integrand returns a rank of 0 while the second:

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

The grad(grad(rho_2)) generates a tensor of rank 2 which can’t be summed with your other term of rank 0.

1 Like