diff --git a/MagneticForces/InfiniteBox.geo b/MagneticForces/InfiniteBox.geo
index 0d9cb8523269bf9301be47852689a96739e39592..821ab74e6658a51568233e737ae7dc4466cda322 100644
--- a/MagneticForces/InfiniteBox.geo
+++ b/MagneticForces/InfiniteBox.geo
@@ -1,8 +1,6 @@
-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] };
diff --git a/MagneticForces/magnets.geo b/MagneticForces/magnets.geo
index 89143ab0700d57d224c62d711b7d2cae74503a83..622703f5bfa8288c95abd0720ae1fd8d2f4fb044 100644
--- a/MagneticForces/magnets.geo
+++ b/MagneticForces/magnets.geo
@@ -1,94 +1,75 @@
-// 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[] };
+
+
+
+
+
+
+
diff --git a/MagneticForces/magnets_common.pro b/MagneticForces/magnets_common.pro
index 3c75dce75aa16ece8ea75c1329382fc82947e957..49aaef843dde8b5faa5f7c894001356e43ceabb5 100644
--- a/MagneticForces/magnets_common.pro
+++ b/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}