SetFactory("OpenCASCADE");
Mesh.SurfaceFaces=1;
Mesh.VolumeFaces=1;
Geometry.OldNewReg=0;
Mesh.CharacteristicLengthMax=1;
Mesh.Optimize=1;
Mesh.OptimizeNetgen=1;


//-----------------------------------------
Sphere(0)={-5.76609,1.62854,-7.43851,6.5494};
Sphere(1)={-5.76609,1.62854,17.5615,6.5494};
Sphere(2)={5.46091,-6.09951,-7.02879,3.15384};
Sphere(3)={-11.9781,7.26836,5.27169,3.42841};
Sphere(4)={13.0219,7.26836,5.27169,3.42841};
Sphere(5)={-0.655766,10.7521,10.5308,3.5641};
Sphere(6)={-0.655766,-14.2479,-14.4692,3.5641};
Sphere(7)={-0.655766,10.7521,-14.4692,3.5641};
Sphere(8)={-0.655766,-14.2479,10.5308,3.5641};
Sphere(9)={0.219352,-8.52649,-1.13186,3.14633};
Sphere(10)={4.93294,8.95222,-7.52396,2.9908};
Sphere(11)={2.57157,-0.534455,6.99065,3.32245};
Sphere(12)={-10.4156,-10.5811,-4.79865,3.15397};
Sphere(13)={14.5844,14.4189,-4.79865,3.15397};
Sphere(14)={-10.4156,14.4189,-4.79865,3.15397};
Sphere(15)={14.5844,-10.5811,-4.79865,3.15397};
Sphere(16)={8.98045,-5.4487,-0.159841,2.84704};
Sphere(17)={3.8831,5.96599,-1.79258,1.55282};
Sphere(18)={7.73606,-12.2894,7.63764,2.15519};
Sphere(19)={7.73606,12.7106,7.63764,2.15519};
Sphere(20)={-6.65571,-9.31082,7.86195,2.42687};
Sphere(21)={10.0097,0.240715,2.5998,1.97489};
Sphere(22)={6.77244,-0.5047,-5.27474,1.7112};
Sphere(23)={-9.93568,-4.90776,3.46805,2.06251};
Sphere(24)={-2.79817,2.96525,7.01735,1.79715};
Sphere(25)={-1.73078,5.65928,1.92305,1.67478};
Sphere(26)={-0.378048,-0.857189,0.941877,1.7211};
Sphere(27)={4.19362,1.57819,0.555055,2.29541};
Sphere(28)={-6.60077,9.9517,0.329981,2.02214};
Sphere(29)={-5.24194,1.11354,3.84953,1.81833};
Sphere(30)={9.32455,2.08327,8.06185,1.92495};
Sphere(31)={9.04149,0.841166,-9.83039,1.81429};
Sphere(32)={0.952559,-9.78843,-10.1356,1.53765};
Sphere(33)={9.72282,-3.74119,5.91868,1.68636};
Sphere(34)={-7.81007,-4.43854,9.25207,1.76072};
Sphere(35)={-0.64127,9.67322,3.38717,1.67751};
Sphere(36)={4.72541,8.99266,5.39605,1.89746};
Sphere(37)={-7.4942,10.5392,9.59082,1.55948};
Sphere(38)={-3.21169,-5.01457,7.29278,2.09926};
Sphere(39)={9.47485,4.99321,-1.22951,2.68247};
Sphere(40)={-12.4542,-7.94511,9.22613,2.7072};
Sphere(41)={12.5458,-7.94511,9.22613,2.7072};
Sphere(42)={2.56012,5.88359,2.50824,1.74265};
Sphere(43)={-9.15822,-0.136189,6.5817,1.59692};
Sphere(44)={-6.11553,5.52118,7.8238,1.54842};
Sphere(45)={4.19057,4.67811,11.8797,2.22402};
Sphere(46)={4.19057,4.67811,-13.1203,2.22402};
Sphere(47)={2.01994,9.47562,-2.60361,1.55043};
Sphere(48)={6.20785,-5.39911,9.86167,1.72637};
Sphere(49)={-5.92936,-8.80116,1.38211,1.67242};
Sphere(50)={2.25799,-10.0852,5.53186,2.14357};
Sphere(51)={8.76225,8.09084,12.2589,2.24949};
Sphere(52)={8.76225,8.09084,-12.7411,2.24949};
Sphere(53)={0.35287,-6.1491,11.8896,2.10757};
Sphere(54)={0.35287,-6.1491,-13.1104,2.10757};
Sphere(55)={10.3004,4.59418,-7.77345,1.69173};
Sphere(56)={3.36047,-5.79737,3.1064,1.71049};
Sphere(57)={-2.33734,-4.38818,2.35336,1.63613};
Sphere(58)={9.58243,-10.7704,2.87828,1.51118};
Sphere(59)={-6.37494,5.7783,2.31674,1.64981};
Rotate {{1,0,0},{-0.655766,10.7521,10.5308}, 1.570796} {Volume{5};}
Rotate {{1,0,0},{-0.655766,-14.2479,-14.4692}, 1.570796} {Volume{6};}
Rotate {{1,0,0},{-0.655766,10.7521,-14.4692}, 1.570796} {Volume{7};}
Rotate {{1,0,0},{-0.655766,-14.2479,10.5308}, 1.570796} {Volume{8};}
Rotate {{1,0,0},{-10.4156,-10.5811,-4.79865}, 1.570796} {Volume{12};}
Rotate {{1,0,0},{14.5844,14.4189,-4.79865}, 1.570796} {Volume{13};}
Rotate {{1,0,0},{-10.4156,14.4189,-4.79865}, 1.570796} {Volume{14};}
Rotate {{1,0,0},{14.5844,-10.5811,-4.79865}, 1.570796} {Volume{15};}
Rotate {{1,0,0},{7.73606,-12.2894,7.63764}, 1.570796} {Volume{18};}
Rotate {{1,0,0},{7.73606,12.7106,7.63764}, 1.570796} {Volume{19};}

