Laboratoire RDP » Biophysics Team » Dandelion pappus morphing

include "getARGV.idp"

//include "bib_meca2d.cpp"

//usage :

//Freefem++ pulvinus.edp [-fE fEvalue] [-Nit Nit] [-geometry g]

//arguments:

//-fE fEvalue:

//-Nit Nit: number of iterations

//-geometry g: differnet geometries

//1: 

//2: 

//3: 

func real lineD(real xA, real xB, real yA, real yB)

{        return (yB-yA)/(xB-xA);

}

func real lineC(real xA, real xB, real yA, real yB)

{        return (yA*xB-yB*xA)/(xB-xA);

}

//real PI=3.14159265;

func real angleToVertical(real Mx, real My, real Nx, real Ny)

{ /*  Computes the sinus of angle between the MN vector and Oy   */

        return asin((Nx-Mx)/sqrt((Mx-Nx)^2+(My-Ny)^2))*180/pi;

}

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

 *                  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 fE = getARGV("-fE", 0.2);

real sm = -getARGV("-sm", 0.95);// measured value

real fs = getARGV("-fs", 0.02);

int Nit = getARGV("-Nit", 1);

string geom = getARGV("-geometry", "1");

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

int hyp=getARGV("-hyp", 1);

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

cout<<"Hypothesis no. "<<hyp<< ":   sm="<<sm<< endl;

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

// (lengths are measured in micrometers)

real units=200;

// pulvinus dimensions

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

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

// nectary height

real hnect = 38.8/units;   //sd=9.6

real R=363.0/units;    //sd=49.4

// side dimensions

real hside=91.4/units;     //12.6

real lside=47.6/units/2;     //6.6

// vasculature dimensions

real vthickness=34.3/units;//5.5

real cwidth=74.8/units;    //15.2

real cwidthdry=44.2/units; //8.4

// vasculature displacement (from vascular bundle distance)

real vd=(cwidth-cwidthdry)/2;

real vposition=lx/2-vthickness-hnect/ly*(lx/2-cwidth/2-vthickness); //cout<<"vposition="<<vposition<<endl;

// mesh parameters

int nvertex=50; cout << "nvertex="<<nvertex<<endl;

real pinned=lx/nvertex/2;

// ---------------- Material properties

real nu = 0.29;    // Poisson's ratio

// Young modulus

real Ev;   // vasculature Young modulus

real fEm; // mesophyl Em/Ev 

real fEn;  // nectary  En/Ev 

real fEs;  // side     Es/Ev 

// hydrophobicity as swelling ability

//   sm        // mesophyl swelling property

real fsn;   // nectary     sn/sm

real fss;   // side        ss/sm

real fsv;   // vasculature sv/sm

// measured values for the density

//Region,Mean_tissue_density,SD_tissue_density

real dm=0.668060839; //sd=0.04091249

real dn=1.01896008; //sd=0.015464575

real ds=0.918032482; //sd=0.075097509

real dv=0.882171532; //sd=0.066651037

real nEd=1.0;

// Young modulus

Ev=1;   // vasculature Young modulus

fEm=(dm/dv)^nEd; // mesophyl Em/Ev 

fEn=(dn/dv)^nEd;  // nectary  En/Ev 

fEs=(ds/dv)^nEd;  // side     Es/Ev 

// hydrophobicity as swelling ability

cout << "hyp="<<hyp<<"============================="<<endl;

if (hyp==1) 

{

fsn=fs;   // nectary     sn/sm

fss=fs;   // side        ss/sm

fsv=fs;   // vasculature sv/sm

}

if (hyp==2) 

{

fsn=fs;   // nectary     sn/sm

fss=1;   // side        ss/sm

fsv=fs;   // vasculature sv/sm

}

if (hyp==3) 

{

fsn=1;   // nectary     sn/sm

fss=1;   // side        ss/sm

fsv=fs;   // vasculature sv/sm

}

if (hyp==4) 

{

fsn=1;   // nectary     sn/sm

fss=fs;   // side        ss/sm

fsv=fs;   // vasculature sv/sm

}

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

 *                  GEOMETRY and MESH

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

string geomfilename="geometry"+geom+".cpp";

cout<<"including "<<geomfilename<<endl;

include "geometry1.cpp";

// -------------------- 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 strain, stress;

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

func bb=[[-lx/2*1.5,-ly*1],[lx/2*1.5,ly*.1]];

real coef=1;

cout << "Coefficent of amplification:"<<coef<<endl;

//plot(Th, fill=1, ps=output+"/original_geometry.png", bb=bb); -------------------------------

// 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

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

 *                  MECHANICS

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

// definition of integer functions for the structural domains

func vasculature=int(vasculatureBool(x,y));

func nectary=int(nectaryBool(x,y));

func side=int(sideBool(x,y));

func mesophyl=int(mesophylBool(x,y));

func geometry = 1*vasculature + 3*nectary + 2*side + 4*mesophyl;

// definition of FE variables for the structural domains

Sh0 vasculatureh=vasculature;

//plot(vasculatureh,wait=1,value=true,fill=1, ps=output+"/vasculature-original.png", bb=bb);

Sh0 nectaryh=nectary;

//plot(nectaryh,wait=1,value=true,fill=1, ps=output+"/nectary-original.png", bb=bb);

Sh0 sideh=side;

//plot(sideh,wait=1,value=true,fill=1, ps=output+"/side-original.png", bb=bb);

