Re: implementing linear viscoelasticity in getfem

[Konstantinos Poulios]

Dear Anne-Cecile

To implement Eq. (12) from your notes (btw. t on the left hand side should be t_{n+1}), you need the following:

1) define an im_data structure for storing the components of a 3x3 tensor

mimdb3x3 = gf.MeshImData(mimu3b, -1, [3,3])

(this is a bit of wasting of memory in your case because for small strains, only 6 out of the 9 components are needed, but it makes things simpler)

2) storage for \sigma^{visc}(t_n) and H(t_n)

md.add_im_data("Svprev", mimdb3x3)

md.add_im_data("Hprev", mimdb3x3)

3) macros corresponding to Eq. (9) (I guess it should be the symmetric part of the gradient of u) and Eq. (12)

md.add_macro("H(u)", "Sym(Grad(u))+1/(1-2*nub)*Div(u)*Id(3)") # note the Sym() and Id(3) added here

md.add_macro("Sv(u)", "Svprev-beta*(G0-Ginf)*dt*0.5*(Hprev+exp(-beta*dt)*H(u))")

4) use the Sv macro in your virtual work _expression_

md.add_linear_term(mimu3b, 'Sv(ub):Grad(Test_ub)')

this already contains the volumetric part. The explanation is because Id(3):Grad(Test_ub)= Div(Test_ub)

5) an update mechanism for your stored quantities after each iteration

md.set_variable("Svprev", md.interpolation("Sv(ub)", mimdb3x3, -1)) # it is important that this is done before the next line, because the result depends on Hprev

md.set_variable("Hprev", md.interpolation("H(ub)", mimdb3x3, -1))

are memory you can also avoid storing Hprev and Sprev and recalculate them ev

If you want to spare memory you can also avoid storing Hprev and Sprev and recalculate them every time from ub_prev, but if memory is not an issue I would recommend you the solution above.

Extra notes:
- decide about which symbol to use for shear modulus G or mu and stick to one symbol

- for small strains, your \nabla u in Eq. (2) should be \nabla_s u, so that your stress tensor is symmetric

- for your Zener model you do not need to split to volumetric and deviatoric parts (H1 and H2) as was the case in my notes

Best regards

Kostas

On Fri, Mar 11, 2022 at 2:44 AM Lesage,Anne Cecile J <AJLesage@mdanderson.org> wrote:

Dear all

 

I prepared those notes to adapt Pr Poulios’s script on Quasi nonlinear viscoelastic to linear viscoelastic

Please find some important part of the script

I do not know how to define Sv2 in my script it is the second part of the viscous stress

 

Thank you

BR

Anne-Cecile

 

Please find the script enclosed in the email

                    

 

mfub = gf.MeshFem(meshb, 3)

mfub.set_classical_fem(disp_fem_order)

mimu3b = gf.MeshIm(meshb, 5)   # integration displacement on tetrahedra brain

 

mfpb = gf.MeshFem(meshb, 1)

mfpb.set_classical_fem(press_fem_order)

mimp3b = gf.MeshIm(meshb, 5)   # integration pressure on tetrahedra brain

 

mimd3x3 = gf.MeshImData(mimp3b, -1., [3,3])                                      

 

md.add_fem_variable("ub", mfub)      # displacements field brain. poroelastic

md.add_fem_data("ub_prev", mfub)

md.add_fem_data("Svprev2", mfpb)

md.add_fem_variable("Sv2",mfpb)

 

md.add_initialized_data("G0", G0)

md.add_initialized_data("Ginf", Ginf)

md.add_initialized_data("nub", Nub)

md.add_initialized_data("beta", beta)

 

…

 

md.add_macro("H1","Grad_ub")

md.add_macro("H2","1.0/(1-2*nub)*Div(ub)")

md.add_macro("Hprev1","Grad_ub_prev")

md.add_macro("Hprev2","1.0/(1-2*nub)*Div(ub_prev)")

 

md.add_im_data("Svprev1", mimd3x3)   # Viscous  stress at previous timestep

