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:
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
Code:
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 )
Code:
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