Sh0 mesophylh=mesophyl;

//plot(mesophylh,wait=1,value=true,fill=1, ps=output+"/mesophyl-original.png", bb=bb);

Sh0 geometryh=geometry;

string fname=output+"/geometry-original.png";

plot(geometryh,fill=1, ps=fname, bb=bb);

// spatial dependence of Young modulus

func E=Ev*(vasculature + fEn*nectary + fEs*side + fEm*mesophyl);

Sh0 Eh=E;

// spatial dependence of the swelling property

func s=sm*(fsv*vasculature + fsn*nectary + fss*side + mesophyl);

Sh0 sh=s;

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

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

func K=lambda+2*mu/3;

// ------ iterations

{

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

 *                  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(32,u1=vd,u2=0)

  + on(33,u1=-vd,u2=0) 

  ;

stress=2*K*(strain-s);

Sh0 e11=dx(u1)+1.;

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

Sh0 e22=dy(u2)+1.;

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;

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

 *                  VISUALISATION

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

real voltotal0=int2d(Th)(1);

real volvasculature0=int2d(Th)(vasculatureh);

real volnectary0=int2d(Th)(nectaryh);

real volmesophyl0=int2d(Th)(mesophylh);

real volside0=int2d(Th)(sideh);

cout<<"Original volume="<<voltotal0<<endl;

cout<<"Original vasculature volume="<<volvasculature0<<endl;

// compute mean strain per region

real straintotal=int2d(Th)(strain)/voltotal0;

real strainvasculature=int2d(Th)(strain*vasculatureh)/volvasculature0;

real strainnectary=int2d(Th)(strain*nectaryh)/volnectary0;

real strainside=int2d(Th)(strain*sideh)/volside0;

real strainmesophyl=int2d(Th)(strain*mesophylh)/volmesophyl0;

// compute mean strain anisotropy per region

real satotal=int2d(Th)(strainanisotropy)/voltotal0;

real savasculature=int2d(Th)(strainanisotropy*vasculatureh)/volvasculature0;

real sanectary=int2d(Th)(strainanisotropy*nectaryh)/volnectary0;

real saside=int2d(Th)(strainanisotropy*sideh)/volside0;

real samesophyl=int2d(Th)(strainanisotropy*mesophylh)/volmesophyl0;

cout<<"Mean strain and strain anisotropy per region :"<<endl;

cout<<"Mesophyll:"<<strainmesophyl<<"   "<<samesophyl<<endl;

cout<<"Nectary:"<<strainnectary<<"   "<<sanectary<<endl;

cout<<"Side:"<<strainside<<"   "<<saside<<endl;

cout<<"Vasculature:"<<strainvasculature<<"   "<<savasculature<<endl;

cout<<"Total:"<<straintotal<<"   "<<satotal<<endl;

mesh Th0=Th;

Th=movemesh(Th,[x+u1*coef,y+u2*coef]);

plot(Th, Th0, ps=output+"/geometry-deformed-superposed_hyp"+string(hyp)+"_sm"+string(sm)+"_fs"+string(fs)+".png",bb=bb);

redefineVariable(strain);

redefineVariable(stress);

plot(strain, fill=1,wait=1,value=true, ps="strain.png",bb=bb);

plot(stress, fill=1,wait=1,value=true, ps="stress.png",bb=bb);

redefineVariable(geometryh);

plot(geometryh, fill=1, ps=output+"/geometry-deformed_hyp"+string(hyp)+"_sm"+string(sm)+".png",bb=bb);

redefineVariable(vasculatureh);

redefineVariable(nectaryh);

redefineVariable(sideh);

redefineVariable(mesophylh);

// compute structure volumes

real voltotal=int2d(Th)(1);

real volvasculature=int2d(Th)(vasculatureh);

real volnectary=int2d(Th)(nectaryh);

real volside=int2d(Th)(sideh);

real volmesophyl=int2d(Th)(mesophylh);

cout<<"Deformed total volume="<<voltotal<<endl;

cout<<"Deformed vasculature volume="<<volvasculature<<endl;

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

 *         LOOKING FOR ANGLE AND NECKHEIGHT

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

// displaced upper corner

real Nx=lx/2+u1(lx/2.,0), Ny=u2(lx/2.,0);

// compute tangentangle

real height=ly/nvertex/100;

real Mx=lx/2+u1(lx/2.,-height), My=-height+u2(lx/2.,-height);

real tangentangle=angleToVertical(Mx, My, Nx, Ny);

cout<<"Tangent angle="<<tangentangle<<endl;

// compute sideangle

height=hside;

Mx=lx/2+u1(lx/2.,-height); My=-height+u2(lx/2.,-height);

real sideangle=angleToVertical(Mx, My, Nx, Ny);

cout<<"Side angle="<<sideangle<<endl;

real am=volmesophyl0/volmesophyl;

real an=volnectary0/volnectary;

real as=volside0/volside;

real av=volvasculature0/volvasculature;

ofstream textfile(output+"/results.csv", append);

/*

hyp;fs;tangentangle;sideangle;am;an;as;av;sam;san;sas;sav

*/

textfile<<hyp<<";"<<fs<<";"<<tangentangle/38<<";"<<sideangle/38

<<";"<<am<<";"<<an<<";"<<as<<";"<<av

<<";"<<samesophyl<<";"<<sanectary<<";"<<saside<<";"<<savasculature

<<endl;

}//end for iteration loop

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

 *         Writing results

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