|
From: | Konstantinos Poulios |
Subject: | Re: implementing linear viscoelasticity in getfem |
Date: | Fri, 11 Mar 2022 09:20:38 +0100 |
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.
[Prev in Thread] | Current Thread | [Next in Thread] |