/* Code for Newmark-beta method */ load "msh3" load "tetgen" include "cube.idp" load "lapack" load "UMFPACK64" load "iovtk" load "MUMPS" load "scotch" load "shell" string filename; int[int] Nxyz=[20,4,10]; real [int,int] Bxyz=[[0.,1.],[0.,0.04],[0.,0.1]]; int [int,int] Lxyz=[[1,2],[3,4],[5,6]]; //mesh3 Th=Cube(Nxyz,Bxyz,Lxyz); mesh3 Th; if(mpirank==0) { Th=Cube(Nxyz,Bxyz,Lxyz); //Th=buildmesh(fr1(5*n)+fr2(n)+fr3(5*n)+fr4(n)+fr5(-n*3)); int[int] nupart(Th.nt); nupart=0; if(mpisize>1) scotch(nupart, Th, mpisize); Th=change(Th,fregion= nupart[nuTriangle]); //savemesh(Th,fTh); plot(Th,wait=0); } broadcast(processor(0),Th); real E = 1000;//1000; real nu =0.3; real rho = 1.; real mu = E/(2*(1+nu)); real lambda = E*nu/((1+nu)*(1-2*nu)); // Fespace fespace Vh(Th,[P2,P2,P2]); Vh [u1,u2,u3],[uu1,uu2,uu3], [vv1,vv2,vv3], [v1,v2,v3]; fespace Uh(Th,P2); Uh force; Uh uu,vv,ww; Uh uud,vvd,wwd; 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 div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM //real time=clock(); real timeI=mpiWtime(); real time1=mpiWtime(); varf stiffness([u1,u2,u3],[v1,v2,v3])= int3d(Th,mpirank)( lambda*div(u1,u2,u3)*div(v1,v2,v3) +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) ) + on(1,u1=0,u2=0,u3=0) ; varf mass([u1,u2,u3],[v1,v2,v3])= int3d(Th,mpirank)( rho*(u1*v1 + u2*v2 + u3*v3) ) ; real c1=0.; real c2=0.; varf damping([u1,u2,u3],[v1,v2,v3])= int3d(Th,mpirank)( c1*rho*(u1*v1+u2*v2+u3*v3) +c2*lambda*div(u1,u2,u3)*div(v1,v2,v3) +c2*2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) ) ; matrix K= stiffness(Vh,Vh,solver=GMRES, tgv=1e50); //stiffness matrix if(mpirank==0) cout << "size of matrix " << K.n << " x " << K.m << " nn nzero coef = " << K.nbcoef << endl; time1=mpiWtime()-time1; real time2=mpiWtime(); matrix M= mass(Vh,Vh,solver=CG, eps=1e-12); // mass matrix time2=mpiWtime()-time2; real time3=mpiWtime(); matrix C=damping(Vh,Vh,solver=CG, eps=1e-12); //Damping coefficient matrix time3=mpiWtime()-time3; /******************* Newmark-Beta method****************************/ real dt = 0.08;//1e-5; // time step in seconds real tt = 4;//0.4096; // total time in seconds int nsbt = int(tt/dt); // Number of time-steps real alpha = 0.25; // gamma for the test case real beta = 0.5; // alpha for the test case cout<< " nsbt " << nsbt << endl; cout << " K.n "<< Vh.ndof << " size array :"<< K.n*(nsbt+1) << endl; int n = Vh.ndof; real[int] F(n); real[int] U(n), Ud(n), Udd(n); //Displacement, Velocity and Acceleration real[int] Up(n), Upd(n), Updd(n); //Previous value of Displacement, Velocity and Acceleration real[int] dummy(n),dummy4(n),dummy5(n),q(n),r(n),xP(n); //real-valued arrays for dummy vaiables real[int]Ekin(n), Eel(n), Etot(n), Edamp(n); /* Calculation of Initial Acceleration */ real time4=mpiWtime(); matrix a1 = (1/(alpha*dt*dt))*M + (beta/(alpha*dt))*C + K; time4=mpiWtime()-time4; real time5=mpiWtime(); set(a1, solver=sparsesolver,master=-1); time5=mpiWtime()-time5; real[int] U1(K.m/3),U2(K.m/3),U3(K.m/3); real[int] U1d(K.m/3),U2d(K.m/3),U3d(K.m/3); real coef = 1.; //Amplification co-efficeint to see the plots int[int] ref2 = [1, 0, 0, 2]; /* Initial Conditions : Displacement and Velocity */ //savevtk("beam.vtu",Th); U = 0; Ud = 0; /********************* Time loop start ***************************/ //ofstream file("linearNM.dat"); for (int nt=0; nt<=nsbt; nt+=1){ //nsbt = 50 real t= nt*dt; real tc = tt/5; //tc = 0.8 if(nt>0){ Up=U; // nt