hybrid extruded/unstructured from ruth

benchmarks/3d/coil_crack/KDS70-2_dat.pro

+micra = 1e-6;
+
+// Differential pencil probe KDS 70-2
+// Producer: Rohmann GmbH, Germany
+// Coils with ferrite core and shielding
+// 2 coils with 35 turns (wire diameter 0.05 mm ) on 'half D' core (0.5 mm)
+// Primary coil with 21 turns (wire diameter 0.1 mm)- global diameter 1.6 mm
+// lift-off: 0.92 mm 
+// (distance between the end of probe - in contact with sample and first turn)
+
+Rcore = 450 * micra ;
+w = 50 * micra ;
+turns = 11 ;
+turns_in = turns ;
+turns_out = turns;
+w_in = w ;
+w_out_ = w ;
+d_coils = 0 ;
+
+hc = 620 * micra ; 
+sc = w;
+
+dist = 380 * micra ;
+
+w_air = (230+170+70)*micra/2 ;
+w_blindage = (400+290)*micra/2 ;
+
+liftoff = 100 * micra ; 
+
+//-----------------------------------------------------------------------------------
+COILOUT = 1000 ;
+COILIN1 = 1001 ;
+COILIN2 = 1002 ;
+
+SKINCOILIN1 = 1011 ;
+SKINCOILIN12 = 1012 ;
+
+SKINCOILIN2 = 1021 ;
+SKINCOILIN22 = 1022 ;
+
+SKINCOILOUT = 1100;
+SKINCOILOUT2 = 1200;
+
+CUTCOILIN1 = 1111 ; 
+CUTCOILIN2 = 1112 ;
+CUTCOILOUT = 1113 ;
+ 
+CORE1 = 2000 ;
+CORE2 = 2001 ;
+
+CORE1_0 = 2010 ;
+CORE2_0 = 2011 ;
+
+CORE1_1 = 2012 ;
+CORE2_1 = 2013 ;
+
+SKINCORE1 = 2110 ;
+SKINCORE2 = 2111 ;
+
+GAP = 2222 ;
+
+SHIELD = 5000;
+

benchmarks/3d/coil_crack/KDS70-2z1.geo