md.add_macro("Sv1","Svprev1-beta*(G0-Ginf)*dt*0.5*(Hprev1+exp(-beta*dt)*H1)")

 

#md.add_im_data("Svprev2", mimp3b)   # Viscous  stress at previous timestep

md.add_macro("Sv2","Svprev2-beta*(G0-Ginf)*dt*0.5*(Hprev2+exp(-beta*dt)*H2)")

 

steps = ma.floor(tmax/dt)

steps = 4

 

# brain elastic

md.add_initialized_data('cmub', Mub)

md.add_initialized_data('clambdab', Lambdab)

 

assemelas = 3

 

if assemelas ==1:

   md.add_isotropic_linearized_elasticity_brick(mimu3b, 'ub', 'clambdab', 'cmub')

elif assemelas==2:

   md.add_linear_term(mimu3b, 'G0*Grad(ub):Grad(Test_ub)+G0/(1-2*nub)*Div(ub)*Div(Test_ub)')

elif assemelas==3:

   md.add_linear_term(mimu3b, 'G0*H1:Grad(Test_ub)+G0*H2*Div(Test_ub)')

  

viscous = 1

if viscous ==1:

   md.add_linear_term(mimu3b, 'Sv1:Grad(Test_ub)+Sv2*Div(Test_ub)')

 

 

if gravdir=='-X':

   print('Assembling buoyancy term along -ex\n');

   md.add_linear_term(mimu3b, 'Test_ub(1)*g*(rhot-rhoa*Heaviside(X(1)-csflevel)-rhow*Heaviside(-X(1)+csflevel))')

elif gravdir=='X':

   print('Assembling buoyancy term along ex\n');

   md.add_linear_term(mimu3b, 'Test_ub(1)*g*-1.0*(rhot-rhoa*Heaviside(-X(1)+csflevel)-rhow*Heaviside(X(1)-csflevel))')

elif gravdir=='-Z':

   print('Assembling buoyancy term along -ez\n');

   md.add_linear_term(mimu3b, 'Test_ub(3)*g*(rhot-rhoa*Heaviside(X(3)-csflevel)-rhow*Heaviside(-X(3)+csflevel))')

     

      

                                        

print('Solve problem with ', md.nbdof(), ' dofs')

time_elapsed(timer)

 

mass_mat = gf.asm_mass_matrix(mimp3b, mfpb)

 

######################################

# track landmarks

pos_tracked = gf.Slice("points", meshb, preopf)

 

 

 

for step in range(steps):

   nit, conv = md.solve('noisy', 'max_iter', 25, 'max_res', 1e-10, 'lsolver', 'mumps',

                                'lsearch', 'simplest', 'alpha max ratio', 100, 'alpha min', 0.2, 'alpha mult', 0.6,

                        'alpha threshold res', 1e9)

   time_elapsed(timer)    

   it = step+1

     

   mfoutb.export_to_vtu("LVE%sbrain_%i.vtu" % (patient,it),

                      mfub, md.variable("ub"), "Displacements")

   #np.savetxt("brain_displacements.txt",md.variable("ub"))

   #mfouth.export_to_vtu("PoroEResect25head_%i.vtu" % step,

   #                   mfuh, md.variable("uh"), "Displacements")

   md.set_variable("ub_prev", md.variable("ub"))

   md.set_variable("Svprev1", md.interpolation("Sv1", mimd3x3, -1))

   md.set_variable("Svprev2", md.variable("Sv2"))

  u_tracked = gf.compute_interpolate_on( mfub, md.variable("ub"), pos_tracked)

….

 

 

The information contained in this e-mail message may be privileged, confidential, and/or protected from disclosure. This e-mail message may contain protected health information (PHI); dissemination of PHI should comply with applicable federal and state laws. If you are not the intended recipient, or an authorized representative of the intended recipient, any further review, disclosure, use, dissemination, distribution, or copying of this message or any attachment (or the information contained therein) is strictly prohibited. If you think that you have received this e-mail message in error, please notify the sender by return e-mail and delete all references to it and its contents from your systems.