xmin=-12.5;
xmax=12.5;
ymin=-12.5;
ymax=12.5;
zmin=-12.5;
zmax=12.5;
maxDiameter=16;

maxSphere=59;
For i In {1:6}
	boxNumbers(i-1)=maxSphere+i;
EndFor

Box(boxNumbers(0))={xmin-1.5*maxDiameter,ymin-1.5*maxDiameter,zmin,1.5*maxDiameter,ymax-ymin+2*1.5*maxDiameter,zmax-zmin};
Box(boxNumbers(1))={xmax,ymin-1.5*maxDiameter,zmin,1.5*maxDiameter,ymax-ymin+2*1.5*maxDiameter,zmax-zmin};
Box(boxNumbers(2))={xmin,ymin-1.5*maxDiameter,zmin,xmax-xmin,1.5*maxDiameter,zmax-zmin};
Box(boxNumbers(3))={xmin,ymax,zmin,xmax-xmin,1.5*maxDiameter,zmax-zmin};
Box(boxNumbers(4))={xmin-1.5*maxDiameter,ymin-1.5*maxDiameter,zmin-1.5*maxDiameter,xmax-xmin+2*1.5*maxDiameter,ymax-ymin+2*1.5*maxDiameter,1.5*maxDiameter};
Box(boxNumbers(5))={xmin-1.5*maxDiameter,ymin-1.5*maxDiameter,zmax,xmax-xmin+2*1.5*maxDiameter,ymax-ymin+2*1.5*maxDiameter,1.5*maxDiameter};

tool=BooleanUnion{Volume {boxNumbers(0)};Delete;} {Volume {boxNumbers()};Delete;};
Printf("Boolean Operation: %g",tool);

