Index: README.txt
===================================================================
--- README.txt	(revision 0)
+++ README.txt	(revision 22)
@@ -0,0 +1,79 @@
+----------------------------------------------------------------------
+			        - Dandelion -
+    Dandelion pappus morphing is actuated by radially patterned material swelling
+----------------------------------------------------------------------
+
+FreeFem++ simulation scripte used in the following publication :
+Madeleine Seale, Annamaria Kiss, Simone Bovio, Ignazio Maria Viola, Enrico
+Mastropaolo, Arezki Boudaoud, Naomi Nakayama, "Dandelion pappus morphing is actuated by radially patterned material swelling", doi: https://doi.org/10.1101/2021.08.23.457337
+----------------------------------------------------------------------
+
+Author :
+--------
+Annamaria Kiss (annamaria.kiss@ens-lyon.fr, Laboratoire RDP, ENS Lyon)
+
+Download
+--------
+svn --username $USER checkout http://forge.cbp.ens-lyon.fr/svn/dandelion
+
+Dependency
+--------
+The scripte is interpreted and executed by the FreeFem++ software:
+https://freefem.org/
+
+Usage
+-----
+FreeFemm++ dandelion.cpp [-param value]
+
+If no options are given, the reference model is simulated, namely all parametervalues are chosen to have their reference values. Parameter values can be specified by using the parameter's name as option. Accepted options are the following (see Table 1 and Table 2 of the publication for parameter description and reference values).
+- geometrical parameters: D, H, Hpod, R, Hside, Wside, Wvasculature, Dcavity
+- vasculature displacement: dvasc
+- swelling factors: scort, spod, sside, svasc
+
+
+Example 1 : the reference model
+================================
+FreeFemm++ dandelion.cpp
+
+
+nvertex=40
+  --  mesh:  Nb of Triangles =   8774, Nb of Vertices 4629
+  -- Solve : 
+          min -0.292037  max 0.292132
+----------------------
+Areachanges (Awet/Adry) :
+----------------------
+cortex -> 1.67524
+floral podium -> 1.61194
+side -> 1.79592
+vasculature -> 1.2898
+----------------------
+Generated side angle = 19.5215 degrees
+----------------------
+times: compile 0.008144s, execution 9.69662s,  mpirank:0
+ CodeAlloc : nb ptr  5686,  size :562720 mpirank: 0
+Ok: Normal End
+
+Example 2 : higher swelling abilities for the cortex generates higher side angle
+================================
+FreeFem++ dandelion.cpp -scort 0.7
+
+
+nvertex=40
+  --  mesh:  Nb of Triangles =   8774, Nb of Vertices 4629
+  -- Solve : 
+          min -0.404528  max 0.404761
+----------------------
+Areachanges (Awet/Adry) :
+----------------------
+cortex -> 2.21836
+floral podium -> 1.54942
+side -> 1.54601
+vasculature -> 1.25095
+----------------------
+Generated side angle = 31.9336 degrees
+----------------------
+times: compile 0.007776s, execution 10.0292s,  mpirank:0
+ CodeAlloc : nb ptr  5686,  size :562720 mpirank: 0
+Ok: Normal End
+
Index: dandelion.cpp
===================================================================
--- dandelion.cpp	(revision 0)
+++ dandelion.cpp	(revision 22)
@@ -0,0 +1,338 @@
+/***************************
+ * ----------------------------------------------------------------------
+ * Dandelion pappus morphing is actuated by radially patterned material swelling
+ * ----------------------------------------------------------------------
+ * 
+ * FreeFem++ simulation scripte used in the following publication :
+ * 
+ * Madeleine Seale, Annamaria Kiss, Simone Bovio, Ignazio Maria Viola, Enrico
+ * Mastropaolo, Arezki Boudaoud, Naomi Nakayama, 
+ * "Dandelion pappus morphing is actuated by radially patterned material swelling",
+ * doi: https://doi.org/10.1101/2021.08.23.457337
+ * 
+ * Author :
+ * Annamaria Kiss 
+ * (annamaria.kiss@ens-lyon.fr, Laboratoire RDP, ENS Lyon)
+ * 
+ * ***************************/
+
+include "getARGV.idp"
+
+// reading script arguments
+for (int i=0;i<ARGV.n;++i)
+{ cout << ARGV[i] << " ";}
+cout<<endl;
+
+/***************************
+*      GLOSSARY
+*   nectary = floral podium
+*   mesophyl = cortex
+/****************************/
+
+/*************************************************************
+ *                  PARAMETERS
+ * **********************************************************/
+
+string output = getARGV("-out", ".");
+
+// ---------------- geometrical parameters ----------------
+
+// (lengths are measured in micrometers)
+real units=200;
+
+// actuator dimensions
+real lx=getARGV("-D", 491.0)/units;       //sd =34.4
+real ly=getARGV("-H", 240.0)/units;      //sd=36.0
+
+// floral podium
+real hnect = getARGV("-Hpod", 38.8)/units;   //sd=9.6
+real Rmeasured=getARGV("-R", 363.0)/units;    //sd=49.4
+real R=Rmeasured-hnect;
+
+// side
+real hside=getARGV("-Hside", 91.4)/units;     //12.6
+real lside=getARGV("-Wside", 47.6)/units;     //6.6
+
+// vasculature
+real vthickness=getARGV("-Wvasculature", 34.3)/units;//5.5
+real cwidth=getARGV("-Dcavity", 74.8)/units;    //15.2
+
+// vasculature displacement (from vascular bundle distance)
+real vd=getARGV("-dvasc", 15.3)/units; //15.3*0.22
+real vposition=lx/2-vthickness-hnect/ly*(lx/2-cwidth/2-vthickness);
+
+// mesh parameters
+int nvertex=40;
+cout << "nvertex="<<nvertex<<endl;
+
+// ---------------- Material properties ----------------
+
+// Poisson's ratio
+real nu = 0.29;
+
+// Young modulus
+real Ev;   // vasculature Young modulus
+real fEm; // cortex Em/Ev 
+real fEn;  // floral podium En/Ev 
+real fEs;  // side     Es/Ev 
+	
+// measured values for the density (Mean, SD)
+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
+
+// Young modulus is linked to measured densities
+Ev=1;   // vasculature Young modulus
+fEm=dm/dv; // cortex Em/Ev 
+fEn=dn/dv;  // floral podium  En/Ev 
+fEs=ds/dv;  // side     Es/Ev 
+
+// -----   swelling ability of each region is a parameter
+real ani = getARGV("-ani", 1.);// anisotropy of swelling property
+real sm = -getARGV("-scort", 0.464888);// cortex swelling property
+real sn = -getARGV("-spod", 0.443697);// floral podium swelling property
+real ss = -getARGV("-sside", 0.573482);// side swelling property
+real sv = -getARGV("-svasc", 0.236087);// vasculature swelling property
+
+
+/*************************************************************
+ *                  GEOMETRY and MESH
+ * **********************************************************/
+
+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 sinus of angle between the MN vector and Oy   */
+	return asin((Nx-Mx)/sqrt((Mx-Nx)^2+(My-Ny)^2))*180/pi;
+}
+
+// ----------------- External boundary/ -----------------
+
+real rnect=lx/2-vthickness;
+real tB=2*pi-acos(rnect/R);
+real tA=pi+acos(rnect/R);
+real xC=0;
+real yC=sqrt(R*R-rnect*rnect);
+
+// top
+border a1(t=0,1){x=lx/2+t*(lx/2-vthickness-lx/2);y=0;label=11;};
+border a2(t=tB,tA){x=xC+R*cos(t);y=yC+R*sin(t);label=12;};
+border a3(t=0,1){x=-lx/2+vthickness+t*(-lx/2-(-lx/2+vthickness));y=0;label=13;};
+// left
+border b1(t=0,1){x=-lx/2;y=t*(-ly);label=2;};
+// bottom
+border c1(t=0,1){x=-lx/2+t*(-cwidth/2-vthickness+lx/2);y=-ly;label=31;};
+border c2(t=0,1){x=-cwidth/2-vthickness+t*vthickness;y=-ly;label=32;};
+border c3(t=0,1){x=cwidth/2+t*vthickness;y=-ly;label=33;};
+border c4(t=0,1){x=cwidth/2+vthickness+t*(lx/2-cwidth/2-vthickness);y=-ly;label=34;};
+// right
+border d1(t=0,1){x=lx/2;y=-ly+t*(ly);label=4;};
+
+// -------------------- Internal boundary/ -----------------
+
+real xA=-rnect; real yA=0;
+real xB=-cwidth/2; real yB=-ly;
+real D=lineD(xA, xB, yA, yB);
+real C=lineC(xA, xB, yA, yB);
+real delta=D^2*(C-yC)^2-(1+D^2)*((C-yC)^2-(R+hnect)^2);
+real xM=(-D*(C-yC)+sqrt(delta))/(D^2+1);
+real yM=D*xM+C;
+real tN=2*pi+asin(-sqrt((R+hnect)^2-xM^2)/(R+hnect));
+real tM=3*pi-tN;
+
+// hole
+border vide1(t=0,1){x=-cwidth/2+t*(xM+cwidth/2);y=-ly+t*(yM+ly);label=51;};
+border vide2(t=tM,tN){x=xC+(R+hnect)*cos(t);y=yC+(R+hnect)*sin(t);label=52;};
+border vide3(t=0,1){x=-xM+t*(cwidth/2+xM);y=yM+t*(-ly-yM);label=53;};
+
+// ------------------- Define the mesh -------------------
+
+mesh Th = buildmesh (a1(nvertex*vthickness)+a2(nvertex*2*rnect)+a3(nvertex*vthickness) // top
+                     + b1(nvertex*ly) + d1(nvertex*ly) // sides
+                     + c1(nvertex*lx/2) + c2(nvertex*vthickness) // left bottom
+                     + vide1(nvertex*ly) + vide2(nvertex*(lx-2*vthickness)) + vide3(nvertex*ly)  //hole
+                     + c3(nvertex*vthickness) + c4(nvertex*lx/2)  // right bottom
+                     );
+
+/*************************************************************
+ *                  DEFINITION OF the REGIONS
+ * **********************************************************/
+
+real xA1=-(vposition+vthickness), yA1=-hnect;
+real xB1=-cwidth/2-vthickness, yB1=-ly;
+real xA2=(vposition+vthickness), yA2=-hnect;
+real xB2=cwidth/2+vthickness, yB2=-ly;
+real D1=lineD(xA1, xB1, yA1, yB1);
+real C1=lineC(xA1, xB1, yA1, yB1);
+real D2=lineD(xA2, xB2, yA2, yB2);
+real C2=lineC(xA2, xB2, yA2, yB2);
+
+func bool vasculatureBool(real xx, real yy)
+{	return (((yy>D1*xx+C1) && (yy<D1*(xx-vthickness)+C1)) || ((yy>D2*xx+C2)&& (yy<D2*(xx+vthickness)+C2)) );
+}
+
+func bool nectaryBool(real xx, real yy)
+{	return (yy>D1*(xx-vthickness)+C1)&&(yy>D2*(xx+vthickness)+C2);
+}
+
+func bool sideBool(real xx, real yy)
+{	return   (yy>-hside) && (((xx<-lx/2+lside)&&(yy<D1*xx+C1)) || ( (xx>lx/2-lside)&& (yy<D2*xx+C2))  );
+}
+
+func bool mesophylBool(real xx, real yy)
+{	return (!vasculatureBool(xx,yy)) && (!nectaryBool(xx,yy)) && (!sideBool(xx,yy));
+}
+
+
+/*************************************************************
+ *                  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, ps=output+"/mesh-original.eps", bb=bb, wait=1);
+
+// 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
+
+/*************************************************************
+ *                 SPATIAL DEPENDENCE OF MATERIAL PROPERTIES
+ * **********************************************************/
+
+// define region geometry as a finite element field
+
+//string mutname="mnsv";
+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.eps";
+plot(geometryh,fill=1, ps=fname, bb=bb, wait=1 );
+
+
+// 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=sv*vasculature + sn*nectary + ss*side + sm*mesophyl;
+Sh0 sh=s;
+
+// spatial dependence of the anisotropy of the swelling property
+Sh0 anih=ani;
+
+func mu=E/(2*(1+nu));
+func lambda=E*nu/((1+nu)*(1-2*nu));
+
+/*************************************************************
+ *                  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) ( sh/(1.+anih)*((2.*mu+lambda)*dx(v1)+lambda*dy(v2)) )
+  - int2d(Th) ( ani*sh/(1.+anih)*(lambda*dx(v1)+(2.*mu+lambda)*dy(v2)) )
++ on(32,u1=vd,u2=0) + on(33,u1=-vd,u2=0) 
+  ;
+
+
+/*************************************************************
+ *                  RESULTS
+ * **********************************************************/
+
+// compute original region 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);
+
+
+// move the mesh and plot the deformed mesh
+mesh Th0=Th;
+Th=movemesh(Th,[x+u1,y+u2]);
+
+plot(Th, Th0,  ps=output+"/mesh-original+deformed-superposed.eps",bb=bb, wait=1);
+plot(Th,  ps=output+"/mesh-deformed.eps",bb=bb, wait=1);
+
+redefineVariable(geometryh);
+plot(geometryh, fill=1, ps=output+"/geometry-deformed.eps",bb=bb, wait=1);
+
+// compute deformed region volumes
+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);
+
+// compute areachanges
+real am=volmesophyl0/volmesophyl;
+real an=volnectary0/volnectary;
+real as=volside0/volside;
+real av=volvasculature0/volvasculature;
+
+// compute the generated side angle
+real height=hside;
+real Nx=lx/2+u1(lx/2.,0), Ny=u2(lx/2.,0);
+real Mx=lx/2+u1(lx/2.,-height), My=-height+u2(lx/2.,-height);
+real sideangle=angleToVertical(Mx, My, Nx, Ny);
+
+/*************************************************************
+ *        PRINT RESULTS
+ * **********************************************************/
+cout<<"----------------------"<<endl;
+cout<<"Areachanges (Awet/Adry) :"<<endl;
+cout<<"----------------------"<<endl;
+cout<<"cortex -> "<<am<<endl;
+cout<<"floral podium -> "<<an<<endl;
+cout<<"side -> "<<as<<endl;
+cout<<"vasculature -> "<<av<<endl;
+cout<<"----------------------"<<endl;
+cout<<"Generated side angle = "<<sideangle<<" degrees"<<endl;
+cout<<"----------------------"<<endl;
