Laboratoire RDP » Biophysics Team » Dandelion pappus morphing

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

 *                  GEOMETRY and MESH

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

// ----------------- external boundary

//border a1(t=0,1){x=lx/2+t*(-lx);y=0;label=1;}; //ok

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

border a1(t=0,1){x=lx/2+t*(lx/2-vthickness-lx/2);y=0;label=11;}; //ok

border a2(t=tB,tA){x=xC+R*cos(t);y=yC+R*sin(t);label=12;}; //ok

//border a2(t=0,1){x=rnect+t*(-2*rnect);y=0;label=12;}; //ok

border a3(t=0,1){x=-lx/2+vthickness+t*(-lx/2-(-lx/2+vthickness));y=0;label=13;}; //ok

border b1(t=0,1){x=-lx/2;y=t*(-ly);label=2;}; //ok

border c1(t=0,1){x=-lx/2+t*(-cwidth/2-vthickness+lx/2);y=-ly;label=31;}; //ok

border c2(t=0,1){x=-cwidth/2-vthickness+t*vthickness;y=-ly;label=32;}; //ok

border cc(t=0,1){x=-cwidth/2+t*cwidth; y=-ly; label=35;}; // bouchon

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

border d1(t=0,1){x=lx/2;y=-ly+t*(ly);label=4;}; //ok

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

/*

cout<<"xM="<<xM<<" yM="<<yM<<endl;

cout<<"xxM="<<xC+(R+hnect)*cos(tM)<<" yyM="<<yC+(R+hnect)*sin(tM)<<endl;

cout<<"xxN="<<xC+(R+hnect)*cos(tN)<<" yyN="<<yC+(R+hnect)*sin(tN)<<endl;

*/

border vide1(t=0,1){x=-cwidth/2+t*(xM+cwidth/2);y=-ly+t*(yM+ly);label=51;};//ok

border vide2(t=tM,tN){x=xC+(R+hnect)*cos(t);y=yC+(R+hnect)*sin(t);label=52;};

//border vide2(t=0,1){x=xM+(-2*xM)*t;y=yM;label=52;};

border vide3(t=0,1){x=-xM+t*(cwidth/2+xM);y=yM+t*(-ly-yM);label=53;};//ok

// ------------------- define the mesh

//border circle(t=2*pi,0){x=xC+R*cos(t);y=yC+R*sin(t);label=12;};

//mesh Th = buildmesh(circle(-nvertex*lx));

/*

plot(a1(nvertex*vthickness)+a2(nvertex*2*rnect)+a3(nvertex*vthickness)

                     + b1(nvertex*ly)

                     + c1(nvertex*lx/2) + c2(nvertex*vthickness)

                     + vide1(nvertex*ly) + vide2(nvertex*(lx-2*vthickness)) + vide3(nvertex*ly)  //hole

                     //+ cc(nvertex*cwidth)

                     + c3(nvertex*vthickness) + c4(nvertex*lx/2)

                     + d1(nvertex*ly), wait=1

                     );

*/

mesh Th = buildmesh (a1(nvertex*vthickness)+a2(nvertex*2*rnect)+a3(nvertex*vthickness)

                     + b1(nvertex*ly)

                     + c1(nvertex*lx/2) + c2(nvertex*vthickness)

                     + vide1(nvertex*ly) + vide2(nvertex*(lx-2*vthickness)) + vide3(nvertex*ly)  //hole

                     //+ cc(nvertex*cwidth)

                     + c3(nvertex*vthickness) + c4(nvertex*lx/2)

                     + d1(nvertex*ly)

                     );

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

 *                  DEFINITION OF DOMAINS

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

//vposition=-xM;

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>-hnect)&&(yy>D1*(xx-vthickness)+C1)&&(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));

}