getfem-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

implementing linear viscoelasticity in getfem


From: Lesage,Anne Cecile J
Subject: implementing linear viscoelasticity in getfem
Date: Fri, 11 Mar 2022 01:43:24 +0000

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.

Attachment: NotesLinearViscoElasticGFem.pdf
Description: NotesLinearViscoElasticGFem.pdf


reply via email to

[Prev in Thread] Current Thread [Next in Thread]