Laboratoire RDP » Biophysics Team » Dandelion pappus morphing

include "getARGV.idp"

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);

}

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

{ /*  Computes the angle between the MN vector and Oy (in degrees) */

        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);

// hydrophobicity as swelling ability

real sm = -getARGV("-sm", 0.568354);// mesophyl swelling property

real fsn = getARGV("-fsn", 0.899975);// nectary     sn/sm

real fss = getARGV("-fss", 1.19112);// side        ss/sm

real fsv = getARGV("-fsv", 0.491088);// vasculature sv/sm

string outputfile = getARGV("-outfile", "test.csv");

// ---------------- 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;     //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=20; 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 

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

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

 *                  GEOMETRY and MESH

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

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

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

Sh0 nectaryh=nectary;

Sh0 sideh=side;

Sh0 mesophylh=mesophyl;

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;

func K=lambda+mu;

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

 *                  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) 

  ;

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

 *         Measuring angles

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

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);

// compute sideangle

height=hside;

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

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

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

 *         Measuring original areas

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

// original volumes

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);

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

 *         Measuring strain anisotropies

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

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;

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

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

 *         Measuring area changes

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

// deformed volumes

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

redefineVariable(vasculatureh);

redefineVariable(nectaryh);

redefineVariable(sideh);

redefineVariable(mesophylh);

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);

real am=volmesophyl0/volmesophyl;

real an=volnectary0/volnectary;

real as=volside0/volside;

real av=volvasculature0/volvasculature;

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

 *         Writing results

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

cout<<"sm;fsn;fss;fsv;sideangle;am;an;as;av;sam;san;sas;sav;"<<endl;

string toprint=string(sm)+";"+string(fsn)+";"+string(fss)+";"+string(fsv)+";"

+string(sideangle/38)+";"

+string(am)+";"+string(an)+";"+string(as)+";"+string(av)+";"

+string(samesophyl)+";"+string(sanectary)+ ";"+string(saside)+";"+string(savasculature);

cout<<toprint<<endl;

ofstream textfile(outputfile);

textfile<<toprint<<endl;