Skip to content
Snippets Groups Projects
Select Git revision
  • c9b0274822c865135ba4c27184c6d01466db0718
  • master default protected
2 results

InfiniteBox.geo

Blame
  • 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