improved top-down methodology

MagneticForces/InfiniteBox.geo

-Function Cuboid
-// This function receives a vector pnt[] with 8 points
-// and generates 12 lines, 6 line loops, 1 volume,
-// 1 physical volume and one physicale surface
+SetFactory("OpenCASCADE");
 
+Function Cuboid
 l12 = newl; Line(l12) = {pnt[0], pnt[1]};
 l23 = newl; Line(l23) = {pnt[1], pnt[2]};
 l34 = newl; Line(l34) = {pnt[2], pnt[3]};
@@ -30,11 +28,13 @@ s4 = news; Plane Surface(s4) = { ll4 } ;
 s5 = news; Plane Surface(s5) = { ll5 } ;
 s6 = news; Plane Surface(s6) = { ll6 } ;
 
-sl = newsl; Surface Loop(sl) = { -s1, s2, s3, s4, s5, s6 };
+sl = newsl; Surface Loop(sl) = { s1, s2, s3, s4, s5, s6 };
 v = newv; Volume(v) = { sl };
+Printf("Cuboid v=%g", v);
 
 If( Flag_TransfInf )
-//Mesh.Algorithm3D = 4; // (4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)
+  //Mesh.Algorithm3D = 4; 
+  // (4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)
   For num In { l12:l85 } 
     Transfinite Line{ num } = 10;
   EndFor
