Thanks, the MatrixValued(L2(mesh)) helped me out. I attach the minimal example for other users!
Code:
[code]# -*- coding: utf-8 -*-
from ngsolve import *
from netgen.geom2d import SplineGeometry
# geometry
geo = SplineGeometry()
geo.AddRectangle(p1 = (-10,-5), p2 = (10, 5), leftdomain = 1, rightdomain = 0, bc = "outer")
geo.AddRectangle(p1 = (-5,-3), p2 = (5,3), leftdomain = 2, rightdomain = 1, bc = 'gamma')
geo.SetMaterial(1, 'mat1')
geo.SetMaterial(2, 'mat2')
ngmesh = geo.GenerateMesh(maxh = 5)
mesh = Mesh(ngmesh)
#%% Jacobian
G = specialcf.JacobianMatrix(2)
print('Jacobian Dims: ', G.dims.NumPy())
G_inv = Inv(G)
print('Jacobian Inv Dims: ', G_inv.dims.NumPy())
G_det = Det(G)
print('Jacobian det dims: ', G_det.dim)
G_det_inv = 1/G_det
mesh.UnsetDeformation()
gfu_G = GridFunction(MatrixValued(L2(mesh)))
gfu_G.Set(G)
gfu_G_det = GridFunction(L2(mesh))
gfu_G_det.Set(G_det)
#% Finite displacement
vdisplacement = 0.5
fes_structure = H1(mesh, order=1, dim=mesh.dim)
gfu_struc = GridFunction(fes_structure)
gfu_struc.Set((vdisplacement,0), definedon=mesh.Materials('mat2'))
mesh.SetDeformation(gfu_struc)
gfu_G2 = GridFunction(MatrixValued(L2(mesh)))
gfu_G2.Set(G)
gfu_G2_det = GridFunction(L2(mesh))
gfu_G2_det.Set(G_det)
mesh.UnsetDeformation()
Gdx = (gfu_G2-gfu_G)/vdisplacement
G_det_dx = (gfu_G2_det-gfu_G_det)/vdisplacement
gfu_Gdet_dx = GridFunction(L2(mesh))
gfu_Gdet_dx.Set(G_det_dx)
#
vtk = VTKOutput(ma = mesh,coefs=[gfu_G_det, gfu_G2_det, gfu_Gdet_dx, gfu_G, Gdx],names=['G_det','G_det2','G_detdx','G','Gdx'],filename = "simple_rectange",subdivision=0, legacy=True)
vtk.Do()