Is there a way to access the JacobianInverse as specialcf in python?
Hi,
you can invert the specialcf.JacobianMatrix to obtain the inverse
F_inv = Inv(specialcf.JacobianMatrix(3))
Best,
Michael
Perfect. Is there a way to Diff a JacobianMatrix?
The goal is to implement the virtual work principle: An original stationary method using local jacobian derivative for direct finite element computation of electromagnetic force, torque and stiffness - ScienceDirect
It should be possible to interpolate the Jacobian (or some components of it) into an L2/VectorL2/MatrixValued(L2) GridFunction and then compute the gradient.
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()
[/code][/code]
I implemented the Local Jacobian Virtual Work Method, but the results are not clear to me.
Does anybody have an idea where things go wrong? Is it the finite differences around line 100? It seemed reasonable good in the minimal example.
# -*- coding: utf-8 -*-
# Johannes Maierhofer
# 2022-11-14
# J.L. Coulomb Example Two Rings
# "Finite Element Implementation of Virtual Work Principle for Magnetic or Electric Force and Torque Computation"
# J.L. Coulomb and G. Meunier
# 10.1109/TMAG.1984.1063232
from ngsolve import *
from netgen.csg import *
# Configuration
R1 = 50e-3 #m
R2 = 60e-3 #m
h = 10e-3 #m
d = 40e-3 #m
j = 3e6 #A/m^2
Fexact = 78e-3 #N
#%%
geom = CSGeometry()
Ring1 = ((Cylinder(Pnt(0,0,d/2), Pnt(0,0,d/2+h), R2)-Cylinder(Pnt(0,0,d/2), Pnt(0,0,d/2+h), R1)) * Plane(Pnt(0,0,d/2), Vec(0,0,-1)) * Plane(Pnt(0,0,d/2+h), Vec(0,0,1))).bc("gamma_ring1")
Ring2 = ((Cylinder(Pnt(0,0,-d/2), Pnt(0,0,-d/2-h), R2)-Cylinder(Pnt(0,0,-d/2), Pnt(0,0,-d/2-h), R1)) * Plane(Pnt(0,0,-d/2), Vec(0,0,1)) * Plane(Pnt(0,0,-d/2-h), Vec(0,0,-1))).bc("gamma_ring2")
r = 0.1
box = OrthoBrick (Pnt(-r,-r,-r), Pnt(r,r,r)).bc('outer')
geom.Add (Ring1.mat("mat1").maxh(0.005))
geom.Add (Ring2.mat("mat2").maxh(0.005))
geom.Add ( (box-Ring1-Ring2).mat("air"))
mesh = Mesh(geom.GenerateMesh(maxh=0.01))
mesh.Curve(3)
domain=mesh.MaterialCF({'air':0, 'mat1':1, 'mat2':2})
fes = HCurl(mesh, order=1, nograds=True, dirichlet='outer')
print ("ndof = ", fes.ndof)
u,v = fes.TnT()
#%%
mu0 = 1.257e-6
mur = {"air" : 1, "mat1" : 1, "mat2":1 }
nu = 1.0 / mu0 / mesh.MaterialCF(mur)
a = BilinearForm(fes)
a += SymbolicBFI (nu*curl(u)*curl(v) + nu*1e-6*u*v)
c = Preconditioner(a, "bddc")
a.Assemble()
f = LinearForm(fes)
tau = CoefficientFunction((y, -x, 0))
tau = 1.0/Norm(tau)*tau
curdens_cf = mesh.MaterialCF({'mat1':j*tau, 'mat2':j*tau}, default=(0,0,0))
f += SymbolicLFI(curdens_cf*v)
f.Assemble()
gfu = GridFunction(fes)
solvers.CG(mat=a.mat, pre=c.mat, rhs=f.vec, sol=gfu.vec)
# Reference-Force from Result-Sheet
print(f'\n Reference Force (attracting): {Fexact:.3g} N')
#%% virtual work with derivatives
G = specialcf.JacobianMatrix(3)
G_inv = Inv(G) # --> same as fem.Inv, officially possible for JacobianMatrix
G_det = Det(G) # --> same as fem.Det
G_det_inv = 1/G_det
B = curl(gfu)
L2_B = GridFunction(VectorL2(mesh))
L2_B.Set(B) # not necessary, nearly same result by directly using B in line 105
mesh.UnsetDeformation()
gfu_G_det = GridFunction(L2(mesh))
gfu_G_det.Set(G_det)
gfu_G = GridFunction(MatrixValued(L2(mesh)))
gfu_G.Set(G)
#% Virtual displacement
vdisp = 0.0001
vdisplacement = CoefficientFunction((0,0,vdisp))
# Virtual displacement
gfu_struc = GridFunction(VectorH1(mesh))
gfu_struc.Set(vdisplacement, definedon=mesh.Materials('mat1'))
mesh.SetDeformation(gfu_struc)
gfu_G2_det = GridFunction(L2(mesh))
gfu_G2_det.Set(G_det)
gfu_G2 = GridFunction(MatrixValued(L2(mesh)))
gfu_G2.Set(G)
mesh.UnsetDeformation()
# Find |G|/ds and G/ds by finite differences
G_det_ds = (gfu_G2_det-gfu_G_det)/Norm(vdisplacement)
G_ds = (gfu_G2-gfu_G)/Norm(vdisplacement)
F_vw = Integrate(-nu*L2_B*(G_inv*G_ds*L2_B) + 0.5*nu*L2_B*L2_B*G_det_inv*G_det_ds, mesh=mesh, definedon=mesh.Materials('air'))
print(f'Force - virtual work - local Jacobian: {F_vw:.3g} N')
Fe = Integrate(-nu*B*(G_inv*G_ds*B) + 0.5*nu*B*B*G_det_inv*G_det_ds, mesh=mesh, element_wise=True)
gfu_Fe_L2 = GridFunction(L2(mesh))
gfu_Fe_L2.vec.FV().NumPy()[:]=Fe
gfu_Fe = GridFunction(VectorL2(mesh))
gfu_Fe.Set(gfu_Fe_L2*vdisplacement/Norm(vdisplacement))
#%%
vtk = VTKOutput(ma = mesh,coefs=[domain,curdens_cf,B, L2_B,G_det_ds,gfu_Fe],names=['domain','curdens','B','L2_B','G_det_ds','F_element'],filename = "JLCoulombRings",subdivision=0, legacy=True)
vtk.Do()
#%% Virtual work Finite Differences
sysE = BilinearForm(fes)
sysE += SymbolicBFI(nu*curl(u)*curl(v))
W1 = sysE.Energy(gfu.vec)
# Virtual displacement
vdisplacement = 0.000001
gfu_struc = GridFunction(VectorH1(mesh))
gfu_struc.Set((0,0,vdisplacement), definedon=mesh.Materials('mat1'))
mesh.SetDeformation(gfu_struc)
# second energy calculation
W2 = sysE.Energy(gfu.vec)
mesh.UnsetDeformation()
F = (W2-W1)/vdisplacement
print(f'Force - virtual work - finite differences: {F:.3g} N')