ElastoHydro.edp

//This code computes the solution of the elasto-hydrodynamic problem for myxobacteria gliding 
//Author: Joel Tchoufag (joel.tchoufag@gmail.com)

verbosity=0;
load "medit"
load "iovtk" // for ParaView postprocessing
load "lapack"
load "fflapack"
load "msh3"
load "UMFPACK64"
defaulttoUMFPACK64();
						
real tgv=1e30;
real yc=0;

///////////////////// Mesh setup ///////////////////////////////////
real lim=2.5;
int n,k;
real L=2.*lim;//Length of the bacterial membrane
border a1(t=0,L){x=t;y=0;label=1;}
border a2(t=0,0.02){x=L;y=t;label=2;}
border a3(t=L,0){x=t;y=0.02; label=3;}//patch
border a4(t=0.02,0){x=0;y=t;label=4;}

//for (n=10;n<=10;n=n+1) //the convergence was done for n=4,8,...,24 with a reference at n=28
//{
n=6;
mesh Th = square(n*6*2.*lim,2*pi*1./lim*lim/5,[0+2*lim*x,y/60]);  

//plot (Th, wait=1, value=1);
//////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////// Finite element space////////////////////////////////
fespace Vh(Th,P1,periodic=[[1,x],[3,x]]);
////////////////////////////////////////////////////////////////////////

///////////////////// Definition of fields ///////////////////////////////
Vh pn,ppn,p0,dpn; // pressure 
Vh dtn,dt0, ddtn, qn; // substrate deformation
Vh hn, qnn,qnn2; // membrane height

//////// other variables ////
Vh dHn, dhn, Hn;
Vh xp;

xp=0;
real lambda;

int nuse=pn.n ;
cout << " N DOF " << nuse << endl;
cout << " NV " << Th.nt << endl;

Vh u=0, v=0; 

real ng=1;
int ncel=1;

real Dpn, dDpn; // for pressure boundary condition
real Dt=0;

/////////////// Parameters ////////////////////

real eta; // softness number
real amp=0.3; //membrane height
int  nbp=20; //number of time steps over a period
real nnp=1.; 
real per=1./nnp; // number of time periods over which to run the simulation
real dg=per/(nbp+1); //time step
real eps=0.01; //lubrication parameter

real wavespeed=0.32*1e-6; // wave/motors speed

//real Cas=15/0.32*wavespeed/(60e-3*0.5)/(eps^3); //Warning: this is Ca_s/epsilon^3 (see theory)!

//real alpha=1e-9*pi/1e-6; // slime-air interface thickness

//real tref=0.16;// surface tension ratio

//real Ca=Cas/(tref)*eps^3; // slime-air capillary number

//////////////////////////////////////////

int nbpos=1500;

	Vh delta,ddelta, delta1, deltaf, delstar, dpx, dpnstar, pstar, shear, Thrust,gp;
	
	real gn,dgn; // Lagrange multiplier to enforce zero lift
	
	real Un,Unn,Uav, dUn,alf,dalf, Dpav,Gnav, ellav, Gap;
	
	real Dragtemp, Thrustemp; 
	
	real Torqav, Powav, Thrustav, Pout; // time average quantities
	
	real[int] Unarray(nbp+1), Torqarray(nbp+1), Powarray(nbp+1), Poutarray(nbp+1), Dpnarray(nbp+1);
	real[int] Gnarray(nbp+1), ellarray(nbp+1);
	
	
	real[int] pnav(nbpos), pstarav(nbpos), shearav(nbpos), deltav(nbpos), xpos(nbpos), propulsav(nbpos),gpav(nbpos);
	
	//ofstream Vlout("Slope0p25Tref0p1Eta1.txt");
	
delta=0;
Un=-1;
pn=0;
real xc=lim;

Dpn=0;

real gam=0;
/////////////////// Sinusoidal membrane ///////////////////////
Vh hn0=amp*sin(2*pi*nnp*(x+gam)); // shape b(x)