+Include "KDS70-2_dat.pro";
+Include "Probe_position.pro";
+
+ncore1 = 3 ;
+ncore2 = 4 ;
+
+s0 = 1 ;
+
+lcc = Rcore/3/s0 ;
+
+cte = news-1 ;
+
+cen0 = newp ; Point(cen0) = {xp0, yp0  , 0,lcc};
+cen1 = newp ; Point(cen1) = {xp0, yp0+w, 0,lcc};
+
+pc1 = newp ; Point(pc1) = {xp0, yp0+w, Rcore, lcc};
+pc2 = newp ; Point(pc2) = {xp0, yp0+w+Rcore, 0, lcc };
+pc3 = newp ; Point(pc3) = {xp0, yp0+w, -Rcore, lcc };
+
+c1 = newl ; Circle(c1) = {pc1,cen1,pc2};
+c2 = newl ; Circle(c2) = {pc2,cen1,pc3};
+
+pc4 = newp ; Point(pc4) = {xp0, yp0+w,Rcore+w,lcc};
+pc5 = newp ; Point(pc5) = {xp0, yp0+2*w+Rcore, 0,lcc};
+pc6 = newp ; Point(pc6) = {xp0, yp0+w, -(Rcore+w),lcc};
+
+c3 = newl ; Circle(c3) = {pc4,cen1,pc5};
+c4 = newl ; Circle(c4) = {pc5,cen1,pc6};
+
+pc7 = newp ; Point(pc7) = {xp0, yp0, Rcore+w,lcc};
+pc8 = newp ; Point(pc8) = {xp0, yp0,-(Rcore+w),lcc};
+
+ln1 = newl ; Line(ln1) ={pc1,pc3}; 
+ln2 = newl ; Line(ln2) ={pc7,pc8}; 
+ln3 = newl ; Line(ln3) ={pc4,pc7}; 
+ln4 = newl ; Line(ln4) ={pc6,pc8}; 
+
+surfcore1 = newl ;
+Line Loop(surfcore1) = {ln1,-c2,-c1};
+Plane Surface(surfcore1) = {surfcore1};
+
+surfcoil1 = newl ;
+Line Loop(surfcoil1) ={ln2,-ln4,-c4,-c3,ln3};
+Plane Surface(surfcoil1) = {surfcoil1,surfcore1};
+
+surf[]=Symmetry {0.0,1.0,0.0,-yp0}{Duplicata { Surface{surfcore1}; }};
+surfcore2 = surf[0];
+surf[]=Symmetry {0.0,1.0,0.0,-yp0}{Duplicata { Surface{surfcoil1}; }};
+surfcoil2 = surf[0];
+
+
+pc1_ = newp ; Point(pc1_) = {xp0, yp0+w, Rcore + w + w,lcc};
+pc2_ = newp ; Point(pc2_) = {xp0, yp0+2*w + Rcore + w,0,lcc};
+pc3_ = newp ; Point(pc3_) = {xp0, yp0+w,-(Rcore+w +w),lcc};
+
+c1_ = newl ; Circle(c1_) = {pc1_,cen1,pc2_};
+c2_ = newl ; Circle(c2_) = {pc2_,cen1,pc3_};
+
+lin[]=Symmetry {0.0,1.0,0.0,-yp0}{Duplicata { Line{c1_}; }};
+c3_ = lin[0];
+lin[]=Symmetry {0.0,1.0,0.0,-yp0}{Duplicata { Line{c2_}; }};
+c4_ = lin[0];
+
+ln5 = newl ; Line(ln5)={pc1_,pc1_+3};
+ln6 = newl ; Line(ln6)={pc3_,pc3_+6};
+
+pa1 = newp ; Point(pa1) = {xp0, yp0, Rcore +2*w+w_air, lcc};
+pa2 = newp ; Point(pa2) = {xp0, yp0+Rcore +2*w+w_air, 0, lcc };
+pa3 = newp ; Point(pa3) = {xp0, yp0, -(Rcore +2*w+w_air), lcc};
+pa4 = newp ; Point(pa4) = {xp0, yp0-(Rcore +2*w+w_air), 0, lcc };
+
+ca1 = newl ; Circle(ca1) = {pa1,cen0,pa2};
+ca2 = newl ; Circle(ca2) = {pa2,cen0,pa3};
+ca3 = newl ; Circle(ca3) = {pa3,cen0,pa4};
+ca4 = newl ; Circle(ca4) = {pa4,cen0,pa1};
+
+pb1 = newp ; Point(pb1) = {xp0, yp0, Rcore +2*w+w_air+w_blindage, lcc};
+pb2 = newp ; Point(pb2) = {xp0, yp0+Rcore +2*w+w_air+w_blindage, 0, lcc };
+pb3 = newp ; Point(pb3) = {xp0, yp0, -(Rcore +2*w+w_air+w_blindage), lcc};
+pb4 = newp ; Point(pb4) = {xp0, yp0-(Rcore +2*w+w_air+w_blindage), 0, lcc };
+
+cb1 = newl ; Circle(cb1) = {pb1,cen0,pb2};
+cb2 = newl ; Circle(cb2) = {pb2,cen0,pb3};
+cb3 = newl ; Circle(cb3) = {pb3,cen0,pb4};
+cb4 = newl ; Circle(cb4) = {pb4,cen0,pb1};
+
+scoilout = news ;
+Line Loop(scoilout) = {c3,c4,ln4,17+cte,18+cte,19+cte,20+cte,-ln3};
+scoilout_ = news;
+Line Loop(scoilout_) = {c1_,c2_,ln6,-c3_,-c4_,-ln5};
+Plane Surface(scoilout) = {scoilout_,scoilout};
+
+lair = news;
+Line Loop(lair) = {ca1,ca2,ca3,ca4};
+sair = news;
+Plane Surface(sair) = {lair,scoilout_};
+
+sblin=news;
+Line Loop(sblin) = {cb1,cb2,cb3,cb4};
+Plane Surface(sblin) = {sblin,lair};
+
+//__________________________________________________________________________
+//__________________________________________________________________________
+
+//Core
+
+dx = -dist ; dy = 0.; dz = 0.;
+
+vol[]=Extrude Surface{surfcore1, {dx,dy,dz}} {Layers {{ncore1},{1}};};;
+vCore10 = vol[1];
+vol[]=Extrude Surface{surfcore2, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vCore20 = vol[1];
+
+dx = hc ;  
+vol[]=Extrude Surface{surfcore1, {dx,dy,dz}} {Layers {{ncore2},{1}};};; 
+vCore11 = vol[1];
+surfcore11 = news-1;
+vol[]=Extrude Surface{surfcore2, {dx,dy,dz}} {Layers {{ncore2},{1}};};; 
+vCore21 = vol[1];
+
+dx = dist ;  
+surfcore21 = news-1;
+vol[]=Extrude Surface{surfcore11, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vCore12 = vol[1];
+vol[]=Extrude Surface{surfcore21, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vCore22 = vol[1];
+
+
+//Coils_in
+dx = hc ;
+vol[]=Extrude Surface{surfcoil1, {dx,dy,dz}} {Layers {{ncore2},{1}};};; 
+vCoilIn1 = vol[1]; surfcoil1up = news-1;
+vol[]=Extrude Surface{surfcoil2, {dx,dy,dz}} {Layers {{ncore2},{1}};};; 
+vCoilIn2 = vol[1]; surfcoil2up = news-1;
+
+dx = -dist;
+vol[]=Extrude Surface{surfcoil1, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vAir0CoilIn1 = vol[1];
+vol[]=Extrude Surface{surfcoil2, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vAir0CoilIn2 = vol[1];
+
+dx = dist;
+vol[]=Extrude Surface{surfcoil1up, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vAir1CoilIn1 = vol[1];
+vol[]=Extrude Surface{surfcoil2up, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vAir1CoilIn2 = vol[1];
+
+//Coil_out
+dx = -dist ; 
+vol[]=Extrude Surface{scoilout, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vAir0CoilOut0 = vol[1];
+dx = hc ;
+vol[]=Extrude Surface{scoilout, {dx,dy,dz}} {Layers {{ncore2},{1}};};; 
+vCoilOut = vol[1];
+dx=dist;
+vol[]=Extrude Surface{news-1, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vAir1CoilOut0 = vol[1];
+
+//Air in
+dx = -dist;
+vol[]=Extrude Surface{sair,{dx,dy,dz} } {Layers {{ncore1},{1}};};;
+vAirIn0 = vol[1];
+dx=hc;
+vol[]=Extrude Surface{sair, {dx,dy,dz}} {Layers {{ncore2},{1}};};; 
+vAirIn = vol[1];
+dx=dist;
+vol[]=Extrude Surface{news-1, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vAirIn1 = vol[1];
+
+//Blindage
+dx = -dist;
+vol[]=Extrude Surface{sblin, {dx,dy,dz}} {Layers {{ncore1},{1}};};;
+vBlin0 = vol[1];
+dx=hc;
+vol[]=Extrude Surface{sblin, {dx,dy,dz}} {Layers {{ncore2},{1}};};; 
+vBlin1 = vol[1];
+
+sblin1 = news-1;
+dx=dist;
+vol[]=Extrude Surface{sblin1, {dx,dy,dz}} {Layers {{ncore1},{1}};};; 
+vBlin2 = vol[1];
+
+
+//-----------------------------------------------------------------------------------
+//-----------------------------------------------------------------------------------
+
+Physical Volume(COILIN1) = {vCoilIn1};
+Physical Volume(COILIN2) = {vCoilIn2};
+Physical Volume(COILOUT) = {vCoilOut};
+Physical Volume(CORE1) = {vCore11};
+Physical Volume(CORE2) = {vCore21};
+Physical Volume(CORE1_0) = {vCore10};
+Physical Volume(CORE2_0) = {vCore20};
+Physical Volume(CORE1_1) = {vCore12};
+Physical Volume(CORE2_1) = {vCore22};
+
+Physical Volume(GAP)={vAirIn,vAirIn0,vAirIn1,
+		      vAir0CoilIn1,vAir0CoilIn2,vAir1CoilIn1,vAir1CoilIn2,
+		      vAir0CoilOut0,vAir1CoilOut0};
+Physical Volume(SHIELD) = {vBlin0,vBlin1,vBlin2};
+
+Physical Surface(SKINCOILIN1) = 
+{89+cte,85+cte,81+cte,158+cte,170+cte,166+cte,162+cte,154+cte,10+cte,183+cte};
+Physical Surface(SKINCOILIN12) ={89+cte,85+cte,81+cte};
+Physical Surface(SKINCOILIN2) = 
+{106+cte,102+cte,204+cte,208+cte,98+cte,154+cte,212+cte,200+cte,15+cte,225+cte};
+Physical Surface(SKINCOILIN22) ={106+cte,102+cte,98+cte};
+
+Physical Surface(CUTCOILIN1) = {90+cte};//up
+Physical Surface(CUTCOILIN2) = {107+cte};//up
+
+Physical Surface(SKINCOILOUT) = 
+{488+cte,484+cte,504+cte,500+cte,496+cte,492+cte,162+cte,166+cte,208+cte,204+cte,
+200+cte,158+cte,170+cte,212+cte,537+cte,35+cte};
+Physical Surface(SKINCOILOUT2) = 
+{162+cte,166+cte,170+cte,212+cte,208+cte,204+cte,200+cte,158+cte};
+
+
+Physical Surface(CUTCOILOUT) = {9+cte,10+cte,11+cte,15+cte};//down
+
+Physical Surface(SKINCORE1)=
+{124+cte,123+cte,119+cte,115+cte,89+cte,85+cte,81+cte,47+cte,55+cte,51+cte,56+cte};
+Physical Surface(SKINCORE2)= 
+{141+cte,132+cte,136+cte,140+cte,98+cte,102+cte,106+cte,64+cte,68+cte,72+cte,73+cte};
+
+
+//-----------------------------------------------------------------------------------
+//-----------------------------------------------------------------------------------
+Color Gold {Surface{89+cte,85+cte,81+cte,158+cte,170+cte,166+cte,162+cte,154+cte,10+cte,183+cte};}
+Color Gold {Surface{106+cte,102+cte,204+cte,208+cte,98+cte,154+cte,212+cte,200+cte,15+cte,225+cte};}
+Color Gold {Volume{vCoilIn1,vCoilIn2};}
+
+Color Red {Surface{488+cte,484+cte,504+cte,500+cte,496+cte,492+cte,162+cte,166+cte,208+cte,204+cte,
+		   200+cte,158+cte,170+cte,212+cte,537+cte,35+cte};}
+Color Red {Volume{vCoilOut};}
+
+Color DeepSkyBlue1 {Surface{124+cte,123+cte,119+cte,115+cte,89+cte,85+cte,81+cte,47+cte,55+cte,51+cte,56+cte};}
+Color DeepSkyBlue1 {Surface{141+cte,132+cte,136+cte,140+cte,98+cte,102+cte,106+cte,64+cte,68+cte,72+cte,73+cte};}
+Color DeepSkyBlue1 {Volume{vCore10,vCore11,vCore12,vCore20,vCore21,vCore22};}
+
+Color MidnightBlue {Volume{vBlin0,vBlin1,vBlin2};}
+
+surftotprobe = news ;
+Surface Loop(surftotprobe) = 
+{891+cte,862+cte,820+cte,-778-cte,-782-cte,824+cte,828+cte,-786-cte,-790-cte,
+832+cte,874+cte,870+cte,866+cte,-807-cte,-661-cte,-465-cte,267+cte,-309-cte,
+-73-cte,56+cte,765+cte,609+cte,-351-cte,393+cte,141+cte,-124-cte};

benchmarks/3d/coil_crack/Piece1_4z1b.geo

+Include "Piece_dat.pro";
+//Defect
+
+LAYERS = 1 ;
+
+//   length x width x depth (mm)
+//B  1 x 0.2 x 0.1
+//D  1 x 0.4 x 0.2
+length = 1 * mm ;
+depth = 0.1 * mm ;//Shape: half disk
+width = 2 * depth ;
+
+
+//Mesh.Algorithm = 2 ; // 2D mesh algorithm (1=isotropic, 2=anisotropic, 3=triangle)
+//Mesh.Algorithm3D = 1; // 3D mesh algorithm (1=isotropic, 4=netgen, 5=tetgen)
+Mesh.RandomFactor = 0.0002 ;
+
+Geometry.AutoCoherence = 1;
+
+s = 1. ;
+
+Freq = 100.e3 ;
+mu0 = 4.e-7 * Pi ;
+sigma = 1.798e6 ;//Titanium
+delta =  Sqrt[1/(mu0*Freq*sigma*Pi)] ;// skin depth
+
+Printf("Freq %g delta %g",Freq,delta);
+
+pdwn = h1/2/s ;
+
+pflaw = width/4/s;
+pflaw2 = width/2/s;
+
+ph1 = h1/3/s ;
+ph2 = 14*pflaw ;
+ph3 = 14*pflaw ;
+ph3_ = 14*pflaw ;
+
+pBox = 4*pdwn/s;
+
+AngRot = Pi/16 ;
+AngRot_ = -2*AngRot;
+
+nDefect = Ceil[length/pflaw]; //20 ;
+nAngle  = Ceil[AngRot*Rdwn/ph1] ;
+
+Rlimup  = 0.8*Rup ;
+Rlimdwn = 0.75*Rdwn;
+
+// Box
+xBox1 = Rlimdwn*1.2 ;
+xBox2 = 0.4*Rup ; //-xBox1/3 ;
+yBox1 = 2.*h2 ;
+yBox2 = -0.8*h2 ;
+
+zB[] = {0.35*xBox1,-0.35*xBox1};
+
+//---------------------------------------------------------------------
+//---------------------------------------------------------------------
+
+//Defect Zone 1, centered around (h3+h2)/2
+rr = Rup ;
+aa = 2*Asin[width/2./rr];
+zz = rr*Sin[aa/2];
+xx = rr*Cos[aa/2];
+
+//Piece
+cen[0] = newp ; Point(cen[0])={0, 0, 0, pdwn};
+cen[1] = newp ; Point(cen[1])={Rc,h2,0, ph2};//curve
+
+cen[2] = newp ; Point(cen[2])={0,h1,0,ph1};
+cen[3] = newp ; Point(cen[3])={0,h2,0,ph2};
+cen[4] = newp ; Point(cen[4])={0,h3,0,ph3};
+
+p[0] = newp ; Point(p[0])={Rlimdwn, 0, 0,pdwn};
+p[1] = newp ; Point(p[1])={Rlimdwn, h1, 0,ph1};
+p[2] = newp ; Point(p[2])={Rc, h1, 0, ph1};
+
+p[3] = newp ; Point(p[3])={Rc-R10,h2,0,ph2};
+
+p[4] = newp ; Point(p[4])={Rup, h3-length, 0, ph3};
+p[5] = newp ; Point(p[5])={Rup, h3, 0, ph3};
+
+p[6] = newp ; Point(p[6])={Rlimup, h3, 0, ph3_};
+p[7] = newp ; Point(p[7])={Rlimup,  0, 0, pdwn};
+
+lp[0] = newl ; Line(lp[0])={p[6],p[7]};//axis
+
+lp[1] = newl ; Line(lp[1])={p[0],p[1]};
+lp[2] = newl ; Line(lp[2])={p[1],p[2]};
+lp[3] = newl ; Circle(lp[3])={p[2],cen[1],p[3]};
+lp[4] = newl ; Line(lp[4])={p[3],p[4]};
+lp[5] = newl ; Line(lp[5])={p[4],p[5]};
+lp[6] = newl ; Line(lp[6])={p[5],p[6]};
+lp[7] = newl ; Line(lp[7])={p[7],p[0]};
+
+surf0 = news ; 
+Line Loop(surf0) = {lp[1],lp[2],lp[3],lp[4],lp[5],lp[6],lp[0],lp[7]};
+Plane Surface(surf0) = {surf0};
+
+surf[]= Rotate {{0.0,1.0,0.0},{0.0,0.0,0.0},aa/2} {Surface{surf0};};
+sRight =surf[0]; 
+surf[]= Rotate {{0.0,1.0,0.0},{0.0,0.0,0.0},-aa} { Duplicata{Surface{sRight};} };
+sLeft = surf[0];
+
+//Defect
+cend1 = newp ; Point(cend1)={Rup,h3-length,0,pflaw};
+cend2 = newp ; Point(cend2)={0,h3-length,0,pflaw};
+
+pd[0] = 10 ;
+pd[1] = 28 ;
+
+/*
+pd[0] = newp ; Point(pd[0])={xx,h3-length,zz,pflaw};
+pd[1] = newp ; Point(pd[1])={xx,h3-length,-zz,pflaw};
+*/
+
+pd[2] = newp ; Point(pd[2])={Rup-depth,h3-length,0,pflaw};
+
+
+cd[0]= newl ; Circle(cd[0])={pd[0],cend1,pd[2]};
+cd[1]= newl ; Circle(cd[1])={pd[2],cend1,pd[1]};
+cd[2] = newl ;Circle(cd[2]) = {pd[0],cend2, pd[1]} ;
+
+sDefect = news ;
+Line Loop(sDefect) = {cd[2],-cd[1],-cd[0]};
+Plane Surface(sDefect) = {sDefect};
+
+If (LAYERS)
+vol[]=Extrude Surface{sDefect, {0,length,0}} 
+//;;
+{Layers {{nDefect},{1}};};;
+vDefect = vol[1];
+
+vol[] = Extrude Surface {sRight, {0.0,1.0,0.0}, {0.0,0.0,0.0},AngRot-aa/2}// ;;
+  {Layers {{nAngle},{1}};}; ;
+vPieceR=vol[1];
+
+vol[] = Extrude Surface {sLeft, {0.0,1.0,0.0}, {0.0,0.0,0.0},-(AngRot-aa/2)}// ;;
+  {Layers {{nAngle},{1}};}; ;
+vPieceL=vol[1];
+EndIf
+
+If (!LAYERS)
+vol[]=Extrude Surface{sDefect, {0,length,0}} ;;
+vDefect = vol[1];
+
+vol[] = Extrude Surface {sRight, {0.0,1.0,0.0}, {0.0,0.0,0.0},AngRot-aa/2} ;;
+vPieceR=vol[1];
+
+vol[] = Extrude Surface {sLeft, {0.0,1.0,0.0}, {0.0,0.0,0.0},-(AngRot-aa/2)};;
+vPieceL=vol[1];
+EndIf
+
+
+
+
+
+
+
+Printf("vPieceR %g vPieceL %g",vPieceR,vPieceL);
+
+Circle(124) = {14,1,6};
+Circle(125) = {15,3,7};
+Circle(126) = {19,3,8};
+Circle(127) = {24,4,9};
+Circle(128) = {36,5,12};
+Circle(129) = {40,1,13};
+Line Loop(130) = {124,2,-125,-11};
+Ruled Surface(131) = {130};
+Line Loop(132) = {126,4,-127,-13};
+Ruled Surface(133) = {132};
+Line Loop(134) = {127,5,21,-14};
+Ruled Surface(135) = {134};
+Line Loop(136) = {129,-1,-128,17};
+Ruled Surface(137) = {136};
+Line Loop(138) = {125,3,-126,-12};
+Plane Surface(139) = {138};
+Line Loop(140) = {124,-8,-129,18};
+Plane Surface(141) = {140};
+Line Loop(142) = {16,128,-7,-26,-25};
+Plane Surface(143) = {142};
+
+Surface Loop(144) = {131,-141,-9,139,133,135,-22,34,38,-143,10,-137};
+Volume(145) = {144};
+vPieceC=145;
+
+//Include "Probe.geo";
+Include "KDS70-2z1.geo";
+
+For k In {0:1}
+Point(newp)={ xBox2, yBox2,  zB[k], pBox};
+Point(newp)={ xBox1, yBox2,  zB[k], pBox};
+Point(newp)={ xBox1, yBox1,  zB[k], pBox};
+Point(newp)={ xBox2, yBox1,  zB[k], pBox};
+
+lbox[4*k] = newl ; Line(lbox[4*k])={newp-4,newp-3};
+lbox[4*k+1] = newl ; Line(lbox[4*k+1])={newp-3,newp-2};
+lbox[4*k+2] = newl ; Line(lbox[4*k+2])={newp-2,newp-1};
+lbox[4*k+3] = newl ; Line(lbox[4*k+3])={newp-1,newp-4};
+
+sbox[k] = news ; Line Loop(sbox[k]) = {newreg-4:newreg-1};
+Plane Surface(sbox[k]) = {sbox[k]};
+EndFor
+
+
+lbv1 =newl; Line(lbv1)={newp-8,newp-4};
+lbv2 =newl; Line(lbv2)={newp-7,newp-3};
+lbv3 =newl; Line(lbv3)={newp-6,newp-2};
+lbv4 =newl; Line(lbv4)={newp-5,newp-1};
+
+sbox1 = news ; Line Loop(news) = {lbox[0],lbv2,-lbox[4],-lbv1};
+Plane Surface(sbox1) = {sbox1};
+sbox2 = news ; Line Loop(news) = {lbox[1],lbv3,-lbox[5],-lbv2};
+Plane Surface(sbox2) = {sbox2};
+sbox3 = news ; Line Loop(news) = {lbox[2],lbv4,-lbox[6],-lbv3};
+Plane Surface(sbox3) = {sbox3};
+sbox4 = news ; Line Loop(news) = {lbox[3],lbv1,-lbox[7],-lbv4};
+Plane Surface(sbox4) = {sbox4};
+
+Surface Loop(news) =  {sbox2,-sbox[0],sbox1,sbox[1],sbox3,sbox4};
+Surface Loop(news) = 
+{52,-131,141,80,81,56,-139,-133,60,64,-135,30,-110,-114,143,137,76,72,68,-118,-122,-94,-98,-102,-106,-123,39}; 
+vAir= newv;
+Volume(vAir) = {news-2,news-1,surftotprobe};//Rest Air
+
+Characteristic Length{10,11,28,32}= pflaw ;//Watch out! Changing characteristic lenths of the flaw at the end
+//Characteristic Length{66,93}= ph3 ;
+
+//Transfinite Line{cd[0],cd[1],cd[2]} = 5 ;
+//Transfinite Line{15,6} = nDefect+1 ;
+
+//-------------------------------------------------------------------
+
+Physical Volume (CRACK) = {vDefect};
+Physical Volume(PIECE)={vPieceR,vPieceL,vPieceC};
+
+Physical Surface(SKINPIECE) =
+{72,143,-114,-118,137,141,-131,52,56,-139,-133,60,64,-135,22,-34,-38,68,81,76,80,-110,-123,-94,-98,-102,-106,-122};
+Physical Surface(SKINPIECEWCRACK) =
+{52,-131,141,80,81,56,-139,-133,60,64,-135,30,-110,-114,143,137,76,72,68,-118,-122,-94,-98,-102,-106,-123,39}; 
+
+Physical Volume(AIR) = {vAir};
+Physical Surface(SURFACEGH0) = {sbox[0],-sbox1,-sbox2,-sbox3,-sbox4,-sbox[1]};
+
+

benchmarks/3d/coil_crack/Piece_dat.pro

+mm = 1e-3;
+
+//Piece
+Rup = 64*mm/2 ;
+Rdwn = 121*mm/2 ;
+Rc = 84*mm/2;
+R10 = 10 * mm ;
+
+h1 = 5*mm ;
+h2 = 15*mm ;
+h3 = 20*mm ;
+
+//Defect
+//   length x width x depth (mm)
+//A  0.4 x 0.2 x 0.1
+//B  1 x 0.2 x 0.1
+//C  2 x 0.2 x 0.1
+//D  1 x 0.4 x 0.2
+//E  2 x 0.4 x 0.2
+
+length = 1. * mm ;
+depth = 0.1 * mm ;//Shape: half disk
+width = 2 * depth ;
+
+//Box
+xBox1 = Rdwn*1.5 ;
+xBox2 = -xBox1 ;
+yBox1 = 2.5*h3 ;
+yBox2 = -h3 ;
+
+L = (Rdwn-Rc)/2 ;
+
+//--------------------------------------------------------
+//Physical geometry
+
+PIECE = 3000; 
+SKINPIECE = 3001 ;
+SKINPIECEWCRACK = 3002 ;
+CRACK = 3333; 
+
+AIR = 4000 ;
+SURFACEGH0 = 4001 ;

benchmarks/3d/coil_crack/Probe_position.pro

+//Piece1_4z1.geo
+/*
+xp0 = dist+liftoff + Rup ; 
+yp0 = (h2+h3)/2 ;
+*/
+
+//Piece1_4z1b.geo
+
+xp0 = dist+liftoff + Rup ; 
+yp0 = h3 ; //Edge
+//yp0 = h3-length/2 ;//Center defect
+
+//Piece1_4z3.geo
+/*
+xp0 = Rc+length/2 ;
+yp0 = dist + liftoff + h1 ;
+*/

benchmarks/3d/lam3d_shell2_indxz.geo

+
+//BEGIN PARAMS
+
+//Include "dimensions_dat.pro";
+
+// Dimensions lamination stack
+
+//Include "freq3d.dat";
+
+u = 1e-3;
+
+nlam = 10 ;
+
+d = 0.5*u ; //thickness of laminations
+lambda = 0.9 ;//90%-98%
+
+e = d*(1-lambda)/lambda ;//thickness isolation
+
+W2 = 50*u/4 ; //Half thickness of stack
+D2 = 1. * u ; //Half gap 
+
+T = nlam*(d+e) ;
+
+wind = 4*u ;
+
+//W2 = T ;
+//dd = W2*Sqrt(2)+2*d ;
+
+
+dd = W2*1.5 ;//Radius inductor
+dd2 = W2*1.2 ;//indxz
+
+rind = dd;
+rind2 = dd2;
+
+xBox = dd*1.4 ;
+
+Rint = xBox ;
+Rext = xBox*1.2 ;
+
+x0 = 0 ; y0 = 0 ; z0 = 0; 
+z1 = (nlam-1)*(d+e)+d/2 ;
+
+AXISY = 100;
+AXISX = 101;
+LIM_UP    = 102;
+LIM_RIGHT = 103;
+
+SKININD = 200;
+IND = 500 ;
+CUTIND = 400 ;
+CUTINDAIR = 300;
+
+ELEC0 = 600 ;
+ELEC1 = 601 ;
+
+LAMINATION = 1000;
+ISOLATION  = 2000;
+
+AIR = 3000;
+AIRINF = 3001;
+
+SURFACEGINF  = 4000;
+SURFACEGE0 = 4001;//Side
+SURFACEGH0 = 4002;//Symmetry 
+SURFACEAIR = 4003;
+
+SURFBACK = 4010;
+SURFLEFT = 4011;
+SURFFRONT = 4012;
+SURFRIGHT = 4013;
+
+SURFBACKI = 4020;
+SURFLEFTI = 4021;
+SURFFRONTI = 4022;
+SURFRIGHTI = 4023;
+
+SURFUP = 4030;
+SURFDWN = 4031;
+
+SKINLAM = 5000;
+
+//END PARAMS
+
+//-----------------------------------------------------------------
+
+HOMO = 0 ;
+GAP = 1 ;
+
+// Some gmsh parameters
+
+Geometry.OldNewReg = 0;
+Geometry.AutoCoherence = 1;
+
+//Mesh.Algorithm = 1 ;
+
+TRANSIND = 1 ;
+
+s = 1.5 ;
+
+// L a m i n a t i o n s :
+//------------------------------------------------------------------
+
+lc     = s * W2/16;//W2/20;
+lcind  = s * wind/5;
+lcind2 = lcind ; 
+
+lcbox    = lc*2;
+lcboxi    = lcind ;
+lcboxout = lc*2;
+
+ndive = 1 ;
+
+If(HOMO)
+ndivlam = 2 ;
+ndivlam_half = 1 ;
+EndIf
+If(!HOMO)
+ndivlam = 8 ;
+EndIf
+
+
+
+//nind = 20;
+nind = Ceil[Pi*dd/2/lcind]+1 ;
+
+valProgression = 1. ;
+valBump = 0.1 ;
+//valBump = 1 ;
+
+// Print Vector (for checking numbers)
+//  IN:
+//    vector[] : list of quantities to be printed
+// OUT:
+//    /
+////////
+Function PrintVector
+vector_size_ = #vector[];
+Printf("Vector: size: %g", vector_size_);
+For i In {0:vector_size_-1}
+  Printf(" %g", vector[i]);
+EndFor 
+Return
+/* use e.g.:
+vector[]=yourlist[]; Call PrintVector;
+ */
+
+
+//---------------------------------------------------------------------
+
+x0 = 0 ; y0 = 0 ; z0 = 0; 
+
+p0 = newp ; Point(p0)={x0,y0,z0,lc};
+p1 = newp ; Point(p1)={x0+W2,y0,z0,lc};
+
+p2 = newp ; Point(p2)={x0+W2,y0+W2,z0,lc};
+p3 = newp ; Point(p3)={x0,y0+W2,z0,lc};
+
+l0 = newl ; Line(l0) = {p0,p1};  
+l1 = newl ; Line(l1) = {p1,p2};  
+l2 = newl ; Line(l2) = {p2,p3};  
+l3 = newl ; Line(l3) = {p3,p0};  
+
+surflam0 = news;
+Line Loop (surflam0) = {l0,l1,l2,l3};
+Plane Surface(surflam0) = {surflam0};
+
+Transfinite Line{l0,l3} = Ceil[W2/lc]+1 Using Bump valBump ;
+Transfinite Line{l1,l2} = Ceil[W2/lc]+1 Using Bump valBump ;
+Transfinite Surface{surflam0} = {p0,p1,p2,p3};
+
+vollam[]={}; voliso[]={};
+lback[]={}; lleft[]={};lfront[]={};lright[]={};
+lbacki[]={};llefti[]={};lfronti[]={};lrighti[]={};
+lup[] = {};ldwn[] = {};
+lupi[]={};
+
+If(!HOMO)
+bumpvector[] = {}; bumpvector1[] = {};
+lvector[] = {}; lvector1[] = {};
+
+alpha = 1.5 ; cte  = 0.;
+
+For i In{0:ndivlam/2-1}
+cte += alpha^i; 
+EndFor
+
+a = T/2/cte;
+
+For i In {0:ndivlam/2-1}
+lvector1[] += 1 ;
+testv = a*alpha^i/(T/2)/2;
+If(i>0)
+  testv += bumpvector1[i-1];
+EndIf
+bumpvector1[] += testv   ;
+EndFor
+
+For i In {0:#bumpvector1[]-1}//Half lamination
+bumpvector2[i] = 2*bumpvector1[i];
+Printf("%g %g",lvector1[i],bumpvector2[i]);
+EndFor
+
+bumpvector[]={bumpvector1[]};//Complete lamination
+For i In {0:ndivlam/2-1}
+If(i<ndivlam/2-1)
+bumpvector[ndivlam/2+i] = 
+bumpvector[ndivlam/2+i-1]+ bumpvector[ndivlam/2-1-i]-bumpvector[ndivlam/2-2-i] ;
+EndIf
+
+If(i==ndivlam/2-1)
+bumpvector[ndivlam/2+i] = 1;
+EndIf
+EndFor
+lvector[] ={lvector1[],lvector1[]};
+
+ldwn[]+= surflam0;//Watch Out: this is half a lamination
+vol[]= Extrude Surface { surflam0, {0,0,d/2}} {Layers {ndivlam/2,1} ; };;
+vollam[] += vol[1];
+lup[] += news-1;
+lfront[]  += news-5 ;
+lright[]  += news-4 ;
+lback[] += news-3 ;
+lleft[] += news-2 ;
+
+For i In {1:nlam-1}
+vol[] = Extrude Surface {news-1, {0,0,e}} {Layers {{ndive},{1}} ; };;
+voliso[]  += vol[1];
+lupi[]    += news-1 ;
+lfronti[] += news-5 ;
+lrighti[] += news-4 ;
+lbacki[]  += news-3 ;
+llefti[]  += news-2 ;
+
+ldwn[]+= news-1 ;
+vol[] = Extrude Surface { news-1, {0,0,d}} {Layers {lvector[],bumpvector[]} ; };;
+vollam[] += vol[1];
+lup[]    += news-1 ;
+lfront[] += news-5 ;
+lright[] += news-4 ;
+lback[]  += news-3 ;
+lleft[]  += news-2 ;
+EndFor
+
+vol[]=Extrude Surface { news-1, {0,0,e/2}} {Layers {{ndive},{1}} ; };;
+voliso[] += vol[1];
+lupi[]    += news-1 ;
+lfronti[]  += news-5 ;
+lrighti[]  += news-4 ;
+lbacki[] += news-3 ;
+llefti[] += news-2 ;
+
+borderRight[] = {};
+borderLeft[] = {};
+
+For i In{0:2*nlam-1}
+borderRight[] += 11+i*16;
+borderLeft[]  += 17+i*16;
+EndFor
+paxetop = 5 + 40*(2*nlam-1);
+EndIf
+
+
+
+If(HOMO)
+ldwn[]+= surflam0;//Watch Out: this is half a lamination
+vol[]= Extrude Surface { surflam0, {0,0,d/2+e}} {Layers {ndivlam_half,1} ; };;
+vollam[] += vol[1];
+lup[] += news-1;
+lfront[]  += news-5 ;
+lright[]  += news-4 ;
+lback[] += news-3 ;
+lleft[] += news-2 ;
+
+lvector[] = {ndivlam};
+bumpvector[] = {1};
+
+For i In {1:nlam-2}
+ldwn[]+= news-1 ;
+vol[] = Extrude Surface { news-1, {0,0,d+e}} {Layers {lvector[],bumpvector[]} ; };;
+vollam[] += vol[1];
+lup[]    += news-1 ;
+lfront[] += news-5 ;
+lright[] += news-4 ;
+lback[]  += news-3 ;
+lleft[]  += news-2 ;
+EndFor
+
+ldwn[]+= news-1 ;
+vol[]=Extrude Surface { news-1, {0,0,d+e/2}} {Layers {{lvector[]},{bumpvector[]}} ; };;
+vollam[] += vol[1];
+lup[]    += news-1 ;
+lfront[] += news-5 ;
+lright[] += news-4 ;
+lback[]  += news-3 ;
+lleft[]  += news-2 ;
+
+borderRight[] = {};
+borderLeft[] = {};
+
+For i In{0:nlam-1}
+borderRight[] += 11+i*16;
+borderLeft[]  += 17+i*16;
+EndFor
+paxetop = 5 + 40*(nlam-1);
+
+EndIf
+
+
+// I N D U C T O R
+
+pi1 = newp ; Point(pi1)= {x0+rind2,      y0,        z0, lcind};
+pi2 = newp ; Point(pi2)= {x0+rind2+wind, y0,        z0, lcind};
+pi3 = newp ; Point(pi3)= {x0+rind2+wind, y0+wind/2, z0, lcind2};
+pi4 = newp ; Point(pi4)= {x0+rind2,      y0+wind/2, z0, lcind2};
+
+li1 = newl ; Line(li1)= {pi1,pi2};
+li2 = newl ; Line(li2)= {pi2,pi3};
+li3 = newl ; Line(li3)= {pi3,pi4};
+li4 = newl ; Line(li4)= {pi4,pi1};
+
+lineloopind = newl;
+Line Loop (lineloopind) = {li1,li2,li3,li4};
+surfind = news;
+Plane Surface(surfind) = {lineloopind};
+
+If (TRANSIND)
+Transfinite Line{li1,li3} = Ceil[wind/lcind] ;
+Transfinite Line{li2,li4} = Ceil[wind/2/lcind] ;
+Transfinite Surface{surfind} = {pi1,pi2,pi3,pi4};
+vol[] = Extrude{{0,1,0},{0,0,0}, -Pi/2}{Surface{surfind}; Layers {{nind},{1}}; };
+volind = vol[1];
+EndIf
+
+If (!TRANSIND)
+vol[] = Extrude{{0,1,0},{0,0,0}, -Pi/2}{Surface{surfind}; };
+volind = vol[1];
+EndIf
+
+Transfinite Line{li1,li3} = Ceil[wind/lcind] ;
+Transfinite Line{li2,li4,li2+4,li4+4} = Ceil[wind/2/lcind]+1 ;
+
+symind   = surfind + 2 ;
+skinind1 = surfind + 3 ;
+skinind2 = surfind + 4 ;
+skinind3 = surfind + 5 ;
+surfind_ = surfind + 6 ;
+
+// Air
+pb0 = newp ; Point(pb0)={Rint,y0,z0,lcboxi};
+pb1 = newp ; Point(pb1)={0,y0,Rint,lcboxi};
+pb2 = newp ; Point(pb2)={0,Rint,z0,lcbox};
+
+// center p0
+cs1 = newl; Circle(cs1)={pb0,p0,pb1};
+cs2 = newl; Circle(cs2)={pb1,p0,pb2};
+cs3 = newl; Circle(cs3)={pb0,p0,pb2};
+
+lcutr1 = newl ; Line(lcutr1) = {p1,pi1};
+lcutr2 = newl ; Line(lcutr2) = {pi2,pb0};
+laxe0 = newl ;  Line(laxe0)  = {paxetop,pi1+4};
+laxe1 = newl ;  Line(laxe1)  = {pi2+4,pb1};
+lcutl = newl ;  Line(lcutl)  = {p3,pb2};
+
+pb0_ = newp ; Point(pb0_)={Rext,y0,z0,lcboxout};
+pb1_ = newp ; Point(pb1_)={0,y0,Rext,lcboxout};
+pb2_ = newp ; Point(pb2_)={0,Rext,z0,lcboxout};
+
+// center p0
+cs1_ = newl; Circle(cs1_)={pb0_,p0,pb1_};
+cs2_ = newl; Circle(cs2_)={pb1_,p0,pb2_};
+cs3_ = newl; Circle(cs3_)={pb0_,p0,pb2_};
+
+lairinf1 = newl ; Line(lairinf1) = {pb0,pb0_};
+lairinf2 = newl ; Line(lairinf2) = {pb1,pb1_};
+lairinf3 = newl ; Line(lairinf3) = {pb2,pb2_};
+
+surfge0_1 = news ;
+Line Loop(surfge0_1) = 
+{lcutr1,-li4,-li3,-li2,lcutr2,cs3,-lcutl,-l2,-l1};
+Plane Surface(surfge0_1) = {surfge0_1};
+
+surfge0_inf1 = news ;
+Line Loop(surfge0_inf1) = {lairinf1,cs3_,-lairinf3,-cs3};
+Plane Surface(surfge0_inf1) = {surfge0_inf1};
+
+surfge0_2 = news ;
+Line Loop(surfge0_2) = 
+{laxe0,-(li4+4),-(li3+4),-(li2+4),laxe1,cs2,-lcutl,
+borderLeft[],(borderLeft[#borderLeft[]-1]-9)};
+Plane Surface(surfge0_2) = {surfge0_2};
+
+surfge0_inf2 = news ;
+Line Loop(surfge0_inf2) = {lairinf3,-cs2_,-lairinf2,cs2};
+Plane Surface(surfge0_inf2) = {surfge0_inf2};
+
+
+surfgh0_1 = news ;
+Line Loop(surfgh0_1) = 
+{borderRight[],-(borderRight[#borderRight[]-1]-6),laxe0,-(li4+6),-lcutr1};
+Plane Surface(surfgh0_1) = {surfgh0_1};
+
+surfgh0_2 = news ;
+Line Loop(surfgh0_2) = {lcutr2, cs1,-laxe1,-(li1+10)};
+Plane Surface(surfgh0_2) = {surfgh0_2};
+
+surfgh0_inf = news ;
+Line Loop(surfgh0_inf) = {lairinf1,cs1_,-lairinf2,-cs1};
+Plane Surface(surfgh0_inf) = {surfgh0_inf};
+
+surfinf_in = news ;
+Line Loop(surfinf_in) = {cs3,-cs2,-cs1};
+Ruled Surface(surfinf_in) = {surfinf_in};
+
+surfinf_out = news;
+Line Loop(surfinf_out) = {cs1_,cs2_,-cs3_};
+Ruled Surface(surfinf_out) = {surfinf_out};
+
+
+If(!HOMO)
+Surface Loop(newsl) = 
+{skinind3,-surfge0_2,surfgh0_1,surfge0_1,
+ lright[],lrighti[],lback[],lbacki[],lupi[#voliso[]-1],
+-surfinf_in,surfgh0_2,skinind1,skinind2};
+EndIf
+
+If(HOMO)
+Surface Loop(newsl) = 
+{skinind3,-surfge0_2,surfgh0_1,surfge0_1,
+ lright[],lback[],lup[#vollam[]-1],
+-surfinf_in,surfgh0_2,skinind1,skinind2};
+EndIf
+volair = newv;
+Volume(volair) = {newsl-1};
+
+Surface Loop(newsl) = {surfge0_inf2,surfge0_inf1,-surfgh0_inf,surfinf_out,surfinf_in};
+volairinf = newv;
+Volume(volairinf) = {newsl-1};
+
+
+//--------------------------------------------------------
+
+If(!HOMO)
+For i In {0:#vollam[]-1}
+  Physical Volume (LAMINATION+i)= {vollam[i]};
+  If(i==0)
+    //Physical Surface(SKINLAM+i)= {lup[i],lfront[i],lright[i],lback[i],lleft[i]};
+    Physical Surface(SKINLAM+i)= {lup[i],lright[i],lback[i]};
+  EndIf
+  If(i>0)
+    //Physical Surface(SKINLAM+i)= {lup[i],ldwn[i],lfront[i],lright[i],lback[i],lleft[i]};
+    Physical Surface(SKINLAM+i)= {lup[i],ldwn[i],lright[i],lback[i]};
+  EndIf
+EndFor
+Physical Volume(ISOLATION)  = {voliso[]};
+
+Physical Surface(SURFACEGE0) = 
+{surfge0_1,surfge0_2,surfge0_inf1,surfge0_inf2, surflam0, lleft[],llefti[],surfind_,surfind};
+Physical Surface(CUTIND) = {lfront[],lfronti[],surfgh0_1};
+Physical Surface(SURFACEGH0) = {lfront[],lfronti[],symind,surfgh0_1,surfgh0_2,surfgh0_inf};
+EndIf
+
+If(HOMO)
+  Physical Volume (LAMINATION) = {vollam[]};
+  Physical Surface(SKINLAM)    = {lup[#vollam[]-1],lright[],lback[]};
+  Physical Surface(SURFACEGE0) = 
+{surfge0_1,surfge0_2,surfge0_inf1,surfge0_inf2,surflam0,lleft[],surfind_,surfind};
+Physical Surface(CUTIND) = {lfront[],surfgh0_1};
+Physical Surface(SURFACEGH0) = {lfront[],symind,surfgh0_1,surfgh0_2,surfgh0_inf};
+EndIf
+
+Physical Volume(AIR)  = {volair};
+Physical Volume(AIRINF)  = {volairinf};
+
+Physical Volume(IND) = {volind};
+Physical Surface(SKININD) = {skinind1,skinind2,skinind3};
+Physical Surface(ELEC0) = {surfind};
+ 
+Physical Surface(SURFACEGINF) = {surfinf_out};
+
+
+POST = 0;
+If(POST)
+db= 1e-4;
+
+//perpendicular to the laminations
+pp0 = newp ; Point(pp0)={W2/2,db, z0+1e-6,lc};
+pp1 = newp ; Point(pp1)={W2/2,db, z1-1e-6,lc};
+Line(newl)={pp0,pp1};
+
+pp2 = newp ; Point(pp2)={W2-db,W2/2, z0+1e-6,lc};
+pp3 = newp ; Point(pp3)={W2-db,W2/2, z1-1e-6,lc};
+Line(newl)={pp2,pp3};
+
+//paralel to the laminations
+pp0_ = newp ; Point(pp0_)={W2-db, db   , 2*(d+e)+e/2,lc};
+pp1_ = newp ; Point(pp1_)={W2-db, W2-db, 2*(d+e)+e/2,lc};
+Line(newl)={pp0_,pp1_};
+
+pp2_ = newp ; Point(pp2_)={db,    db, 2*(d+e)+e/2,lc};
+pp3_ = newp ; Point(pp3_)={W2-db, db, 2*(d+e)+e/2,lc};
+Line(newl)={pp2_,pp3_};
+
+EndIf