@@ -52,8 +52,10 @@ Return
 DefineConstant[
   Flag_InfiniteBox = {1, Choices{0,1}, Name "Infinite box/Add infinite box"}
   Flag_TransfInf = {0, Choices{0,1}, Name "Infinite box/Transfinite mesh", Visible 0}
-  ratioInf = {1.5, Name "Infinite box/Ratio ext-int", Visible Flag_InfiniteBox}
-  //ratioLc = {25, Name "Infinite box/ration int-Lc", Visible Flag_InfiniteBox}
+  ratioInf = {2, Name "Infinite box/Ratio ext-int", Visible Flag_InfiniteBox}
+  ratioBox = {1.25, Name "Infinite box/Ratio int-content", Visible Flag_InfiniteBox}
+  ratioLc = {10, Name "Infinite box/Ratio int-Lc", Visible Flag_InfiniteBox}
+
   xInt = {1, Name "Infinite box/xInt", Visible 0}
   yInt = {1, Name "Infinite box/yInt", Visible 0}
   zInt = {1, Name "Infinite box/zInt", Visible 0}
@@ -65,22 +67,36 @@ DefineConstant[
   zCnt = {0, Name "Infinite box/zCenter", Visible 0}
 ];
 
-If(!Flag_InfiniteBox)
-  Abort; // exit
-EndIf
-
 // Compute parameters related to the Infinite box
-
 BoundingBox;
 xCnt = (General.MinX + General.MaxX) / 2.;
 yCnt = (General.MinY + General.MaxY) / 2.;
 zCnt = (General.MinZ + General.MaxZ) / 2.;
-xInt = (General.MaxX - General.MinX)/2. * ratioBox;
-yInt = (General.MaxY - General.MinZ)/2. * ratioBox;
-zInt = (General.MaxZ - General.MinZ)/2. * ratioBox;
-xExt = xInt * ratioInf;
-yExt = yInt * ratioInf;
-zExt = zInt * ratioInf;
+
+xInt = (General.MaxX - General.MinX)/2.;
+yInt = (General.MaxY - General.MinZ)/2.;
+zInt = (General.MaxZ - General.MinZ)/2.;
+// If BoundingBox is empty, replace it with a unit cube.
+// This makes this file callable in isolation, e.g., for testing
+If (xInt == 0 )
+  xInt=1;
+EndIf
+If (yInt == 0 )
+  yInt=1;
+EndIf
+If (zInt == 0 )
+  zInt=1;
+EndIf
+
+// Compute aspect ratio of the Infinite Box
+diagInt = Sqrt[ xInt^2 + yInt^2 + zInt^2 ];
+pp = 1.2; // pp=1 square InfBox, pp>3 Box and content have the same aspect ratio
+xInt *= ratioBox*(diagInt/xInt)^(1/pp); 
+yInt *= ratioBox*(diagInt/yInt)^(1/pp); 
+zInt *= ratioBox*(diagInt/zInt)^(1/pp); 
+
+
+// FIXME is there another (better) way to convey the values calculated here to getDP
 SetNumber("Infinite box/xInt",xInt);
 SetNumber("Infinite box/yInt",yInt);
 SetNumber("Infinite box/zInt",zInt);
@@ -91,9 +107,7 @@ SetNumber("Infinite box/xCenter",xCnt);
 SetNumber("Infinite box/yCenter",yCnt);
 SetNumber("Infinite box/zCenter",zCnt);
 
-lc1inf = lc2;
-lc2inf = lc1inf;
-
+lc1inf = Sqrt[ xInt^2 + yInt^2 + zInt^2] / ratioLc;
 p1 = newp; Point (p1) = {xCnt-xInt, yCnt-yInt, zCnt-zInt, lc1inf};
 p2 = newp; Point (p2) = {xCnt+xInt, yCnt-yInt, zCnt-zInt, lc1inf};
 p3 = newp; Point (p3) = {xCnt+xInt, yCnt+yInt, zCnt-zInt, lc1inf};
@@ -103,47 +117,51 @@ p6 = newp; Point (p6) = {xCnt+xInt, yCnt-yInt, zCnt+zInt, lc1inf};
 p7 = newp; Point (p7) = {xCnt+xInt, yCnt+yInt, zCnt+zInt, lc1inf};
 p8 = newp; Point (p8) = {xCnt-xInt, yCnt+yInt, zCnt+zInt, lc1inf};
 
-pp1 = newp; Point (pp1) = {xCnt-xExt, yCnt-yExt, zCnt-zExt, lc2inf};
-pp2 = newp; Point (pp2) = {xCnt+xExt, yCnt-yExt, zCnt-zExt, lc2inf};
-pp3 = newp; Point (pp3) = {xCnt+xExt, yCnt+yExt, zCnt-zExt, lc2inf};
-pp4 = newp; Point (pp4) = {xCnt-xExt, yCnt+yExt, zCnt-zExt, lc2inf};
-pp5 = newp; Point (pp5) = {xCnt-xExt, yCnt-yExt, zCnt+zExt, lc2inf};
-pp6 = newp; Point (pp6) = {xCnt+xExt, yCnt-yExt, zCnt+zExt, lc2inf};
-pp7 = newp; Point (pp7) = {xCnt+xExt, yCnt+yExt, zCnt+zExt, lc2inf};
-pp8 = newp; Point (pp8) = {xCnt-xExt, yCnt+yExt, zCnt+zExt, lc2inf};
-
-SurfInf[] = {} ;
-SurfExt[] = {} ;
-VolInf[] = {};
 
-pnt[]={p1,p2,p3,p4,pp1,pp2,pp3,pp4};
-Call Cuboid; 
-SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
-
-pnt[]={p5,p6,p7,p8,pp5,pp6,pp7,pp8};
-Call Cuboid;
-SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+If(Flag_InfiniteBox)
+  xExt = xInt * ratioInf;
+  yExt = yInt * ratioInf;
+  zExt = zInt * ratioInf;
+
+  lc2inf = lc1inf;
+  pp1 = newp; Point (pp1) = {xCnt-xExt, yCnt-yExt, zCnt-zExt, lc2inf};
+  pp2 = newp; Point (pp2) = {xCnt+xExt, yCnt-yExt, zCnt-zExt, lc2inf};
+  pp3 = newp; Point (pp3) = {xCnt+xExt, yCnt+yExt, zCnt-zExt, lc2inf};
+  pp4 = newp; Point (pp4) = {xCnt-xExt, yCnt+yExt, zCnt-zExt, lc2inf};
+  pp5 = newp; Point (pp5) = {xCnt-xExt, yCnt-yExt, zCnt+zExt, lc2inf};
+  pp6 = newp; Point (pp6) = {xCnt+xExt, yCnt-yExt, zCnt+zExt, lc2inf};
+  pp7 = newp; Point (pp7) = {xCnt+xExt, yCnt+yExt, zCnt+zExt, lc2inf};
+  pp8 = newp; Point (pp8) = {xCnt-xExt, yCnt+yExt, zCnt+zExt, lc2inf};
+
+  pnt[]={p1,p2,p3,p4,pp1,pp2,pp3,pp4}; Call Cuboid; 
+  VolInf[] = v;
+  // FIXME seems to be forbidden to have the command "VolInf[] = v;" 
+  // on the same line as "Call Cuboid" 
+  pnt[]={p5,p6,p7,p8,pp5,pp6,pp7,pp8}; Call Cuboid; 
+  VolInf[] += v;
+  pnt[]={p1,p2,p6,p5,pp1,pp2,pp6,pp5}; Call Cuboid; 
+  VolInf[] += {v};
+  pnt[]={p3,p4,p8,p7,pp3,pp4,pp8,pp7}; Call Cuboid; 
+  VolInf[] += {v};
+  pnt[]={p2,p3,p7,p6,pp2,pp3,pp7,pp6}; Call Cuboid; 
+  VolInf[] += {v};
+  pnt[]={p4,p1,p5,p8,pp4,pp1,pp5,pp8}; Call Cuboid; 
+  VolInf[] += {v};
+
+  Physical Volume("InfiniteX", 1) = { VolInf[4], VolInf[5] };
+  Physical Volume("InfiniteY", 2) = { VolInf[2], VolInf[3] };
+  Physical Volume("InfiniteZ", 3) = { VolInf[0], VolInf[1] };
+EndIf
 
-pnt[]={p1,p2,p6,p5,pp1,pp2,pp6,pp5};
-Call Cuboid;
-SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+pnt[]={p1,p2,p3,p4,p5,p6,p7,p8}; Call Cuboid;
+InteriorInfBox = v;
 
-pnt[]={p3,p4,p8,p7,pp3,pp4,pp8,pp7};
-Call Cuboid;
-SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+For num In {0:#VolInf()-1}
+   Printf("VolInf %5g", VolInf[num]);
+EndFor
 
-pnt[]={p2,p3,p7,p6,pp2,pp3,pp7,pp6};
-Call Cuboid;
-SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
 
-pnt[]={p4,p1,p5,p8,pp4,pp1,pp5,pp8};
-Call Cuboid;
-SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
 
-Coherence;
 
-Physical Volume("InfiniteX", 1) = { VolInf[4], VolInf[5] };
-Physical Volume("InfiniteY", 2) = { VolInf[2], VolInf[3] };
-Physical Volume("InfiniteZ", 3) = { VolInf[0], VolInf[1] };
 
 

MagneticForces/magnets.geo

-// include parameters common to geometry and solver
+SetFactory("OpenCASCADE");
 
 Include "magnets_common.pro";
 
-// define geometry-specific parameters
-DefineConstant[
-  lc1 = {TotalMemory <= 2048 ? 10*mm : 5*mm, Name "Parameters/2Mesh size on magnets [m]"}
-  lc2 = {TotalMemory <= 2048 ? 20*mm : 10*mm, Name "Parameters/2Mesh size at infinity [m]"}
-  ratioBox = {2.0, Name "Parameters/2Ratio Box to model size"}
-];
+mm = 1.e-3;
 
-// change global Gmsh options
+SetFactory("OpenCASCADE");
+
+
+// set some global Gmsh options
 Mesh.Optimize = 1; // optimize quality of tetrahedra
-Mesh.VolumeEdges = 0; // hide volume edges
-Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
+Mesh.VolumeEdges = 0; // Toggle mesh display
+Mesh.SurfaceEdges = 0; 
 Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)
 
-// create magnets
+VolMagnets[] = {};
 For i In {1:NumMagnets}
+  iLabel = newv;
   If(M~{i} == 0) // cylinder
-    p1 = newp; Point(p1) = {X~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
-    p2 = newp; Point(p2) = {X~{i}+R~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
-    p3 = newp; Point(p3) = {X~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
-    p4 = newp; Point(p4) = {X~{i}+R~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
-    l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p4};
-    l3 = newl; Line(l3) = {p4,p3}; l4 = newl; Line(l4) = {p3,p1};
-    ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
-    s1 = news; Plane Surface(s1) = {ll1};
-    e1[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{s1}; };
-    e2[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e1[0]}; };
-    e3[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e2[0]}; };
-    e4[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e3[0]}; };
-    Magnet~{i}[] = {e1[1], e2[1], e3[1], e4[1]};
-  EndIf
-  If(M~{i} == 1) // parallelepiped
-    p1 = newp; Point(p1) = {X~{i}-Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
-    p2 = newp; Point(p2) = {X~{i}+Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
-    p3 = newp; Point(p3) = {X~{i}+Lx~{i}/2, Y~{i}+Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
-    p4 = newp; Point(p4) = {X~{i}-Lx~{i}/2, Y~{i}+Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
-    l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p3};
-    l3 = newl; Line(l3) = {p3,p4}; l4 = newl; Line(l4) = {p4,p1};
-    ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
-    s1 = news; Plane Surface(s1) = {ll1};
-    e1[] = Extrude {0, 0, Lz~{i}} { Surface{s1}; };
-    Magnet~{i}[] = {e1[1]};
+    Cylinder(iLabel) = { X~{i}, Y~{i}-L~{i}/2, Z~{i}, 0, L~{i}, 0, R~{i} };
+  ElseIf(M~{i} == 1) // parallelepiped
+    Box(iLabel) = { X~{i}-Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, Lx~{i}, Ly~{i}, Lz~{i} };
   EndIf
-  Rotate { {0,0,1}, {X~{i},Y~{i},Z~{i}}, deg*Rz~{i} } { Volume{Magnet~{i}[]}; }
-  Rotate { {0,1,0}, {X~{i},Y~{i},Z~{i}}, deg*Ry~{i} } { Volume{Magnet~{i}[]}; }
-  Rotate { {1,0,0}, {X~{i},Y~{i},Z~{i}}, deg*Rx~{i} } { Volume{Magnet~{i}[]}; }
 
-  Physical Volume(Sprintf("Magnet_%g",i),10*i) = Magnet~{i}[]; // magnet volume
-  skin~{i}[] = CombinedBoundary{ Volume{Magnet~{i}[]}; };
+  Rotate { {0,0,1}, {X~{i},Y~{i},Z~{i}}, deg*Rz~{i} } { Volume{ iLabel }; }
+  Rotate { {0,1,0}, {X~{i},Y~{i},Z~{i}}, deg*Ry~{i} } { Volume{ iLabel }; }
+  Rotate { {1,0,0}, {X~{i},Y~{i},Z~{i}}, deg*Rx~{i} } { Volume{ iLabel }; }
+
+  Physical Volume(Sprintf("Magnet_%g",i),10*i) = { iLabel };
+  skin~{i}[] = CombinedBoundary{ Volume{ iLabel }; };
   Physical Surface(Sprintf("SkinMagnet_%g",i),10*i+1) = -skin~{i}[]; // magnet skin
+  VolMagnets[] += { iLabel };
 EndFor
 
 
-If (Flag_InfiniteBox )
-   Include "InfiniteBox.geo";
-Else
-// create air box around magnets
-  BoundingBox; // recompute model bounding box
-  cx = (General.MinX + General.MaxX) / 2;
-  cy = (General.MinY + General.MaxY) / 2;
-  cz = (General.MinZ + General.MaxZ) / 2;
-  lx = (General.MaxX - General.MinX) * 2 * ratioBox;
-  ly = (General.MaxY - General.MinZ) * ratioBox; 
-  lz = (General.MaxZ - General.MinZ) * 2 * ratioBox;
-  p1 = newp; Point (p1) = {cx-lx/2, cy-ly/2, cz-lz/2, lc2};
-  p2 = newp; Point (p2) = {cx+lx/2, cy-ly/2, cz-lz/2, lc2};
-  l1 = newl; Line(l1) = {p1, p2};
-  e1[] = Extrude {0, ly, 0} { Line{l1}; };
-  e2[] = Extrude {0, 0, lz} { Surface{e1[1]}; };
-  Delete { Volume{e2[1]}; }
-  SurfInt[] = { e1[1],e2[0],e2[2],e2[3],e2[4],e2[5] };
-  SurfExt[] = SurfInt[];
-EndIf
+Include "InfiniteBox.geo";
 
+// The overall dimensions of the model have been calculated in InfiniteBox.geo
+// So we use to characteristic length set for the infinite box for the whole mesh.
+Mesh.CharacteristicLengthMin = lc1inf;
+Mesh.CharacteristicLengthMax = lc1inf;
 
 
-For num In {0:#SurfInt[]-1}
-  Printf("%g %g", num, SurfInt[num]);
-EndFor
+AirBox[] = BooleanDifference{ 
+  Volume{ InteriorInfBox }; Delete;
+}{ 
+  Volume{ VolMagnets[] };
+};
 
-temp = newsl; Surface Loop(temp) = { SurfInt[] };
-vv[] = { temp};
-For i In {1:NumMagnets}
-  temp = newsl; Surface Loop(temp) = { skin~{i}[] };
-  vv[] += temp;
+Physical Volume("AirBox",4) = { AirBox[] };
+
+Volumes[] = { VolInf[], VolMagnets[], AirBox[] };
+
+vv[] = BooleanFragments{
+  Volume { Volumes[] }; Delete; }{};
+
+If( #vv[] > #Volumes[] )
+  Error("Overlapping magnets");
+  Abort;
+EndIf
+
+Printf("Check whether BooleanFragments has preserved volume numbering:");
+For num In {0:#vv()-1}
+   Printf("Fragment %5g -> %g", Volumes[num], vv(num));
 EndFor
-v1 = newv; Volume(v1) = { vv[] };
-Physical Volume("AirBox", 4) = v1;
 
-Physical Surface("OuterSurface", 5) = { SurfExt[] };
-Physical Surface("InnerSurface", 6) = { SurfInt[] };
+Outer[] = CombinedBoundary{ Volume{ Volumes[] }; };
+Physical Surface("OuterSurface", 5) = { Outer[]  };
+
+
+
+
+
+
+

MagneticForces/magnets_common.pro

@@ -4,8 +4,6 @@ DefineConstant[
   NumMagnets = {2, Min 1, Max 20, Step 1, Name "Parameters/0Number of magnets"}
   Flag_InfiniteBox = {1, Choices{0,1}, Name "Infinite box/Add infinite box"}
   Flag_FullMenu = {0, Choices{0,1}, Name "Parameters/Show all parameters"}
-  //DistToBox = {50*mm, Name "Infinite box/Distance to box [m]", Visible Flag_InfiniteBox}
-  //lcInf = {DistToBox, Name "Infinite box/Mesh size", Visible Flag_InfiniteBox}
 ];
 
 For i In {1:NumMagnets}