Nicer way to write Stress Tensor

Main Forum Old-Forum - do not post here
1

Hello,

I have the following code (working) to calculate forces using the Maxwell Stress Tensor. Is there a more elegant way to write the tensor?

sigma_cf = nu * (OuterProduct(B, B) - 0.5 * (B*B) * Id(3) )

fes = HDiv(object.mesh, order=1)
sigma_1 = GridFunction(fes)
sigma_2 = GridFunction(fes)
sigma_3 = GridFunction(fes)

# columns of the tensor
sigma_1.Set((sigma_cf[0,0],sigma_cf[1,0],sigma_cf[2,0]))
sigma_2.Set((sigma_cf[0,1],sigma_cf[1,1],sigma_cf[2,1]))
sigma_3.Set((sigma_cf[0,2],sigma_cf[1,2],sigma_cf[2,2]))
    
Div_sigma_1 = div(sigma_1)
Div_sigma_2 = div(sigma_2)
Div_sigma_3 = div(sigma_3)   

gfforce = GridFunction(fes)
gfforce.Set((Div_sigma_1,Div_sigma_2,Div_sigma_3))

Force = Integrate(gfforce, mesh, VOL, definedon=mesh.Materials("plunger"))
2

Hi creativeworker,

NGSolve has a symmetric matrix-valued space called HDivDiv, which also has a divergence operator.
A difference to three copies of HDiv is that only the normal-normal component is continuous instead of the normal component. If you just need the space to have the divergence operator at hand not requiring the continuity (which would get enforced by the Set command if the discontinuous Flag is set False) you can simply your code:

sigma_cf = nu * (OuterProduct(B, B) - 0.5 * (B*B) * Id(3) )

fes = HDivDiv(object.mesh, order=1)#,discontinuous=True)
sigma = GridFunction(fes)
sigma.Set(sigma_cf)

Force = Integrate(div(sigma), mesh, VOL, definedon=mesh.Materials("plunger"))

If you really need the normal continuity (like for a-posteriori error estimators) you should stick with your code. However, the last three lines can be simplied to

Force = Integrate(CoefficientFunction( (Div_sigma_1,Div_sigma_2,Div_sigma_3) ), mesh, VOL, definedon=mesh.Materials("plunger"))

Another possibility would be to compute the divergence symbolic by hand and use the components Grad( B )

sigma_cf = nu * (OuterProduct(B, B) - 0.5 * (B*B) * Id(3) )
div_sigma_cf = nu* CoefficientFunction( (Terms of Grad(B)...) )
Force = Integrate(div_sigma_cf, mesh, VOL, definedon=mesh.Materials("plunger"))

Best,
Michael