Select Git revision
SElement.cpp
InfiniteBox.geo 5.85 KiB
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]};
l41 = newl; Line(l41) = {pnt[3], pnt[0]};
l56 = newl; Line(l56) = {pnt[4], pnt[5]};
l67 = newl; Line(l67) = {pnt[5], pnt[6]};
l78 = newl; Line(l78) = {pnt[6], pnt[7]};
l85 = newl; Line(l85) = {pnt[7], pnt[4]};
l15 = newl; Line(l15) = {pnt[0], pnt[4]};
l26 = newl; Line(l26) = {pnt[1], pnt[5]};
l37 = newl; Line(l37) = {pnt[2], pnt[6]};
l48 = newl; Line(l48) = {pnt[3], pnt[7]};
ll1 = newll; Line Loop(ll1) = { l12, l23, l34, l41 };
ll2 = newll; Line Loop(ll2) = { l56, l67, l78, l85 };
ll3 = newll; Line Loop(ll3) = { l12, l26, -l56, -l15 };
ll4 = newll; Line Loop(ll4) = { l23, l37, -l67, -l26 };
ll5 = newll; Line Loop(ll5) = { l34, l48, -l78, -l37 };
ll6 = newll; Line Loop(ll6) = { l41, l15, -l85, -l48 };
s1 = news; Plane Surface(s1) = { ll1 } ;
s2 = news; Plane Surface(s2) = { ll2 } ;
s3 = news; Plane Surface(s3) = { ll3 } ;
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 };
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)
For num In { l12:l85 }
Transfinite Line{ num } = 10;
EndFor
For num In { l15:l48 }
Transfinite Line{ num } = 5;
EndFor
For num In { s1:s6 }
Transfinite Surface{ num } ;
EndFor
Transfinite Volume{ v } ;
EndIf
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 = {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}
xExt = {xInt*ratioInf, Name "Infinite box/xExt", Visible 0}
yExt = {yInt*ratioInf, Name "Infinite box/yExt", Visible 0}
zExt = {zInt*ratioInf, Name "Infinite box/zExt", Visible 0}
xCnt = {0, Name "Infinite box/xCenter", Visible 0}
yCnt = {0, Name "Infinite box/yCenter", Visible 0}
zCnt = {0, Name "Infinite box/zCenter", Visible 0}
];
// 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.;
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);
SetNumber("Infinite box/xExt",xExt);
SetNumber("Infinite box/yExt",yExt);
SetNumber("Infinite box/zExt",zExt);
SetNumber("Infinite box/xCenter",xCnt);
SetNumber("Infinite box/yCenter",yCnt);
SetNumber("Infinite box/zCenter",zCnt);
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};
p4 = newp; Point (p4) = {xCnt-xInt, yCnt+yInt, zCnt-zInt, lc1inf};
p5 = newp; Point (p5) = {xCnt-xInt, yCnt-yInt, zCnt+zInt, lc1inf};
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};
VolInf[]={}; // Define empty array in case Flag_InfiniteBox is not active
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,p3,p4,p5,p6,p7,p8}; Call Cuboid;
InteriorInfBox = v;
For num In {0:#VolInf()-1}
Printf("VolInf %5g", VolInf[num]);
EndFor