Solving a variational form with a third-order tensor

NGSolve
1

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()

GradS = CF( (
    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) )

GradT = CF( (
    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)
a += ( InnerProduct(GradT, GradS) )*dx

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

sol = GridFunction(fes)

sol.vec[:] = 0

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

with TaskManager():
    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!

2

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
3

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()

GradS = CF( (
    grad(S00),
    grad(S01),
    grad(S02),
    grad(S01),
    grad(S11),
    grad(S12),
    grad(S02),
    grad(S12),
    grad(S22)
            ), dims=(3,3,3) )

GradT = CF( (
    grad(T00),
    grad(T01),
    grad(T02),
    grad(T01),
    grad(T11),
    grad(T12),
    grad(T02),
    grad(T12),
    grad(T22)
            ), dims=(3,3,3) )
4

you get the grad of the matrix-entries also via

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

as a third order tensor. This preserves more structure, and NGSolve can better optimize.

Joachim