|
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.
|
NotesLinearViscoElasticGFem.pdf
Description: NotesLinearViscoElasticGFem.pdf
[Prev in Thread] | Current Thread | [Next in Thread] |