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 <
ajles...@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.
>
