
// Axisymmetric model of a simple artery with a membrane inside


//****************
//** Parameters **
//****************

rart=0.5e-3;	//Inner radius of the artery 
Lart=10e-3;	//Length of the artery
hm=0.1e-3;	//Thickness of the external membrane

Lmb=rart*3.0/4.0;	//Length of the membrane from the artery wall.
thmb=hm;		//Thickness of the membrane.


Lc=hm/5;	//Characteristic length for elements. (We will take 2 element on the membrane thickness)


xs[]={0,rart-Lmb,rart,(rart+hm)};	//Positions of the two lines for generating the surfaces.
nxs=#xs[];			// Number of divisions in xs

ys[]={0,(Lart-thmb)/2.0,(Lart+thmb)/2.0,Lart};
nys=#ys[];			// Number of divisions in ys
ndivsy[]={25,4,25};
nnds=#ndivsy[];			// Number of elements in y direction


//**************
//** Geometry **
//**************
//Points - 4
np=newp;
p0=np;
For xval In {0:nxs-1}
	Point(np)={xs[xval],0,0,Lc};
	np=newp;
EndFor

//Lines
For xval In {0:nxs-2}
	nl=newl;
	Line(nl)={p0+xval,p0+xval+1};
EndFor


//Areas
lines0[]={1,2,3};

linesact[]={lines0[]};
areas={};
lines={};
lslatd={};
lslati={};
For dv In {0:nnds-1}
	For ln In {0:#linesact[]-1}
		//ss[]=Extrude {0,ys[dv+1]-ys[dv],0} {Line{linesact[ln]}; Layers{ndivsy[dv]}; Recombine;};
		ss[]=Extrude {0,ys[dv+1]-ys[dv],0} {Line{linesact[ln]};};
		linesact[ln]=ss[0];
		areas[] +={ss[1]};
		lslatd[] += {ss[2]};
		lslati[] += {ss[3]};
		If (ln==0)
			axis[] +={ss[3]};
		EndIf
	EndFor
	membs[] += {ss[1]};
	lext[] += {ss[2]};
EndFor


membs[] +={areas[4]};
blood[]={areas[]};
blood[] -={membs[]};

lsfsi[]={10,8,18,20,34};


//Between two bodies, the normal vector is directed towards the body with lower ID number
Physical Surface(1)={blood[]};	//Blood
Physical Surface(2)={membs[]};	//Membrane


Physical Line(1)={1,2};		//Input fluid
Physical Line(2)={linesact[0],linesact[1]};		//Output fluid
Physical Line(3)={3};		//Fixed position at input
Physical Line(4)={linesact[2]};		//Fixed position at output

Physical Line(5)={axis[]};
Physical Line(6)={lext[]};
Physical Line(7)={lsfsi[]};


Mesh.CharacteristicLengthFromPoints=0;
Mesh.CharacteristicLengthFromCurvature=0;
Mesh.CharacteristicLengthExtendFromBoundary=0;

//Attractor
Field[1]=Attractor;
Field[1].EdgesList={lsfsi[]};
Field[1].NNodesByEdge=1000;

//Minimum mesh + exponentially increasing from the distance to the curve (Attractor Field)
Field[2]=MathEval;
fv=rart/50;
Field[2].F=Sprintf("%g+F1*0.1",fv);

Background Field=2;

RefineMesh;

