# Solving a variational form with a third-order tensor

I’m trying to solve a problem with a third-order tensor.
For the MRE, assume it’s just 9 times the Poisson equation.
MRE:

``````from ngsolve import  *
from ngsolve.meshes import MakeStructured3DMesh

res = 3
mesh = MakeStructured3DMesh(False, nx=res,ny=res,nz=res, mapping=lambda x,y,z : (-1+2*x,-1+2*y,-1+2*z))

fes = MatrixValued(H1(mesh, order=order, dirichlet="left|right"), symmetric = False)
(S), (T) = fes.TnT()

S[0,0].Diff(x), S[0,0].Diff(y), S[0,0].Diff(z),
S[0,1].Diff(x), S[0,1].Diff(y), S[0,1].Diff(z),
S[0,2].Diff(x), S[0,2].Diff(y), S[0,2].Diff(z),
S[1,0].Diff(x), S[1,0].Diff(y), S[1,0].Diff(z),
S[1,1].Diff(x), S[1,1].Diff(y), S[1,1].Diff(z),
S[1,2].Diff(x), S[1,2].Diff(y), S[1,2].Diff(z),
S[2,0].Diff(x), S[2,0].Diff(y), S[2,0].Diff(z),
S[2,1].Diff(x), S[2,1].Diff(y), S[2,1].Diff(z),
S[2,2].Diff(x), S[2,2].Diff(y), S[2,2].Diff(z),
), dims=(3,3,3) )

T[0,0].Diff(x), T[0,0].Diff(y), T[0,0].Diff(z),
T[0,1].Diff(x), T[0,1].Diff(y), T[0,1].Diff(z),
T[0,2].Diff(x), T[0,2].Diff(y), T[0,2].Diff(z),
T[1,0].Diff(x), T[1,0].Diff(y), T[1,0].Diff(z),
T[1,1].Diff(x), T[1,1].Diff(y), T[1,1].Diff(z),
T[1,2].Diff(x), T[1,2].Diff(y), T[1,2].Diff(z),
T[2,0].Diff(x), T[2,0].Diff(y), T[2,0].Diff(z),
T[2,1].Diff(x), T[2,1].Diff(y), T[2,1].Diff(z),
T[2,2].Diff(x), T[2,2].Diff(y), T[2,2].Diff(z),
), dims=(3,3,3) )

F = CF( (1,1,1,
1,1,1,
1,1,1), dims=(3,3) )

a = BilinearForm(fes, symmetric=True, symmetric_storage=True, condense=False)

f = LinearForm(fes)
f += (InnerProduct(T, F))*dx

sol = GridFunction(fes)

sol.vec[:] = 0

r = sol.vec.CreateVector()
w = sol.vec.CreateVector()

f.Assemble()
a.Assemble()
inv = a.mat.Inverse(fes.FreeDofs(a.condense), inverse="umfpack")

a.Apply(sol.vec, r)

r.data -= f.vec
if a.condense:
r.data += a.harmonic_extension_trans * r
w.data = inv * r
if a.condense:
w.data += a.harmonic_extension * w
w.data += a.inner_solve * r
sol.vec.data -= w
``````

For some reason this gives the error:

``````NgException                               Traceback (most recent call last)
<ipython-input-171-2a6ac94913c7> in <module>
2     f.Assemble()
3     a.Assemble()
----> 4     inv = a.mat.Inverse(fes.FreeDofs(a.condense), inverse="umfpack") #umfpack,sparsecholesky
5
6     a.Apply(sol.vec, r)

NgException: UmfpackInverse: Numeric factorization failed.
``````

Clearly, the Poisson problem is well-posed such that I don’t understand why umfpack fails.
I would be grateful for any insight!

Hi,
your “.Diff(x)” evaluates the CoefficientFunction-Tree. As it does not contain “x”, it yields zero. The trial and test functions depend on “x” through the finite elements, but as a `CoefficientFunction` this dependency is not recognized by `.Diff(...)` as the `CoefficientFunction` `x` (or `y` or `z`) are not involved in the definition of `S` or `T`.

What you want to do is taking the gradient `grad` of S and T and access its components to associate these to the components of `GradS` and `GradT`.

Best, Christoph

1 Like

That worked out!, thanks a lot!
Here’s my solution for future inquiries (gradient of a symmetric matrix):

``````fesH1 = H1(mesh, order=order, dirichlet="top")
fes = fesH1*fesH1*fesH1*fesH1*fesH1*fesH1
(S00,S01,S02,S11,S12,S22), (T00,T01,T02,T11,T12,T22) = fes.TnT()

), dims=(3,3,3) )

``````fes = MatrixValued(H1(mesh, order=order, dirichlet="left|right"), symmetric = False)