How to write the weak form for a 3D problem with a product of vectors and matrices?

1

Hello,

I do not know how to approach a simple 3D problem with FEniCS although I have it figured it out on paper.

I found an expression for the weak form of a particular problem, where a term of that weak form must be of the form “a^T B a” where a is a vector, and B a 3x3 matrix (or rank 2 tensor). The ^T denotes the transpose I’d like to implement it correctly into FEniCS. Note that on paper it’s the product of “matrices” of dimension (1x3) x (3x3) x (3x1) = (1x1), i.e. a scalar.

First of all, a is not given, it must be computed as B^-1 grad(V) - C grad(W) where ^-1 denotes the inverse, V and W are scalar valued functions. B and C are explicitly given.

In my code I have

from fenics import *
import numpy as np

(skipped the mesh and boundary conditions)

B = as_matrix(np.array([[10.0, 0, 0], [0, 10.0, 0], [0, 0, 23.0]]))

C = as_matrix(np.array([[18.0, 0, 0], [0, 18.0, 0], [0, 0, 3.0]]))

u, v = TestFunctions(ME)
VW = Function(ME)
V, W = split(VW)

weak_form = (inv(B)* grad(V) - C * grad(W)) * B * (inv(B) * grad(V) - C * grad(W))

But that line yields the error “UFLException: Invalid ranks 1 and 2 in product.”

I have tried many, many variants (including PETsc matrices, but I did something wrong because the program would not halt). For example using as_tensor() instead of as_matri). I tried to use inner() instead of dot(). I also tried transpose(1st term) * the other 2 terms. I tried nesting innner() and dot(), as to compute first a^T * B and then a^T * B * a. I tried to use nabla_grad instead of grad, etc. All my attempts have failed, some of them returned 0 hit in Google.

Examples of errors were “UFLException: Shapes do not match: <ListTensor id=139896187419224> and <Grad id=139896175492936>.” And “AttributeError: ‘Expression2LatexHandler’ object has no attribute ‘visit’”.

I am really at a loss. I do not know how to approach this problem with FEniCS even though on paper I know how each term should be defined and computed.

2

Posting a minimal working example would have been useful.

Anyway, this seems to work

from fenics import *
mesh = UnitSquareMesh(10,10)
ME = VectorFunctionSpace(mesh, 'CG', 1)
B = as_matrix([[10.0, 0, 0], [0, 10.0, 0], [0, 0, 23.0]])
C = as_matrix([[18.0, 0, 0], [0, 18.0, 0], [0, 0, 3.0]])
u, v = TestFunctions(ME)
VW = Function(ME)
V, W = split(VW)
a = inv(B)* grad(V) - C * grad(W)
weak_form = inner(B*a,a)
3

Thanks, indeed I should have included a MWE. In my code P1 = FiniteElement(“Lagrange”, mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

So I tried to replace the ME line with ME = VectorFunctionSpace(mesh, P1*P1, 3), but this yields UFLException: Don’t know how to split tensor valued mixed functions without flattened index space. I tried with a 1 and 2 also, instead of a 3. This always return the same error.

4

This also works

from fenics import *
mesh = UnitSquareMesh(10,10)
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
B = as_matrix([[10.0, 0, 0], [0, 10.0, 0], [0, 0, 23.0]])
C = as_matrix([[18.0, 0, 0], [0, 18.0, 0], [0, 0, 3.0]])
u, v = TestFunctions(ME)
VW = Function(ME)
V, W = split(VW)
a = inv(B)* grad(V) - C * grad(W)
weak_form = inner(B*a,a)
1 Like
5

Thank you very much, it indeed works!
Although I understand why it would work under this form, I do not understand why it didn’t work with several of my attempts.

The way I see it, B*a is a vector multiplied by a matrix, say a (1x3) vector multiplied by a (3x3) matrix. This yields a (1x3) vector, call it “aux”. Then you apply inner(aux, a) where “a” is a (3x1) vector. This yields a scalar which is what it should.

It seems that the product between vectors must be used with inner(), while multiplication between vectors and matrices must be done with *. But when I check the UFL documentation, this is not what I get at all…

So overall, I have no idea why the above code works and why my attempts failed.
Why has nesting 2 inner() failed for example.