# Transient problems in elasticity

**URL:** https://community.freefem.org/t/transient-problems-in-elasticity/863
**Category:** General Discussion
**Created:** [March 18, 2021, 11:08pm UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863 "2021-03-18T23:08:47Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![felric](https://avatars.discourse-cdn.com/v4/letter/f/5f8ce5/32.png) [@felric](https://community.freefem.org/u/felric)
#### Post date: [March 18, 2021, 11:08pm UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863/1 "2021-03-18T23:08:47Z")

</div>

Dear all,

first, many thanks for your wonderful software.  
I’m trying to solve transient problems in elasticity, more in particular simulation of mechanical resonators which have some damping (e.g. the time dependent response of a bell that is stroke at part of the boundary at time t=0, or simulation of ultrasonic transducers).  
I read the manual and looked for examples and also in the community blog but I wasn’t able to find any example of solutions of transient problems in elasticity.

I know that usually there are two ways to perform such task:

- direct time integration
- modal superposition

I guess that Freefem++ is able to solve such problems. Do you have any example or documentation which could help implementing such type of problems in Freefem++?  
Thanks a lot!!

Riccardo

---

<div class="post-metadata">

### Author: ![felric](https://avatars.discourse-cdn.com/v4/letter/f/5f8ce5/32.png) [@felric](https://community.freefem.org/u/felric)
#### Post date: [March 27, 2021, 9:09pm UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863/2 "2021-03-27T21:09:27Z")

</div>

Dear all,

I tried to build a code to solve transient problems in elasticity (a simple bar with a fixed end, gravity force and a periodic force on the tip) but either it gives nonsense results or freefem crashes.  
The time integration algorithm is Newmark (see e.g. Bathe, Finite Element Procedures, paragraph 9.2.3 “The Newmark method”).

I’m struggling to find any working solution.  
Please could you help me?  
Below the code.

Thanks a lot  
Riccardo

load “lapack”

//Parameters  
real dt=0.00001;  
real T=0.001;

real Rho = 8000.; //Density  
real E = 210.e9; //Young modulus  
real Nu = 0.27; //Poisson ratio  
real Gravity = -9.81; //Gravity  
real ForceOnTip = 1e-10;

//Newmark constants  
real delta = 0.5;  
real alpha = 0.25 \* (0.5+delta);  
real alpha0=1/(alpha_dt^2);  
real alpha1=delta/(alpha_dt);  
real alpha2=1/(alpha_dt);  
real alpha3=1/(2.alpha)-1.;  
real alpha4=delta/alpha-1.;  
real alpha5=dt/2.(delta/alpha-2.);  
real alpha6=dt_(1-delta);  
real alpha7=delta\*dt;

//Mesh quality  
real L = 20.; //Beam length  
real H = 1.; //Beam height  
int Fixed = 1; //Beam fixed label  
int Free = 2; //Beam free label  
int Forced = 3; //Beam free label  
border b1(t=0., L){x=t; y=0.; label=Free;};  
border b2(t=0., H){x=L; y=t; label=Forced;};  
border b3(t=0., L){x=L-t; y=H; label=Free;};  
border b4(t=0., H){x=0.; y=H-t; label=Fixed;};

int nnL = 200; //max(2., nnL);  
int nnH = 10; //max(2., nnH);  
mesh Th = buildmesh(b1(nnL) + b2(nnH) + b3(nnL) + b4(nnH));

//Fespace  
func Pk = P2;  
fespace Uh(Th, [Pk, Pk]);  
Uh [ux, uy]=[0,0];  
real[int] UU(Uh.ndof), UUt(Uh.ndof), UUtt(Uh.ndof), temp(Uh.ndof);  
real[int] ElasticityForce(Uh.ndof), UUold(Uh.ndof), UUttold(Uh.ndof);  
UU=ux[]; // cause also the copy of uy  
UUt=0;

//Macro  
real sqrt2 = sqrt(2.);

macro Epsilon(ux, uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt2] //  
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //

//Problem  
real Mu = E/(2._(1.+Nu));  
real Lambda = E_Nu/((1.+ Nu)\*(1.-2.\*Nu));  
real freq = 20000;

varf vElasticity ([ux,uy], [vx, vy]) = int2d(Th)( Lambda \* Divergence(vx, vy) \* Divergence(ux, uy) + 2. \* Mu \* ( Epsilon(vx, vy)’ \* Epsilon(ux, uy) ) );  
varf vBC ([ux,uy], [vx, vy]) = on(Fixed, ux=0, uy=0) ;  
varf vMass ([ux,uy], [vx, vy]) = int2d(Th)( Rho \* (ux \* vx + uy \* vy)/(dt\*dt) );

matrix Mass = vMass(Uh, Uh, solver=sparsesolver);  
matrix Elasticity1 = vElasticity(Uh, Uh);  
matrix Elasticity2 = vBC(Uh, Uh);  
matrix Elasticity = Elasticity1 + Elasticity2; // I’m not sure if boundary conditions should be applied to Elasticity matrix  
matrix Damping = 0_Mass+0_Elasticity;

// modified stiffness matrix for Newmark algorithm  
matrix ElasticityStiffness = Elasticity + alpha0_Mass + alpha1_Damping;

// Solve static problem  
set (Elasticity, solver=sparsesolver);  
real[int] ElasticityBoundary = vElasticity(0, Uh);  
real[int] ElasticityBC = vBC(0, Uh);  
real[int] ElasticityBF = ElasticityForce; //ElasticityBoundary + ElasticityForce;  
real[int] uyy = Elasticity^-1 \* ElasticityBF;

// Prepare dynamic (transient problem)  
// computation of acceleration (UUtt)  
temp = Damping_UUt;  
temp += Elasticity_UU;  
temp = ElasticityForce - temp;  
UUtt = Mass^-1 \* temp;

for (real t=0;t\<T;t+=dt)  
{  
// Newmark algorithm  
varf vForce ([ux,uy], [vx, vy]) = int2d(Th)( Rho \* Gravity \* vy ) + on(Forced, ux=0, uy=ForceOnTip_sin(2_3.14159_freq_(t+dt)));

```
// computation of modified elasticity (stiffness) matrix 
ElasticityForce = vForce(0, Uh); 	
temp = alpha0*UU;
temp += alpha2*UUt;
temp += alpha3*UUtt;
ElasticityForce += Mass * temp;
temp = alpha1*UU;
temp += alpha4*UUt;
temp += alpha5*UUtt;
ElasticityForce += Damping * temp;

UUold = UU;
UUttold = UUtt;
UU = ElasticityStiffness^-1 * ElasticityForce;

temp = UU-UUold;
UUtt = alpha0 * temp;
temp = alpha2*UUt;
temp += alpha3*UUttold;
UUtt -= temp;

UUt += alpha6*UUttold;
UUt += alpha7*UUtt;

```

}  
Uh [uux, uuy] = UU;  
real dxmin = uux[].min;  
real dymin = uuy[].min;  
cout \<\< " - dep. max x = “\<\< dxmin \<\< " y=” \<\< dymin \<\< endl;  
cout \<\< " dep. (20, 0) = " \<\< uux(20, 0) \<\< " " \<\< uuy(20, 0) \<\< endl;

//[ux, uy] = UU;  
//Movemesh  
real coef = 2e11;  
mesh Th1 = movemesh(Th, [x+uux_coef, y+uuy_coef]);

[uux, uuy] = [uux, uuy];  
plot(Th1, [uux, uuy]);

---

<div class="post-metadata">

### Author: ![Markov](https://avatars.discourse-cdn.com/v4/letter/m/f1d935/32.png) [@Markov](https://community.freefem.org/u/Markov)
#### Post date: [July 1, 2021, 10:08am UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863/3 "2021-07-01T10:08:58Z")

</div>

Hello,

I am interested in the above problem, have you found a solution for it?

---

<div class="post-metadata">

### Author: ![felric](https://avatars.discourse-cdn.com/v4/letter/f/5f8ce5/32.png) [@felric](https://community.freefem.org/u/felric)
#### Post date: [July 7, 2021, 9:36pm UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863/4 "2021-07-07T21:36:29Z")

</div>

Dear Markov,

yes I did.  
I tried the Newmark algorithm from a well known reference in the field (the Bathe book “Finite Element Procedures”) but it did not work. However in the same chapter there was the description of an improvement made by Bathe, which he called “Bathe algorithm”.  
It is slightly more complex as it involves computation at half time steps too, but it is much more accurate so according to the Author it proved to be useful not only for the original task (non linear stress, if I remember) but also for general problems.  
Also my implementation of this algorithm was giving faulty results as diverging results were obtained.

However, after inspecting carefully the book, I discovered that there was a “step by step Example” that was supplied short afterwards, and checking carefully I found it had a different value for one constant.  
So I guessed that the algorithm as described had an error (that’s strange, as I was looking in the latest edition of the book). But the example worked like a charm! And in fact it using that value for the constant the algorithm worked really well in any condition!

This is my implementation. It is a 3D fully working example of a tube excited sinusoidally at one end and fixed at the other.

load “msh3”  
load “medit”  
load “lapack”  
load “iovtk”  
load “UMFPACK64”  
//load “PETSc”

real RhoTi=4430;  
real ETi=114e9;  
real NuTi=0.342;

real TuboL=1081.;  
real TuboDi=49.;  
real TuboDe=60.5;  
real TuboRi=TuboDi/2;  
real TuboRe=TuboDe/2;

real freqDelta=30; //hz  
real freqForce=7.500E-03_TuboL_TuboL - 3.141E+01\*TuboL + 4.651E+04;  
real freqstart=freqForce-freqDelta/2.;  
real freqend=freqForce+freqDelta/2.;  
real toffset=-4e-3;

real DampingMass=0;  
real DampingStiffness=25e-9;

int lbCera=110;  
int lbTubo=803;  
int lbTrasdEnd=106;

real MeshAspectRatio=3.;  
real NodesPerMeter=450.;  
real MagnificationOfDeformation=5.;  
real Gravity=0.;  
real ForceOnTip=50e5;  
real dt=0.000001;  
real T=dt\*(18000);  
string OutputFileName =“output\_ThTrasdTube\_” + “Di”+TuboDi_1000.+“De”+TuboDe_1000.+“L”+TuboL\*1000. +  
“dampM”+DampingMass+“dampK”+DampingStiffness+  
“freqStart”+freqstart+“freqEnd”+freqend+“offset”+(-toffset)+".txt";

real work=0.;  
real workmax=0.;  
real defradmax=0.;  
int NumberOfModesToCompute=5;

real Spessore=(TuboDe-TuboDi)/2.0;  
int Fixed = 100; //Beam fixed label  
int Free = 2; //Beam free label  
int Forced = 300; //Beam free label  
int[int] labs=[Fixed, 12, Forced, 13];

int NSteps;  
NSteps=T/dt+1;

border B801(t=0,1){x=(0)+(TuboL)\*t;y=(TuboRi)+(0)\*t;label=801;};  
border B802(t=0,1){x=(TuboL)+(0)\*t;y=(TuboRi)+(5.75)\*t;label=802;};  
border B803(t=0,1){x=(TuboL)+(-TuboL)\*t;y=(TuboRe)+(0)\*t;label=803;};  
border B804(t=0,1){x=(0)+(0)\*t;y=(TuboRe)+(-5.75)\*t;label=804;};

real nn=0.15;  
real h= 1./nn;

mesh Th2=buildmesh(  
B801(TuboL_nn) + B803(TuboL_nn) + B802(20_nn) + B804(20_nn)  
);

int MaxLayersT=8;//(int(2_pi_RR/h)/4)_4;  
func zminT = pi/2.;  
func zmaxT = 0;  
func fx= 0.001_y_sin(z);// / max( abs(cos(z) ), abs(sin(z)));  
func fy= 0.001_y_cos(z);// / max( abs(cos(z) ), abs(sin(z)));  
func fz= 0.001_x;

int[int] r1T=[0,0], r2T=[0,0,2,2];  
int[int] r4T=[0,2];  
int labelUp=999, labelDown=998;  
int[int] lbup = [0, labelUp, 6, labelUp, 8, labelUp, 17, labelUp, 23, labelUp, 25, labelUp, 28, labelUp, 35, labelUp];  
int[int] lbdown = [0, labelDown, 6, labelDown, 8, labelDown, 17, labelDown, 23, labelDown, 25, labelDown, 28, labelDown, 35, labelDown];  
mesh3 Th3T=buildlayers(Th2,coef= max(.1,y/min(max(20.,80-x),40.) ), MaxLayersT,zbound=[zminT,zmaxT],transfo=[fx,fy,fz],facemerge=0, labelup=lbup, labeldown=lbdown);

//  
//Bethe constants  
real delta = 0.5;  
real alpha = 0.25 \* (0.5+delta);  
real alpha0=16./(dt^2);  
real alpha1=4./(dt);  
real alpha2=9./(dt^2);  
real alpha3=3./dt;  
real alpha4=2.\*alpha1;  
real alpha5=12./(dt^2);  
real alpha6=-3./(dt^2);  
real alpha7=-1./(dt);

ofstream Outputfile(OutputFileName);  
Outputfile \<\< “Lunghezza:\t” \<\< TuboL \<\< endl;  
Outputfile \<\< “DiametroInterno:\t” \<\< TuboDi \<\< endl;  
Outputfile \<\< “TuboDe:\t” \<\< TuboDe \<\< endl;

//Problem

func Pk = P2;  
fespace Uh(Th3T, [Pk,Pk,Pk]);  
fespace U0(Th3T, P0);  
Uh [ux, uy, uz]=[0,0,0];

real sqrt2 = sqrt(2.);  
macro Epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM  
macro Divergence(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

U0 fRho = (RhoTi);

U0 fLambda = (ETi_NuTi/((1.+ NuTi)_(1.-2.\*NuTi)));

U0 fMu = (ETi/(2.\*(1.+NuTi)));

varf vElasticity ([ux,uy,uz], [vx,vy,vz], tgv=1e18) = int3d(Th3T)( fLambda \* Divergence(vx, vy, vz) \* Divergence(ux, uy, uz) + 2. \* fMu \* ( Epsilon(vx, vy, vz)’ \* Epsilon(ux, uy, uz) ) );  
varf vBC ([ux,uy,uz], [vx,vy,vz], tgv=1e18) = on(801, uz=0) +on(labelUp, ux=0)+on(labelDown, uy=0) ;  
varf vMass ([ux,uy,uz], [vx,vy,vz]) = int3d(Th3T)( fRho \* (ux \* vx + uy \* vy + uz \* vz) );

matrix Mass = vMass(Uh, Uh, solver=sparsesolver);  
matrix Elasticity1 = vElasticity(Uh, Uh);  
matrix Elasticity2 = vBC(Uh, Uh);  
matrix Elasticity = Elasticity1 + Elasticity2;  
matrix Damping = DampingMass_Mass+DampingStiffness_Elasticity;

real[int] ElasticityForce(Uh.ndof);  
real[int] UU(Uh.ndof), UUt(Uh.ndof), UUtt(Uh.ndof), temp(Uh.ndof);  
real[int] ForceStat(Uh.ndof), ForceHalf(Uh.ndof), ForceFull(Uh.ndof), UUold(Uh.ndof), UUtold(Uh.ndof), UUttold(Uh.ndof);  
real[int] UUhalf(Uh.ndof), UUthalf(Uh.ndof), UUtthalf(Uh.ndof);  
real[int] UUfull(Uh.ndof), UUtfull(Uh.ndof), UUttfull(Uh.ndof);  
UU=ux[]; // cause also the copy of uy  
UU=0;  
UUt=0;  
UUtt=0;

varf vForce0 ([ux,uy,uz], [vx, vy, vz]) = int2d(Th3T,804)(vz \* ForceOnTip);

// Solve static problem  
set (Elasticity, solver=sparsesolver);  
ElasticityForce = vForce0(0, Uh);   
real[int] ElasticityBoundary = vElasticity(0, Uh);  
real[int] ElasticityBC = vBC(0, Uh);  
//real[int] ElasticityBF = ElasticityForce; //ElasticityBoundary + ElasticityForce;  
real[int] ElasticityBF = ElasticityBoundary + ElasticityForce;  
real[int] uu0 = Elasticity^-1 \* ElasticityBF;  
Uh [uux0, uuy0, uuz0] = uu0;  
real coef=1.;

mesh3 Th0 = movemesh3(Th3T, transfo=[x+uux0_coef, y+uuy0_coef, z+uuz0\*coef]);  
mesh3 Th1;  
[uux0, uuy0, uuz0] = [uux0, uuy0, uuz0];  
//plot(Th0, [uux0, uuy0, uuz0],wait = true, cmm="coef amplification = "+coef);  
real dxmin0 = uux0[].min;  
real dymin0 = uuy0[].min;  
real dzmin0 = uuz0[].min;  
cout \<\< " - dep. max x0 = “\<\< dxmin0 \<\< " y0=” \<\< dymin0 \<\< endl;

// modified stiffness matrix for Newmark algorithm  
matrix Stiffness1 = Elasticity + alpha0_Mass + alpha1_Damping;  
matrix Stiffness2 = Elasticity + alpha2_Mass + alpha3_Damping;  
set (Stiffness1, solver=sparsesolver);  
set (Stiffness2, solver=sparsesolver);

// Prepare dynamic (transient problem)  
// computation of acceleration (UUtt)  
//temp = Damping_UUt;  
//temp += Elasticity_UU;  
//temp = ElasticityForce - temp;  
//UUtt = Mass^-1 \* temp;

//ofstream file(“f.txt”);

int fmmcyclen=0;  
real fmmtipvalueold;  
real fmmtipvalue;  
real fmmfreq;  
real fmmt0;  
real pmmforce,pmmdefax,pmmdefrad;  
real amplitude,phase;  
real[int] maxdef(24+2);

NSteps=0;  
amplitude=0;  
phase=0;  
cout \<\< “time,radial\_def,cera\_force,cera\_disp,cera\_speed,cera\_power,cera\_work” \<\< endl;  
Outputfile \<\< “time,radial\_def,cera\_force,cera\_disp,cera\_speed,cera\_power,cera\_work,trasdtip\_disp”;  
real zi;  
int ii;  
for (zi=0;zi\<=TuboL;zi+=TuboL/24.)  
{  
Outputfile \<\< “,def(z=” \<\< zi \<\< “)”;  
maxdef(ii)=0;  
ii++;  
}  
Outputfile \<\< endl;

real cerawork=0;

for (real t=0;t\<T;t+=dt)  
{  
macro fforce(ttt) ((-sin(2_pi_freqstart\*(ttt+toffset))+sin(2_pi_freqend\*(ttt+toffset)))/(pi\*(freqend-freqstart)_2_(ttt+toffset+1e-9)))  
// Bethe algorithm  
varf vForceHalf ([ux,uy,uz], [vx, vy,vz]) = int2d(Th3T,804)( vz \* ForceOnTip \* fforce(t+dt/2)) ;  
varf vForceFull ([ux,uy,uz], [vx, vy,vz]) = int2d(Th3T,804)( vz \* ForceOnTip \* fforce(t+dt)) ;  
pmmforce=ForceOnTip \* fforce(t+dt);  
// computation of modified elasticity (stiffness) matrix  
ForceHalf = vForceHalf(0, Uh);   
temp = alpha0_UU;  
temp += alpha4_UUt;  
temp += 1._UUtt;  
ForceHalf += Mass \* temp;  
temp = alpha1_UU;  
temp += 1.\*UUt;  
ForceHalf += Damping \* temp;

```
UUold = UU;
UUtold = UUt;
UUttold = UUtt;
UUhalf = Stiffness1^-1 * ForceHalf;
UUthalf = UUhalf-UU;
UUthalf *=alpha1;
UUthalf -=UUt;
UUtthalf= UUthalf-UUt;
UUtthalf*=alpha1;
UUtthalf-=UUtt;

// computation of second modified elasticity (stiffness) matrix 
ForceFull = vForceFull(0, Uh); 	
temp = alpha5*UUhalf;
temp += alpha6*UU;
temp += alpha1*UUthalf;
temp += alpha7*UUt;
ForceFull += Mass * temp;
temp = alpha1*UUhalf;
temp += alpha7*UU;
ForceFull += Damping * temp;

UUfull = Stiffness2^-1 * ForceFull;
UUtfull = -alpha7*UU;
UUtfull -=alpha1*UUhalf;
UUtfull +=alpha3*UUfull;

UUttfull = -alpha7*UUt;
UUttfull -=alpha1*UUthalf;
UUttfull +=alpha3*UUtfull;

UU = UUfull;
UUt = UUtfull;
UUtt = UUttfull;

Uh [uux, uuy, uuz] = UU;
Uh [uutx, uuty, uutz] = UUt;

real trasdtipdisp=int2d(Th3T,lbTrasdEnd)((uux*x+uuy*y)/(x^2+y^2)^0.5);
real ceradisp=int2d(Th3T,lbCera)(uuz);
real ceraspeed=int2d(Th3T,lbCera)(uutz);
real cerapower=pmmforce*ceraspeed;
real tubodisp=int2d(Th3T,lbTubo)(abs(uux*x+uuy*y)/(x^2+y^2+1.e-5)^0.5);
cerawork+=cerapower*dt;

//work+=(pmmdefax-pmmdefaxold)*pmmforce;
//if (work>workmax) workmax=work;

cout << t+dt << "," << tubodisp << "," << pmmforce << "," << ceradisp << "," << ceraspeed << "," << cerapower << "," << cerawork << endl;
Outputfile << t+dt << "," << tubodisp << "," << pmmforce << "," << ceradisp << "," << ceraspeed << "," << cerapower << "," << cerawork << "," << trasdtipdisp;

ii=0;
for (zi=0;zi<=TuboL;zi+=TuboL/24.)
{
	Outputfile << "," << uux(TuboRe,0.,zi);
	if (uux(TuboRe,0.,zi)>maxdef(ii)) 
		maxdef(ii)=uux(TuboRe,0.,zi);
	if (maxdef(ii)>defradmax) defradmax=maxdef(ii);
	ii++;
}
Outputfile << endl;
if (1 && (NSteps % 20==0))
{
	//Movemesh 
	coef = 15000;
	Th1 = movemesh3(Th3T, transfo=[x+uux*coef, y+uuy*coef, z+uuz*coef]); 
	[uux, uuy, uuz] = [uux, uuy, uuz];
	plot(Th1, [uux, uuy, uuz], wait = false);
	int[int] fforder = [1];
	savevtk("Defo"+NSteps+".vtu", Th1, [uux, uuy, uuz], dataname = "UU", order = fforder, append = false);
	savemesh(Th1,"mesh"+NSteps+".mesh");
}
NSteps++;

```

}

So I struggled a lot but finally I found a version that worked well

---

<div class="post-metadata">

### Author: ![Markov](https://avatars.discourse-cdn.com/v4/letter/m/f1d935/32.png) [@Markov](https://community.freefem.org/u/Markov)
#### Post date: [July 8, 2021, 1:32pm UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863/5 "2021-07-08T13:32:48Z")

</div>

Dear Riccardo,

Thank you for your detailed reply. I was also trying to implement the Newmark algorithm from Bathe’s book but I was having wrong results (maybe because on how I implemented it but I will double check). I am looking right now at the Bathe method in the book and the code that you sent me. I would like to ask you a question about the static calculation at the beginning, why exactly did you calculate it, is it a specific constraint for the problem?

Mark

---

<div class="post-metadata">

### Author: ![HGH](https://avatars.discourse-cdn.com/v4/letter/h/3be4f8/32.png) [@HGH](https://community.freefem.org/u/HGH)
#### Post date: [July 8, 2021, 1:38pm UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863/6 "2021-07-08T13:38:27Z")

</div>

Hello guys,

I also interested in Newmark-Beta method. I have completed the Newmark-Beta program by myself, and I am very happy to see that your program is also successful. I would like to ask if you know how to parallelize the Newmark-Beta program in FreeFem++?

I have encountered a big problem now, that is, when calculating 250,000 elements, the memory needs 160GB. Although I use a workstation to calculate, it is still too big and it takes a long time. Do you have a good way?

Best  
Jack Huang

---

<div class="post-metadata">

### Author: ![Markov](https://avatars.discourse-cdn.com/v4/letter/m/f1d935/32.png) [@Markov](https://community.freefem.org/u/Markov)
#### Post date: [July 8, 2021, 2:26pm UTC](https://community.freefem.org/t/transient-problems-in-elasticity/863/7 "2021-07-08T14:26:33Z")

</div>

Hello Jack,  
I have not been able to implement the Newmark Method yet, Riccardo was able to implement the Bathe method as you can see above. However, I will also try to parallelize the Newmark-Beta program on a workstation once I figure out the mistakes. I’ll let you know if I will be able to make it. On a side note do you have any references for the Newmark-Beta algorithm that you implemented?

Best  
Mark