/*********************************************************************
 CircularTube:
 Inner diameter: 40 mm
 Outer diameter: 44 mm
 Length:         125 mm  
 *********************************************************************/

unit=0.001; //unit: m
lc = 0.8*unit; //lc: 1 mm: should be smaller than 0.6

//1/4Tube with 45 degree chamfering - Geometry modelling by rotation

//sq2 = Sqrt(2.0);
InnerRadius = 20*unit;//Practical size: 20
PlyThickness = 0.2270*unit;//Thickness: 0.2270
layer = 8;//Layer number
h = 40.*unit;// Short Tube Height

ratio = 5;

//Tube geometry:

//Contact region 

//Main Part

Point(1) = {0., 0., 0., ratio*lc};

r = InnerRadius;

h1 = h-PlyThickness*layer;

//Charmfer

R = r+PlyThickness*layer;

Point(2) = {0., 0., h1, lc};

Point(3) = {R, 0., h1, lc};
Point(4) = {R-PlyThickness, 0., h1, lc};
Point(5) = {R-PlyThickness, 0., h1+PlyThickness, lc};

l1 = newreg; Line(l1) = {3,4};
l2 = newreg; Line(l2) = {4,5};
l3 = newreg; Line(l3) = {5,3};

Lp = newreg; Line loop(Lp)={l1,l2,l3};
S1 = newreg; Surface(S1)={Lp};

Extrude {{0,0,1}, {0,0,0}, Pi/2} {Surface{S1}; Layers{20}; Recombine;}

KP = 0;
KeyPoint[KP] = 4;
KeyPoint[KP+1] = 5;

R=R-PlyThickness;

For n In {2:layer}

R=R-PlyThickness;

p1 = newp; Point(p1) = {R, 0., h1, lc};
p2 = newp; Point(p2) = {R, 0., h1+n*PlyThickness, lc};

l4 = -l2;
l1 = newreg; Line(l1) = {KeyPoint[KP],p1};
l2 = newreg; Line(l2) = {p1,p2};
l3 = newreg; Line(l3) = {p2,KeyPoint[KP+1]};

Lp = newreg; Line loop(Lp)={l1,l2,l3,l4};
S1 = newreg; Surface(S1)={Lp};

Extrude {{0,0,1}, {0,0,0}, Pi/2} {Surface{S1}; Layers{20}; Recombine;}

KP = KP+2;
KeyPoint[KP] = p1;
KeyPoint[KP+1] = p2;

EndFor

//Main Part

p3 = newp; Point(p3) = {r, 0., 0, ratio*lc};
p4 = newp; Point(p4) = {0., r, 0, ratio*lc};

ld1 = newreg; Circle(ld1) = {p3,1,p4};

ld3 = newreg; Line(ld3) = {238,p3};
ld4 = newreg; Line(ld4) = {241,p4};

Ldp = newreg; Line loop(Ldp)={ld1,-ld4,-197,ld3};
Sd1 = newreg; Surface(Sd1)={Ldp};

Transfinite Line{ld1} = 10;

Extrude{ Surface{Sd1}; Layers{{layer},{PlyThickness*layer}}; Recombine;}

Recombine Volume{1,9};
Recombine Volume{2,9};
Recombine Volume{3,9};
Recombine Volume{4,9};
Recombine Volume{5,9};
Recombine Volume{6,9};
Recombine Volume{7,9};
Recombine Volume{8,9};

//Define physical entities

//Define layers
//Layer_0
Physical Volume(1001)={1,3,6,8};

//Layer_90_1
Physical Volume(1002)={2,4};

//Layer_90_2
Physical Volume(1003)={5,7};

//Main Damage Part
Physical Volume(1111)={9};

//Impact Plate

r = InnerRadius+PlyThickness*layer;
ht = 40.*unit;

v = 10; //loading velocity 10 m/s
ImpH = 1.e-11*1*v; //ImpactHeight

X_contact=-2.*r;
Y_contact=-2.*r;
Z_contact=ht+ImpH;
L_contact=4.*r;
Point(11001) = {0,0,Z_contact,r};
Point(11002) = {X_contact+L_contact,0,Z_contact,r};
Point(11003) = {X_contact+L_contact,Y_contact+L_contact,Z_contact,r};
Point(11004) = {0,Y_contact+L_contact,Z_contact,r};
Line(10001) = {11001, 11002};
Line(10002) = {11002, 11003};
Line(10003) = {11003, 11004};
Line(10004) = {11004, 11001};
Line Loop(10001) = {10004, 10001, 10002, 10003};
Physical Point(11001) = {11001};
Plane Surface(10001) = {10001};
Transfinite Line {10001,10002,10003,10004} = 2 Using Progression 1;
Transfinite Surface {10001};
Recombine Surface {10001};
Physical Surface(10001) = {10001};