For t In {0:maxSphere}
	outcome()=BooleanDifference{Volume {t};} {Volume {tool};};
	If (outcome(0)!=t)

		SpheresToBeDeleted()+=t;
		Printf("Boolean Operation succesfull at Sphere %g. Outcome: %g;",t,outcome[0]);
		SphereClass()+=outcome(0);
	Else
		SphereClass()+=t;
	EndIf

EndFor

Recursive Delete{Volume{SpheresToBeDeleted(),tool};}

mainvolume=newv;
Box(mainvolume)={xmin,ymin,zmin,xmax-xmin,ymax-ymin,zmax-zmin};

Physical Volume(1)={SphereClass()};
Physical Volume(2)={mainvolume};

Coherence;
e=1e-3;
Smin() = Surface In BoundingBox{ -12.5-e, -12.5-e, -12.5-e, -12.5+e, 12.5+e, 12.5+e };
For i In {0:#Smin()-1} 
   bb() = BoundingBox Surface { Smin(i) };
   Smax() = Surface In BoundingBox{ bb(0)-e+25, bb(1)-e, bb(2)-e, bb(3)+e+25, bb(4)+e, bb(5)+e };
   For j In {0:#Smax()-1}
       bb2() = BoundingBox Surface { Smax(j) };
       bb2(0)-=25;
       bb2(3)-=25;
      If(Fabs(bb2(0)-bb(0)) < e && Fabs(bb2(1)-bb(1)) < e &&
      Fabs(bb2(2)-bb(2)) < e && Fabs(bb2(3)-bb(3)) < e &&
      Fabs(bb2(4)-bb(4)) < e && Fabs(bb2(5)-bb(5)) < e)
         Periodic Surface {Smax(j)} = {Smin(i)} Translate {25,0,0};
      EndIf
   EndFor
EndFor

Smin() = Surface In BoundingBox{ -12.5-e, -12.5-e, -12.5-e, 12.5+e, -12.5+e, 12.5+e };
For i In {0:#Smin()-1} 
   bb() = BoundingBox Surface { Smin(i) };
   Smax() = Surface In BoundingBox{ bb(0)-e, bb(1)-e+25, bb(2)-e, bb(3)+e, bb(4)+e+25, bb(5)+e };
   For j In {0:#Smax()-1}
       bb2() = BoundingBox Surface { Smax(j) };
       bb2(1)-=25;
       bb2(4)-=25;
      If(Fabs(bb2(0)-bb(0)) < e && Fabs(bb2(1)-bb(1)) < e &&
      Fabs(bb2(2)-bb(2)) < e && Fabs(bb2(3)-bb(3)) < e &&
      Fabs(bb2(4)-bb(4)) < e && Fabs(bb2(5)-bb(5)) < e)
         Periodic Surface {Smax(j)} = {Smin(i)} Translate {0,25,0};
      EndIf
   EndFor
EndFor

Smin() = Surface In BoundingBox{ -12.5-e, -12.5-e, -12.5-e, 12.5+e, 12.5+e, -12.5+e };
For i In {0:#Smin()-1} 
   bb() = BoundingBox Surface { Smin(i) };
   Smax() = Surface In BoundingBox{ bb(0)-e, bb(1)-e, bb(2)-e+25, bb(3)+e, bb(4)+e, bb(5)+e+25 };
   For j In {0:#Smax()-1}
       bb2() = BoundingBox Surface { Smax(j) };
       bb2(2)-=25;
       bb2(5)-=25;
      If(Fabs(bb2(0)-bb(0)) < e && Fabs(bb2(1)-bb(1)) < e &&
      Fabs(bb2(2)-bb(2)) < e && Fabs(bb2(3)-bb(3)) < e &&
      Fabs(bb2(4)-bb(4)) < e && Fabs(bb2(5)-bb(5)) < e)
         Periodic Surface {Smax(j)} = {Smin(i)} Translate {0,0,25};
      EndIf
   EndFor
EndFor


Mesh 3;