Laboratoire RDP » Biophysics Team » Dandelion pappus morphing

include "getARGV.idp"

//usage :

//Freefem++ test-formulation.cpp [-fE fEvalue] [-Nit Nit] [-geometry g]

//arguments:

//-fE fEvalue:

//-Nit Nit: number of iterations

//-geometry g: differnet geometries

//1: 

//2: 

//3: 

/*************************************************************

 *                  PARAMETERS

 * **********************************************************/

// reading script arguments

for (int i=0;i<ARGV.n;++i)

{ cout << ARGV[i] << " ";}

cout<<endl;

verbosity = getARGV("-vv", 0);

int vdebug = getARGV("-d", 1);

real E = getARGV("-E", 0.2);

real s = -getARGV("-s", 0.4);

string outputdir = getARGV("-out", ".");

cout<<"----------------------------------------------"<< endl;

cout<<"E="<<E <<"   s="<<s<< endl;

// ---------------- geometrical parameters

// (lengths are measured in micrometers)

real units=1;

// pulvinus dimensions

real lx=1/units;       //sd =34.4

real ly= 1/units;      //sd=36.0

// mesh parameters

int nvertex=10;

cout << "nvertex="<<nvertex<<endl;

int[int] labs=[10,20,30,40]; // bottom, right, top, left

mesh Th=square(lx*nvertex,ly*nvertex, label=labs, region=0, [x*lx,y*ly]);

// bounding box for the plot (use the same for all images so that they can be superposed)

func bb=[[-1,-1],[lx+1,ly+1]];

real coef=1;

string root=outputdir+"/test-formulation";

// -------------------- define the finite element space

fespace Vh(Th,[P2,P2]);  // vector on the mesh (displacement vector)

Vh [u1,u2], [v1,v2];

fespace Sh(Th, P2);   // scalar on the mesh, P2 elements

fespace Sh0(Th,P0);   // scalar on the mesh, P0 elements

fespace Sh1(Th,P1);   // scalar on the mesh, P1 elements

Sh0 Eh=E;// Young modulus

Sh0 sh=s;// hydrophobicity as swelling ability

real nu = 0.29;    // Poisson's ratio

func mu=E/(2*(1+nu));

func lambda=E*nu/((1+nu)*(1-2*nu));

//func K=lambda+2*mu/3;

func K=lambda+mu;

// macro to redefine variables on the displaced mesh 

macro redefineVariable(vvv)

{        real[int] temp(vvv[].n);

        temp=vvv[]; vvv=0; vvv[]=temp;

        vvv=vvv;

}//

real sqrt2=sqrt(2.);

macro epsilon(u1,u2)  [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] // EOM

//the sqrt2 is because we want: epsilon(u1,u2)'* epsilon(v1,v2) $==  \epsilon(\bm{u}): \epsilon(\bm{v})$

macro div(u1,u2) ( dx(u1)+dy(u2) ) // EOM

/*************************************************************

 *                  SOLVING THE FEM

 * **********************************************************/

solve Lame([u1,u2],[v1,v2])=

  int2d(Th)(

            lambda*div(u1,u2)*div(v1,v2)         

            + 2.*mu*( epsilon(u1,u2)'*epsilon(v1,v2) ) 

              )

  - int2d(Th) ( K*sh*div(v1,v2))

  + on(10,u2=0, u1=0)

  //+ on(40,u1=0) 

  ;

/*************************************************************

 * **********************************************************/

Sh0 e11=dx(u1)+1.;

Sh0 e12=1/2.*(dx(u2) + dy(u1));

Sh0 e22=dy(u2)+1.;

Sh0 strain=e11+e22 ;

Sh0 Det=e11*e22-e12*e12;

Sh0 l1=abs(strain+sqrt(strain*strain-4*Det))/2.;

Sh0 l2=abs(strain-sqrt(strain*strain-4*Det))/2.;

Sh0 lmax=(l1-l2+abs(l1-l2))/2.+l2;

Sh0 lmin=(l1-l2-abs(l1-l2))/2.+l2;

Sh0 strainanisotropy=lmin/lmax;

plot(strainanisotropy, fill=1, value=true, ps="satrainanisotropy-s"+string(s)+".png",  bb=bb, wait=1);

mesh Thm1=movemesh(Th, [x+u1, y+u2]);

plot(Th, Thm1, ps=root+"_geometry-moved.png", wait=1, bb=bb);

// compute original and deformed volume

real vol0=int2d(Th)(1);

real vol=int2d(Thm1)(1);

// compute mean strain

real meanstrain=int2d(Th)(strain/2)/vol0;

// compute mean strain anisotropy

real straina0=int2d(Th)(strainanisotropy)/vol0;

real straina1=int2d(Thm1)(strainanisotropy)/vol;

real strainamixed=int2d(Thm1)(strainanisotropy)/vol0;

cout<<"Initial volume:"<<string(vol0)<<endl;

cout<<"Final volume:"<<string(vol)<<endl;

cout<<"A1/A0="<<vol/vol0<<endl;

cout<<"Mean strain:"<<string(meanstrain)<<endl;

cout<<"Strain anisotropy 0:"<<string(straina0)<<endl;

cout<<"Strain anisotropy 1:"<<string(straina1)<<endl;

cout<<"Strain anisotropy mixed:"<<string(strainamixed)<<endl;

plot(Th, Thm1, ps="deformation-s"+string(s)+".png",  bb=bb, wait=1);