Nicer way to write Stress Tensor

creativeworker
Topic Author
Offline
Senior Member

4 years 1 month ago #3177 by creativeworker

Nicer way to write Stress Tensor was created by creativeworker

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?

Code:

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

mneunteufel
Offline
Premium Member

4 years 1 month ago #3181 by mneunteufel

Replied by mneunteufel on topic Nicer way to write Stress Tensor

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

The following user(s) said Thank You: creativeworker

Time to create page: 0.090 seconds