Vh dho=amp*(2*pi*nnp)*cos(2*pi*nnp*(x+gam)); // slope b'(x)

Vh h1=hn-1.;
Vh htot=1.+hn0; // total gap h(x)=1+b(x)
/////////////////////////////////////////////////////////////

//////////////////// Weak form for the initialization (arbitrary) //////////////////
varf ReyA(pn,qn) = int1d(Th,1)(dx(pn)*qn*(1+hn0)^3)	+ on(2,pn=Dpn);
varf ReyA1(pn,qn) = int1d(Th,1)(-1.*qn);
varf ReyB1(pn,qn) = int1d(Th,1)(6*Un*(x)*qn);//int1d(Th,1)(x*qn);		
varf ReyE(pn,qn) = int1d(Th,1)(1.*qn) ;
varf ReyG(pn,qn) = int1d(Th,1)(x*qn);

varf Reyb(pn,qn)=int1d(Th,1)(-6*Un*(1+hn0)*qn)+ on(2,pn=Dpn)	;

matrix Ag=ReyA(Vh,Vh);
real[int] Ag1=ReyA1(0,Vh);
real[int] Bg1=ReyB1(0,Vh);
real[int] Eg=ReyE(0,Vh);
real[int] Gg=ReyG(0,Vh);

matrix A0=[[Ag,Ag1],
		[Eg',1e-12]];//'
																			
set(A0,solver=sparsesolver);											
real[int] b0(nuse+1),sol0(nuse+1);

b0(0:nuse-1)=Reyb(0,Vh);
b0(nuse)=0;
//b0(nuse+1)=0;

//sol0=(A0^-1)*b0;
	
pn[]=sol0(0:nuse-1);
delta=0;
ddelta=1.;
gn=0;//gn=sol0(nuse);

Vh pn0=pn;
Vh Un0=Un;	
//////////////////////// end initialization ////////////////////////////////////////
	
cout << "initial multipliers " << gn << " " << alf << endl;		
cout << "check zero lift constraint " << int1d(Th,1)(pn) << endl;	
plot(hn0,wait=0,dim=3,fill=1,value=1,cmm="initial membrane");
plot(pn,wait=0,dim=3,fill=1,value=1,cmm="initial pressure");

real X0=1.;

//////////////////////////////Green's function for substrate deformation////////////////
macro K(xx) log((xx)) //

func real CIJ(real xx,real yy, real eta0)
	 			{ 
				real s=0;
				s=int1d(Th,1)( K(abs(xx-x))*u/(pi));  // for elastohydrodynamic
				return s;
				}
///////////////////////////////////////////////////////////////////////////////////
		
ofstream Aout("Meanvalues.txt"); //output for instantaneous scalars: velocities, lagrange multiplier, ...

real etamin=0.01; // min softness parameter  (change if necessary)
real etamax=0.09; // max softness parameter (change if necessary)
real deta=0.01;// parameter increment step for Newton's method (can be changed cautiously)

//real Camin=2E-4/eps^3;
//real Camax=9E-2/eps^3;

int etacount=0;

//eta=500;
ofstream Vloutav("SlopeMean.txt"); // 

// One can loop over different capillary number, or different softness parameter, etc.

for  (eta=etamin; eta<=etamax; eta=eta+deta)//(tref=tmin; tref<=tmax; tref=tref+deta)// (Cas=Camin; Cas<=Camax; Cas=Cas*2)//
{

	ofstream SToutav("ProfileMean"+etacount+".txt"); //output for time averaged fields
	etacount++;
/*	if (etacount==10)
	{
		deta=deta*1;
		etacount=0;
	};
*/
	
//alpha=0.01*eta;			 
real[int,int] Cm(nuse,nuse);
real[int,int] Fm(nuse,nuse), Opm(nuse,nuse);
real[int] Lm(nuse);

/////////////////////////// construct shape function for the non-local substrate field /////////////////

for (int i=0;i<nuse;i++)    //	all degree of freedom
{
u[][i]=1;                       //	the basic function i
qnn=CIJ(x,y,eta);							  
Lm(i)=int1d(Th,1)(u);
for (int j=0;j<nuse;j++)    //	all degree of freedom
 { 
  v[][j]=1; 
  Cm(i,j)=int1d(Th,1)(-v*qnn);
  //Fm(i,j)=-1*Cm(i,j);//int1d(Th,1)(v*u);
  v[][j]=0;
 }
u[][i]=0;               //	reset
};

{
ofstream cmout("Cmatrix.txt");
cmout << Cm;
};
///////////////////////////////////////////////////////////////////////////


///////////////////////Temporal Loop////////////////////////////
gam=0;
hn0=amp*sin(2*pi*nnp*(x+gam)); //shape b(x)
dho=amp*(2*pi*nnp)*cos(2*pi*nnp*(x+gam)); // slope b'(x)
for (int na=0;na <=nbp; na=na+1)//(int nnp=4;nnp<=4;nnp++) //(int nnp=1;nnp<=5;nnp++)//
{
	
cout << " time step " << gam << endl;

func real dF(real xx,real yy,real eta0)
	 			{ 
				real s=0;
				//s=int1d(Th,1)( K(abs(xx-x)+0*2./pi*eta0/Cas)*pn/(pi)); // for elasto-capillary-hydrodynamic
				s=int1d(Th,1)( K(abs(xx-x))*pn/(pi) );  // for elastohydrodynamic                              
			    return s;
				}				
ncel=1;
	
	real xx;

	cout << "U and max p and delta " << Un << " " << pn[].max << " " << delta[].max << endl;

	real Dg0=(-1.*int1d(Th,1)(pn*dho+dx(pn)*(1.+hn0-eta*delta)/2.+1./(1.+hn0-eta*delta)));
	
	cout << "Drag " << Dg0 << " " << 1./(int1d(Th,1)(1.)) *(-1.*int1d(Th,1)(pn*dho+dx(pn)*(1.+hn0-eta*delta)/2.+1./(1.+hn0-eta*delta))) << endl;
	
	real L0;
	
	///////////////////// Variational/Weak forms for Newton jacobian ///////////////////////////////////////////////////////
	varf ReyJacobA(dpn,qn) = int1d(Th,1)(dx(dpn)*qn*(1.+hn0-eta*delta)^3 ) + on(2,dpn=Dpn-pn);
	
	varf ReyJacobB(dtn,qn) = int1d(Th,1)(-3*eta*dtn*dx(pn)*(1.+hn0-eta*delta)^2*qn+6*Un*(-eta*dtn)*qn);
	
	varf ReyJacobD(dtn,qn) = int1d(Th,1)(dtn*qn);
	
	varf ReyJacobA1(dpn,qn) = int1d(Th,1)(-1.*qn);
		
	varf ReyJacobB1(dpn,qn) = int1d(Th,1)(6*(1+hn0-eta*delta)*qn);	

	varf ReyJacobA2(dpn,qn) = int1d(Th,1)(0*qn)+int1d(Th,2)(-1*tgv*qn);	

	varf ReyJacobC1(dtn,qn) = int1d(Th,1)(1e-12*qn);
	
	varf ReyJacobD1(dtn,qn) = int1d(Th,1)(1e-12*qn);

	varf ReyJacobE(dpn,qn) = int1d(Th,1)(1.*qn) ;

	varf ReyJacobF(dtn,qn) = int1d(Th,1)(0*qn);

	varf ReyJacobG(dtn,qn) = int1d(Th,1)(1.*(dho*qn+dx(qn)*(1.+hn0-eta*delta)/2.));

	varf ReyJacobH(dtn,qn) = int1d(Th,1)(1.*(-eta*dx(pn)*qn/2.+eta*(Un+2.)*qn/(1.+hn0-eta*delta)^2));
		
	func real ReyJacobH1(real eta0)
		 			{ 
					  real s;
					  s=int1d(Th,1)(1./(1.+hn0-eta0*delta));                    
					  return s;
					}
	
	Vh delta2=dF(x,y,eta);
	ddelta=dx(delta2); 
	
	varf ReyJacobI(dpn,qn) = 
	int1d(Th,1)(0*qn);
	
	varf Fp(dpn,qn) = int1d(Th,1)(-dx(pn)*qn*(1.+hn0-eta*delta)^3-6*Un*(1.+hn0-eta*delta)*qn+gn*qn)
				+ on(2,dpn=Dpn-pn);
		
	varf Fdelta(dtn,qn)= int1d(Th,1)(-delta*qn+delta2*qn);

	func real Fg(real eta0)
		 			{ 
					real s;
					s=int1d(Th,1)(-pn);                  
					return s;
					}
	
	func real FU(real eta0)
		 			{ 
					real s;
					s=( -1.*int1d(Th,1)(pn*dho+dx(pn)*(1+hn0-eta0*delta)/2.+(Un+2.)/(1.+hn0-eta0*delta)) );               
				    return s;
					}
	
	func real FDp(real eta0)
		 			{ 
					real s;
					s=-Dpn ;
					return s;
					}
	///////////////////////////////////////////////////// End of weak forms ///////////////////////////////////////
/*ofstream DFout("DragAB_Mamp"+na+"C"+ncel+"c.txt");
ofstream LFout("LiftAB_Mamp"+na+"C"+ncel+"c.txt");
*/

	int niter;
	real Drag, Drag2, Lift, Torque, Torque2, Power;
	
	////////////////////////// Matrices to build stiffness matrix //////////////////////// 
	matrix A,B,C,Ce,D,Jn;
	
	int nuse2=2*nuse;

	real[int] sol(nuse2+3),dsol(nuse2+3),Fsol(nuse2+3);
	real[int] Lz(nuse);
	
	Lz=0*Lm;
		
	sol(0:nuse-1)=pn[];
	sol(nuse:2*nuse-1)=delta[];
	sol(2*nuse)=gn;
	sol(2*nuse+1)=Un;
	sol(2*nuse+2)=Dpn;
	
	Fsol(0:nuse-1)=Fp(0,Vh);
	Fsol(nuse:2*nuse-1)=Fdelta(0,Vh);
	Fsol(2*nuse)=Fg(eta);
	Fsol(2*nuse+1)=FU(eta);
	Fsol(2*nuse+2)=FDp(eta);
		
	{
	ofstream cmout("Dmatrix.txt");
	cmout << D;
	};
	
	C=Cm;

	//real CZ=1;
	//cout << "CZ=  " << CZ << endl;
	Vh dsp, dsd, dsH, dsh;
	real  des, de, dlam;
	Vh ptemp, deltemp;
	
	int count;
	
	real Unstar, Unstar2, gnstar, gtemp, Dragstar;
	real sgn;

cout << " max p, delta, gn, Dpn before " << pn[].max << " " << delta[].max <<  " " << gn << " " << Dpn <<  endl;
		
		niter=0;
		A=ReyJacobA(Vh,Vh);
		B=ReyJacobB(Vh,Vh);
		D=ReyJacobD(Vh,Vh);
		real[int] A1=ReyJacobA1(0,Vh);
		real[int] B1=ReyJacobB1(0,Vh);
		real[int] C1=ReyJacobC1(0,Vh);
		real[int] D1=ReyJacobD1(0,Vh);
		real[int] E=ReyJacobE(0,Vh);
		real[int] F=ReyJacobF(0,Vh);
		real[int] G=ReyJacobG(0,Vh);
		real[int] H=ReyJacobH(0,Vh);
		real[int] I=ReyJacobI(0,Vh);
		real[int] J=F; //=0
		real E1=1e-12;//ReyJacobE1(eta);
		real F1=0;
		real G1=0;
		real H1=ReyJacobH1(eta);
		real I1=0;
		real J1=0;
		real I2=1.;
		real[int] A2=ReyJacobA2(0,Vh);
		real[int] C2=F; //=0
		real E2=0;
		real G2=0;	
		
		//Ce=1./ng*C;
		Fsol(0:nuse-1)=Fp(0,Vh);
		Fsol(nuse:2*nuse-1)=Fdelta(0,Vh);
		Fsol(2*nuse)=Fg(eta);
		Fsol(2*nuse+1)=FU(eta);
		Fsol(2*nuse+2)=FDp(eta);
		
	//	cout << " A matrix " << A << endl;
	//	cout << " B matrix " << B << endl;
		
			Jn=[[A,B,A1,B1,A2],
				[C,D,C1,D1,C2],
				[E',F',E1,F1,E2],
				[G',H',G1,H1,G2],
				[I',J',I1,J1,I2]];
	
		set(Jn,solver=sparsesolver);
		dsol=Jn^-1*Fsol;
		sol+=dsol;
		
		dsp[]=dsol(0:nuse-1);
		dsd[]=dsol(nuse:2*nuse-1);
		dgn=dsol(2*nuse);
		dUn=dsol(2*nuse+1);
		dDpn=dsol(2*nuse+2);
		
		pn[]=sol(0:nuse-1);
		delta[]=sol(nuse:2*nuse-1);
		gn=sol(2*nuse);
		Un=sol(2*nuse+1);
		Dpn=sol(2*nuse+2);
		
		delta2=dF(x,y,eta);
		ddelta=dx(delta2);
	
		//cout << " max p, delta, gn, Un  " << pn(2*lim,0) << " " << delta(2*lim,0) <<  " " << gn << " " << Un <<  endl;
		cout << " pressure at leading edge " << pn(2*lim,0) << " " << Dpn << endl;
		
		while ( sqrt(int1d(Th,1)(dsp^2+dsd^2)/int1d(Th,1)(1.)+dgn^2+dUn^2+dDpn^2) >=1e-12)	
		{
			niter++;
			//cout << "in the loop " << endl;
	
			A=ReyJacobA(Vh,Vh);
			B=ReyJacobB(Vh,Vh);
			D=ReyJacobD(Vh,Vh);
			real[int] A1=ReyJacobA1(0,Vh);
			real[int] B1=ReyJacobB1(0,Vh);
			real[int] C1=ReyJacobC1(0,Vh);
			real[int] D1=ReyJacobD1(0,Vh);
			real[int] E=ReyJacobE(0,Vh);
			real[int] F=ReyJacobF(0,Vh);
			real[int] G=ReyJacobG(0,Vh);
			real[int] H=ReyJacobH(0,Vh);
			real[int] I=ReyJacobI(0,Vh);
			real[int] J=F; //=0
			real E1=1e-12;//ReyJacobE1(eta);
			real F1=0;
			real G1=0;
			real H1=ReyJacobH1(eta);
			real I1=0;
			real J1=0;
			real I2=1.;
			real[int] A2=ReyJacobA2(0,Vh);
			real[int] C2=F; //=0
			real E2=0;
			real G2=0;			
			
			Fsol(0:nuse-1)=Fp(0,Vh);
			Fsol(nuse:2*nuse-1)=Fdelta(0,Vh);
			Fsol(2*nuse)=Fg(eta);
			Fsol(2*nuse+1)=FU(eta);
			Fsol(2*nuse+2)=FDp(eta);
			
				
			Jn=[[A,B,A1,B1,A2],
				[C,D,C1,D1,C2],
				[E',F',E1,F1,E2],
				[G',H',G1,H1,G2],
				[I',J',I1,J1,I2]];
		
			set(Jn,solver=sparsesolver);
			dsol=Jn^-1*Fsol;
			sol+=dsol;
			
			dsp[]=dsol(0:nuse-1);
			dsd[]=dsol(nuse:2*nuse-1);
			dgn=dsol(2*nuse);
			dUn=dsol(2*nuse+1);
			dDpn=dsol(2*nuse+2);
			
			pn[]=sol(0:nuse-1);
			delta[]=sol(nuse:2*nuse-1);
			gn=sol(2*nuse);
			Un=sol(2*nuse+1);
			Dpn=sol(2*nuse+2);
			
			delta2=dF(x,y,eta);
			ddelta=dx(delta2);
			
			if (niter>=15)
			{
			cout << "blow-up time = " << gam << endl;
			}
			
			cout << " Gap at the leading edge " << 1.+hn0(2.*lim,0)-eta*delta(2.*lim,0) << endl;
	
		//plot(pn,wait=0,dim=3,fill=1,value=1,cmm="pressure");
		
		Drag=(-1.*int1d(Th,1)(pn*dho+dx(pn)*(1.+hn0-eta*delta)/2.+(Un+2.)/(1+hn0-eta*delta))); //
		pstar= -1.*pn*dho;
		solve shearprob(shear,qn)=int1d(Th,1)(shear*qn)+int1d(Th,1)(qn*(dx(pn)*(1.+hn0-eta*delta)/2.+(Un+2.)/(1+hn0-eta*delta)));
		gp=1+hn0-eta*delta;
		Thrust=-(Un+2.)/(1+hn0-eta*delta);
		Lift=(int1d(Th,1)(pn)); //
		}
		/////////////////////End loop //////////////////////////////////////	
		
		Unn=Un;
		ptemp=pn;
		deltemp=delta;
		gtemp=gn;
		Dragtemp=Drag;
		

			htot=1.+hn0;	
		//	dho=dx(htot);
		ptemp=eta*delta;//dx(pn)*(1+hn0-eta*delta)/2.;
		
		///////////// Visualize different fields ////////////////
		//	plot(pn,wait=1,dim=3,fill=1,value=1,cmm="pressure");
			plot(pn, delta,wait=0,dim=3,fill=1,value=1,cmm="pressure substrate");
		//	plot(htot,ptemp,wait=0,dim=3,fill=1,value=1,cmm="total gap and pressure field");
		//	plot(Hn,wait=0,dim=3,fill=1,value=1,cmm="curvature");
		//	plot(htot,wait=0,dim=3,fill=1,value=1,cmm="membrane");
		//	plot(dho, wait=0, dim=3,fill=1,value=1,cmm="slope membrane");
		//////////////////////////////////////////////////////////////////////////
		
		cout << " number of iterations = " << niter << " for eta =  " << eta << endl;
	

		
	//} // end Un iteration 
htot=1.+hn0;

/*if(!(na%5)) // output for average fields: pressure, substrate deformation, stresses, (For ParaView visualisation)
	{
	ofstream STout("MeshConv/Profile0p1Mesh"+n+"Eta"+na+".txt");
	for (real tx=0; tx<=2.*lim; tx=tx+2.*lim/(nbpos-1))
	{
		STout << tx << " " << dho(tx,yc) << " " << pn(tx,yc) << " " << delta(tx,yc) << " " << pstar(tx,yc) <<  " " << eta*delta(tx,yc) << " " << hn0(tx,yc) << " " << shear(tx,yc) << " " << shear(tx,yc)-Thrust(tx,yc) << endl;
	} 
	}
	*/
	int xcount=0;
	for (real tx=0; tx<=2.*lim; tx=tx+2.*lim/(nbpos-1))
	{
		pnav(xcount)=pnav(xcount)+pn(tx,yc);
		pstarav(xcount)=pstarav(xcount)+pstar(tx,yc);
		shearav(xcount)=shearav(xcount)+shear(tx,yc);
		propulsav(xcount)=propulsav(xcount)+Thrust(tx,yc);
		deltav(xcount)=deltav(xcount)+delta(tx,yc);
		gpav(xcount)=gpav(xcount)+gp(tx,yc);
		xpos(xcount)=tx;
		xcount++;
	}
 
xp=xp+(Un+2.)*dg; // bacterial center of mass position

////////// Output for paraview: speed, position, lagrange multipliers, tip pressure, lift, drag///////////////////////////////
//Vlout << Un+2 << " " << xp << " " << << gam << " " << Un*eta << " " << gn << " " << Dpn << " " << gp << " " << int1d(Th,1)(pstar/(2*lim)) << " " << int1d(Th,1)((shear-Thrust)/(2*lim)) << endl;
//
		func real In(int jj)
		{
			int1d(Th,1)(1./(1+hn0-eta*delta)^jj);
		}
		Unarray(na)=Un;
		Torqarray(na)=Torque;
		Powarray(na)=Power;
		//Thrustarray(na)=Thrust;
		Poutarray(na)=(Un+2.)^2;	
		Dpnarray(na)=Dpn;
		Gnarray(na)=2-(Dpn-0*pn(0,0))/(6.*In(2))*(1-1./(1-1.5*In(2)^2/(In(1)*In(3))))+(1./(1-1.5*In(2)^2/(In(1)*In(3))));
		// 2-Dpn*1./(6.*In(2))*(1-1./(1-1.5*In(2)^2/(In(1)*In(3))))+(1./(1-1.5*In(2)^2/(In(1)*In(3))));//3+(gn*(int1d(Th,1)(1./(1+hn0-eta*delta)^2)))*1./(4.*int1d(Th,1)(1./(1+hn0-eta*delta)^1));
		ellarray(na)=2-(Dpn-pn(0,0))/(6.*In(2))*(1-1./(1-1.5*In(2)^2/(In(1)*In(3))))+(1./(1-1.5*In(2)^2/(In(1)*In(3)))); //eps/Ca;
	//

gam+=dg;
hn0=amp*sin(2*pi*nnp*(x+gam)); //(cos(nnp*x)+sin(nnp*x));
dho=amp*(2*pi*nnp)*cos(2*pi*nnp*(x+gam)); //(-sin(nnp*x)+cos(nnp*x));


}
///////////////////////////////////// end temporal loop///////////////////////

////////////////////////////////////Averaging post-process//////////////////////////////
	func real avr(real[int] ar)
	{
		real uavr=0;
		for (int k=0; k<=nbp; k++)
		{
			uavr=uavr+ar(k)*dg/per;
		}
		return uavr;
	}
	
		Uav=avr(Unarray);
		Torqav=avr(Torqarray);
		Powav=avr(Powarray);
		Pout=avr(Poutarray);
		Dpav=avr(Dpnarray);
		Gnav=avr(Gnarray);
		ellav=avr(ellarray);
		
	func real avr2(real[int] ar)
	{
		real uavr=0;
		for (int k=0; k<=nbpos-1; k++)
		{
			ar=dg/per*ar;
			uavr=uavr+ar(k)*1./nbpos;
		}
		return uavr;
	}
		Gap=avr2(gpav);
		
	
Aout << eta << " " << Uav << " " << ellav << " " << Torqav << " " << Gap << " " << Gnav << " " << Dpav << endl;
	
	cout << nbpos << " " << "size pnav = " << deltav.n << endl;
	
	
	for (int kn=0; kn<=nbpos-1; kn++)
		{
		SToutav << xpos(kn) << " " << pnav(kn)*dg/per << " " <<  pstarav(kn)*dg/per << " " << shearav(kn)*dg/per << " " << propulsav(kn)*dg/per << " " << deltav(kn)*dg/per << endl;
		}
////////////////////////////End averaging and saving //////////////////////////////////////////

}
cout << "gliding speed is " << Uav << endl;

//}// end n