[Gmsh] question about meshing a floating plane inside a cube

Nacho Andres nacho.andres at deos.tudelft.nl
Fri Oct 8 17:36:21 CEST 2004


Hi everybody,

Continuing somehow with the question from Mrs Wang, the only way to
generate a gruyere cheese volume (but with holes over the external
surface) is to declare in precise detail that surface ... and if so, how
to do it correctly for the case of a sphere with many holes over it? So
far I hace tried:

1) Generating first the sphere and then taking the volumes (holes) away
from the volume does not work since coplanarity problems arise when
generating the mesh .. (see lageos.geo, with the rest of the necessary
files)

2) Creating the outer surface first crashes when meshing the sphere with
the spots on it .. besides some of the spots lie over one of the
circumferences that define the sphere ... (see tilted+holes.sph.geo)

I assume that is not an easy problem, and I do not want to take too much
of your time, but I am so close to the solution!!! :P ..  Any help is
highly appreciated!

Many thanks in advance,
Nacho

On Fri, 2004-10-08 at 17:07, Christophe Geuzaine wrote:
> remacle wrote:
> > At 07:25 08/10/2004, C. D. Wang wrote:
> > 
> >> Hi Christophe,
> >>
> >> I am trying to mesh a 3D object (e.g. a cube) where it has a plane
> >> floating inside it. I intended to have elements surrounding the plane
> >> without the plane having to cut through any of them. I tried this, first
> >> in 2D (i.e. a straight line enclosed by a square) with attractor line
> >> and transfinite line, and then 3D with transfinite surface and none
> >> worked. I got elements that were cut through by the plane (or the line).
> >> Can you suggest me any method to do this kind of meshing?
> > 
> > 
> > I'm sorry but gmsh meshing algorithms do not take into account
> > embedded surfaces. 
> 
> More precisely, Gmsh's algorithms are designed in such a way that the
> mesh of an entity is only constrained by the mesh of its boundary--so
> that your floating plane needs to be somehow connected to the
> surrounding volume if you want the mesh of the plane to constrain the
> volume mesh. See
> http://geuz.org/gmsh/doc/texinfo/gmsh_fot.html#FOOT1 for more info.
> 
> Best,
> 
> Christophe
> 

-------------- next part --------------
// This is a extruder generator for getting a tilted sphere
// N at cho Andres, Delft 2004

Function TiltedHoles

// Pre-calculations
     end1= 2.0*Pi/angle1;
     end1= end1-1;
     d1= r*Sin(theta_tilt);
     d2= r*Cos(theta_tilt);
   
// SPHERE

     p1 = newp; Point(p1) = {x,  y,  z_1,  lc1} ;
     p2 = newp; Point(p2) = {x+r,y,  z_1,  lc1} ;
     p3 = newp; Point(p3) = {x,  y+d2,z_1-d1,  lc1} ;
     p4 = newp; Point(p4) = {x,  y+d1,  z_1+d2,lc1} ;
     p5 = newp; Point(p5) = {x-r,y,  z_1,  lc1} ;
     p6 = newp; Point(p6) = {x,  y-d2,z_1+d1,  lc1} ;
     p7 = newp; Point(p7) = {x,  y-d1,  z_1-d2,lc1} ;
   
     c1 = newreg; Circle(c1) = {p2,p1,p7};
     c2 = newreg; Circle(c2) = {p7,p1,p5};
     c3 = newreg; Circle(c3) = {p5,p1,p4};
     c4 = newreg; Circle(c4) = {p4,p1,p2};
     c5 = newreg; Circle(c5) = {p2,p1,p3};
     c6 = newreg; Circle(c6) = {p3,p1,p5};
     c7 = newreg; Circle(c7) = {p5,p1,p6};
     c8 = newreg; Circle(c8) = {p6,p1,p2};
     c9 = newreg; Circle(c9) = {p7,p1,p3};
     c10 = newreg; Circle(c10) = {p3,p1,p4};
     c11 = newreg; Circle(c11) = {p4,p1,p6};
     c12 = newreg; Circle(c12) = {p6,p1,p7};
   
     l1 = newreg; Line Loop(l1) = {c5,c10,c4}; ps1 = newreg;  Ruled Surface(ps1) = {l1};
     l2 = newreg; Line Loop(l2) = {c9,-c5,c1}; ps2 = newreg;  Ruled Surface(ps2) = {l2};
     l3 = newreg; Line Loop(l3) = {-c12,c8,c1}; ps3 = newreg; Ruled Surface(ps3) = {l3};
     l4 = newreg; Line Loop(l4) = {c8,-c4,c11}; ps4 = newreg; Ruled Surface(ps4) = {l4};
     l5 = newreg; Line Loop(l5) = {-c10,c6,c3}; ps5 = newreg; Ruled Surface(ps5) = {l5};
     l6 = newreg; Line Loop(l6) = {-c11,-c3,c7}; ps6 = newreg; Ruled Surface(ps6) = {l6};
     l7 = newreg; Line Loop(l7) = {c2,c7,c12}; ps7 = newreg;  Ruled Surface(ps7) = {l7};
     l8 = newreg; Line Loop(l8) = {-c6,-c9,c2}; ps8 = newreg; Ruled Surface(ps8) = {l8};
  
     theloops[t] = newreg; Surface Loop(theloops[t]) = {l8+1,l5+1,l1+1,l2+1,-(l3+1),-(l7+1),l6+1,l4+1};
  
//  geometrical data (see [Slabinski,1988])
     Theta_CCR[] = {0.000,10.119,19.849,29.579,39.309,49.039,58.769,67.018,76.748,85.135};
     Mu_CCR[] = {0.000,55.577,10.577,5.577,12.658,10.339,8.619,2.812,5.625,0.000};
     n_CCR[] = {1.0,6.0,12.0,18.0,23.0,27.0,31.0,31.0,32.0,32.0};
//   Loop
      i = 0; r = 0; bandmax = 10;
      For i In {0:bandmax-1}
//       Northern Hemisphere
         Theta_c = Theta_CCR[i]*Pi/180.0;
         For j In {0:n_CCR[i]-1}
            Mu_c = Mu_CCR[i]*Pi/180.0+2.0*j*Pi/n_CCR[i];
//         Trigonometrical combinations
            sth = Sin(Theta_c);
            cth = Cos(Theta_c);
            smu = Sin(Mu_c);
            cmu = Cos(Mu_c);
            smsa = smu*Sin(angle1);
            smca = smu*Cos(angle1);
            cmca = cmu*Cos(angle1);
            stsa = sth*Sin(angle1);
            ctcm = cth*cmu;
            ctsm = cth*smu;
            stcm = sth*cmu;
            stsm = sth*smu;
            ctcmsa = cth*cmu*Sin(angle1);
            ctsmsa = cth*smu*Sin(angle1);
            ctsmca = cth*smu*Cos(angle1);
            stcmsa = sth*cmu*Sin(angle1);
            stsmca = sth*smu*Cos(angle1);
//        Points
            p101 = newp; Point(p101) = {stcm*z_2, stsm*z_2, cth*z_2, lc1};
            p102 = newp; Point(p102) = {ctcmsa*ri-smca*ri+stcm*z_2, ctsmsa*ri+cmca*ri+stsm*z_2, -stsa*ri+cth*z_2, lc1};
            p103 = newp; Point(p103) = {-smu*ri+stcm*z_2,cmu*ri+stsm*z_2,cth*z_2,lc1};
//        Circles
            c1 = newreg; Circle(c1) = {p102,p101,p103};
            For x In {1:end1}
               Rotate {{stcm,stsm,cth},{0,0,0},x*angle1} {Duplicata{Line{c1};}}
            EndFor
            circulin[r] = newreg; Line Loop(circulin[r]) = {c1,c1+1,c1+2,c1+3};
            r = r+1;
         EndFor
//       Southern Hemisphere
         Theta_c = Pi-Theta_CCR[i]*Pi/180.0;
         For j In {0:n_CCR[i]-1}
            Mu_c = Mu_CCR[i]*Pi/180.0+2.0*j*Pi/n_CCR[i];
//         Trigonometrical combinations
            sth = Sin(Theta_c);
            cth = Cos(Theta_c);
            smu = Sin(Mu_c);
            cmu = Cos(Mu_c);
            smsa = smu*Sin(angle1);
            smca = smu*Cos(angle1);
            cmca = cmu*Cos(angle1);
            stsa = sth*Sin(angle1);
            ctcm = cth*cmu;
            ctsm = cth*smu;
            stcm = sth*cmu;
            stsm = sth*smu;
            ctcmsa = cth*cmu*Sin(angle1);
            ctsmsa = cth*smu*Sin(angle1);
            ctsmca = cth*smu*Cos(angle1);
            stcmsa = sth*cmu*Sin(angle1);
            stsmca = sth*smu*Cos(angle1);
//        Points
            p101 = newp; Point(p101) = {stcm*z_2, stsm*z_2, cth*z_2, lc1};
            p102 = newp; Point(p102) = {ctcmsa*ri-smca*ri+stcm*z_2, ctsmsa*ri+cmca*ri+stsm*z_2, -stsa*ri+cth*z_2, lc1};
            p103 = newp; Point(p103) = {-smu*ri+stcm*z_2,cmu*ri+stsm*z_2,cth*z_2,lc1};
            //p103 = newp; Point(p103) = {-ctcm*ri+stcm*z_2,-ctsm*ri+stsm*z_2,sth*ri+cth*z_2,lc1};
            p104 = newp; Point(p104) = {stcm*z_3, stsm*z_3, cth*z_3, lc1};
//        Circles
            c1 = newreg; Circle(c1) = {p102,p101,p103};
            For x In {1:end1}
               Rotate {{stcm,stsm,cth},{0,0,0},x*angle1} {Duplicata{Line{c1};}}
            EndFor
            circulin[r] = newreg; Line Loop(circulin[r]) = {c1,c1+1,c1+2,c1+3};
            r = r+1;
         EndFor
      EndFor
      prueba = newreg; Ruled Surface(prueba) = {l8,l5,l1,l2,l3,l7,l6,l4,circulin[]};
      probando = newreg; Physical Surface(probando) = {prueba};
      //surfilltd = newreg; Physical Surface (surfilltd) = {-ps1,ps4,ps5,ps6,circulin[]};
      //surfdark = newreg; Physical Surface (surfdark) = {ps2,ps3,-ps7,ps8,circulin[]};

Return

// Example
// geometrical variables
lc1 = .025;
t = 1;
r_LAGEOS = 0.3; angle1=(Pi/2.0);
x= 0.0; y = 0.0; z_1 = 0.0; re = 0.0238; ri = 0.0193; z_2 = Sqrt(r_LAGEOS*r_LAGEOS-re*re); z_3 = r_LAGEOS; theta_tilt= Pi/6.0; r = r_LAGEOS;
Call TiltedHoles;
-------------- next part --------------
// This is a CCR generator for getting LAGEOS geometry
// N at cho Andres, Delft 2004
//

   Function CCR

// Pre-calculations
     end1= 2.0*Pi/angle1;
     end1= end1-1;
//  Trigonometrical combinations
     sth = Sin(Theta_c);
     cth = Cos(Theta_c);
     smu = Sin(Mu_c);
     cmu = Cos(Mu_c);
     smsa = smu*Sin(angle1);
     smca = smu*Cos(angle1);
     cmca = cmu*Cos(angle1);
     stsa = sth*Sin(angle1);
     ctcm = cth*cmu;
     ctsm = cth*smu;
     stcm = sth*cmu;
     stsm = sth*smu;
     ctcmsa = cth*cmu*Sin(angle1);
     ctsmsa = cth*smu*Sin(angle1);
     ctsmca = cth*smu*Cos(angle1);
     stcmsa = sth*cmu*Sin(angle1);
     stsmca = sth*smu*Cos(angle1);

// LOWER CYLINDER
// Points
     // p1 = newp; Point(p1) = {0,0,0, lc1};
     // p100 = newp; Point(p100) = {x_ref,y_ref,z_ref, lc1};
     p101 = newp; Point(p101) = {stcm*z, stsm*z, cth*z, lc10};
     p102 = newp; Point(p102) = {ctcmsa*ri-smca*ri+stcm*z, ctsmsa*ri+cmca*ri+stsm*z, -stsa*ri+cth*z, lc20};
     p103 = newp; Point(p103) = {-smu*ri+stcm*z,cmu*ri+stsm*z,cth*z,lc20};

// Circles
     c1 = newreg; Circle(c1) = {p102,p101,p103};

// Points
     p104 = newp; Point(p104) = {stcm*z_2, stsm*z_2, cth*z_2, lc10};
     p105 = newp; Point(p105) = {ctcmsa*ri-smca*ri+stcm*z_2, ctsmsa*ri+cmca*ri+stsm*z_2, -stsa*ri+cth*z_2, lc20};
     p106 = newp; Point(p106) = {-smu*ri+stcm*z_2,cmu*ri+stsm*z_2,cth*z_2,lc20};

// Lines
     l1 = newreg; Line(l1) = {p101,p102};
     l2 = newreg; Line(l2) = {p103,p101};

// Plane Surface generator
     l3 = newreg; Line Loop(l3) = {-l2,-c1,-l1}; ps101 = newreg;  Ruled Surface(ps101) = {l3};
     // Printf("Ruled surface has a number %g",ps101);

     For x In {1:end1}
       Rotate {{stcm,stsm,cth},{0,0,0},x*angle1} {Duplicata{Surface{ps101};}}
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     // lowScylCCRLow[i,j] = newreg;  Physical Surface(lowScylCCRLow[i,j]) = {ps101,ps101+1,ps101+4,ps101+7};

// Lines
     l4 = newreg; Line(l4) = {p103,p106};
     l5 = newreg; Line(l5) = {p105,p102};

// Circles
     c2 = newreg; Circle(c2) = {p106,p104,p105};

// Ruled Surface generator
     l6 = newreg; Line Loop(l6) = {l4,c2,l5,c1}; ps103 = newreg;  Ruled Surface(ps103) = {l6};
     // Printf("Ruled surface has a number %g",ps103);

     For x In {1:end1}
       Rotate {{stcm,stsm,cth},{0,0,0},x*angle1} {Duplicata{Surface{ps103};}}
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     // latScylCCRLow[i,j] = newreg; Physical Surface (latScylCCRLow[i,j]) = {ps103,ps103+1,ps103+4,ps103+7};

// MID CYLINDER

// Points
     p107 = newp; Point(p107) = {ctcmsa*re-smca*re+stcm*z_2, ctsmsa*re+cmca*re+stsm*z_2, -stsa*re+cth*z_2, lc30};
     p108 = newp; Point(p108) = {-smu*re+stcm*z_2,cmu*re+stsm*z_2,cth*z_2,lc30};

// Circles
     c3 = newreg; Circle(c3) = {p107,p104,p108};

// Lines
     l7 = newreg; Line (l7) = {p108,p106};
     l8 = newreg; Line (l8) = {p105,p107};

// Ruled Surface generator
     l9 = newreg; Line Loop(l9) = {-l8,-c2,-l7,-c3}; ps104 = newreg;  Ruled Surface(ps104) = {l9};
     // Printf("Ruled surface has a number %g",ps104);

     For x In {1:end1}
       Rotate {{stcm,stsm,cth},{0,0,0},x*angle1} {Duplicata{Surface{ps104};}}
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor
     // lowScylCCRM[i,j] = newreg; Physical Surface (lowScylCCRM[i,j]) = {ps104,ps104+1,ps104+6,ps104+11};

// Points
     p109 = newp; Point(p109) = {stcm*z_3, stsm*z_3, cth*z_3, lc10};
     p110 = newp; Point(p110) = {ctcmsa*re-smca*re+stcm*z_3, ctsmsa*re+cmca*re+stsm*z_3, -stsa*re+cth*z_3, lc20};
     p111 = newp; Point(p111) = {-smu*re+stcm*z_3,cmu*re+stsm*z_3,cth*z_3,lc20};

// Lines
     l10 = newreg; Line(l10) = {p108,p111};
     l11 = newreg; Line(l11) = {p110,p107};

// Circles
     c4 = newreg; Circle(c4) = {p111,p109,p110};

// Ruled Surface generator
     l12 = newreg; Line Loop(l12) = {l10,c4,l11,c3}; ps105 = newreg;  Ruled Surface(ps105) = {l12};
     // Printf("Ruled surface has a number %g",ps105);

     For x In {1:end1}
       Rotate {{stcm,stsm,cth},{0,0,0},x*angle1} { Duplicata{ Surface{ps105}; } }
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     // latScylCCRM[i,j] = newreg; Physical Surface(latScylCCRM[i,j]) = {ps105,ps105+1,ps105+4,ps105+7};

// Upper Spherical Shell

// Points
     p0 = newp; Point(p0) = {0,0,0,lc1};
     p112 = newp; Point(p112) = {stcm*z_4, stsm*z_4, cth*z_4, lc10};
     // Printf("The value of point p112 is %g",p112);

// Circles
     c5 = newreg; Circle(c5) = {p110,p0,p112};
     c6 = newreg; Circle(c6) = {p112,p0,p111};

// Ruled Surface generator
     l13 = newreg; Line Loop(l13) = {-c5,-c6,-c4}; ps106 = newreg;  Ruled Surface(ps106) = {l13};
     // Printf("Ruled surface has a number %g",ps106);

     For x In {1:end1}
       Rotate {{stcm,stsm,cth},{0,0,0},x*angle1} {Duplicata{Surface{ps106};}}
       // Printf("Rotated surface SHELL %g, has number %g",x,newreg-1);
     EndFor

     // upSphShell[i,j] = newreg; Physical Surface(upSphShell[i,j]) = {ps106,ps106+1,ps106+5,ps106+9};

// Volume
     SCCR[t] = newreg; Surface Loop(SCCR[t]) = {ps106+5,ps106+1,ps106,ps106+9,ps105+7,ps105,ps105+1,ps105+4,ps104+6,ps104+1,ps104,ps104+11,ps103+7,ps103,ps103+1,ps103+4,ps101+4,ps101+7,ps101,ps101+1}; 
     VCCR[t] = newreg; Volume(VCCR[t]) = SCCR[t];
     CCRname[t] = newreg; Physical Volume (CCRname[t]) = VCCR[t] ;
     t = t+1;

Return

////  geometrical data (see [Slabinski,1988])
//     // Theta_CCR[] = {0.000,10.119,19.849,29.579,39.309,49.039,58.769,67.018,76.748,85.135};
//     // Mu_CCR[] = {0.000,55.577,10.577,5.577,12.658,10.339,8.619,2.812,5.625,0.000};
//     // n_CCR[] = {1.0,6.0,12.0,18.0,23.0,27.0,31.0,31.0,32.0,32.0};
//     Theta_CCR[] = {0.000,10.119};
//     Mu_CCR[] = {0.000,55.577};
//     n_CCR[] = {1.0,6.0};
//
//     lc1 = 0.005;
//     lc10 = 0.025;
//     lc20 = 0.015;
//     lc30 = 0.005;
//     p1 = newp; Point(p1) = {0,0,0, lc1};
//     Printf("CCR");
//     x_ref = 0.0; y_ref = 0.0; z_ref = 0.0;
//     ri = 0.0193; re = 0.0238; height1 = 0.021; height2 = 0.010; angle1=(Pi/2.0); r_LAGEOS = 0.3;
//     z = Sqrt(r_LAGEOS*r_LAGEOS-re*re)-height1-height2; height3 = r_LAGEOS-z-height1-height2;
//     z_2= z+height1; z_3= z+height1+height2; z_4 = r_LAGEOS; // z_4 = z+height1+height2+height3;
//     i = 0; t = 0; s = 1000;
//     For i In {0:1}
////      Northern Hemisphere
//        Theta_c = Theta_CCR[i]*Pi/180.0;
//        For j In {0:n_CCR[i]-1}
//           t = t+1;
//           s = s+1;
//           Mu_c = Mu_CCR[i]*Pi/180.0+2.0*j*Pi/n_CCR[i];
//           Call CCR;
//           CCRname[t] = s;
//           Physical Volume (CCRname[t]) = VCCR[t] ;
//        EndFor
////      Southern Hemisphere
//        Theta_c = Pi-Theta_CCR[i]*Pi/180.0;
//        For j In {0:n_CCR[i]-1}
//           t = t+1;
//           s = s+1;
//           Mu_c = Mu_CCR[i]*Pi/180.0+2.0*j*Pi/n_CCR[i];
//           Call CCR;
//           CCRname[t] = s;
//           Physical Volume (CCRname[t]) = VCCR[t] ;
//        EndFor
//        // i = i+1; The For sentence already sums 1 to the variable ..
//     EndFor


// General.Axes = 0;
// General.Trackball = 0;
// General.RotationCenterGravity = 1;
// General.RotationCenterX = stcm*z_4;
// General.RotationCenterY = stsm*z_4;
// General.RotationCenterZ = cth*z_4;
-------------- next part --------------
// This is a cylinder generator for getting LAGEOS geometry
// N at cho Andres, Delft 2003
//

   Function JustOneCylinder

// Pre-calculations
     end1= 2.0*Pi/angle1;
     end1= end1-1;
     end2= 2.0*Pi/angle2;
     end2= end2-1;

// LOWER CYLINDER
// Points
     p1 = newp; Point(p1) = {0, 0, -z, lc1};
     p2 = newp; Point(p2) = {ri*Sin(angle1), ri*Cos(angle1), -z, lc1};
     p3 = newp; Point(p3) = {0, ri, -z, lc1};

// Circles
     c1 = newreg; Circle(c1) = {p2,p1,p3};

// Points
     p4 = newp; Point(p4) = {0, 0, -z+height1, lc1};
     p5 = newp; Point(p5) = {ri*Sin(angle1), ri*Cos(angle1), -z+height1, lc1};
     p6 = newp; Point(p6) = {0, ri, -z+height1, lc1};

// Lines
     l1 = newreg; Line(l1) = {p1,p2};
     l2 = newreg; Line(l2) = {p3,p1};

// Plane Surface generator
     l3 = newreg; Line Loop(l3) = {-l2,-c1,-l1}; ps1 = newreg;  Ruled Surface(ps1) = {l3};
     // Printf("Ruled surface has a number %g",ps1);

     For x In {1:end1}
       Rotate {{0,0,1.0},{0,0,0},(x*Pi)/10.0} {Duplicata{Surface{ps1};}}
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     l0 = newreg; Line Loop (l0) = {c1,-(c1+7),-(c1+10),-(c1+13),-(c1+16),-(c1+19),-(c1+22),-(c1+25),-(c1+28),-(c1+31),-(c1+34),-(c1+37),-(c1+40),-(c1+43),-(c1+46),-(c1+49),-(c1+52),-(c1+55),-(c1+58),-(c1+61)};

     ps2 = newreg;  Plane Surface(ps2) = {l0};
     lowScylLow = newreg;  Physical Surface(lowScylLow) = {ps2};

// Lines
     l4 = newreg; Line(l4) = {p3,p6};
     l5 = newreg; Line(l5) = {p5,p2};

// Circles
     c2 = newreg; Circle(c2) = {p6,p4,p5};

// Ruled Surface generator
     l6 = newreg; Line Loop(l6) = {l4,c2,l5,c1}; ps3 = newreg;  Ruled Surface(ps3) = {l6};
     // Printf("Ruled surface has a number %g",ps3);

     For x In {1:end1}
       Rotate {{0,0,1.0},{0,0,0},(x*Pi)/10.0} {Duplicata{Surface{ps3};}}
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     latScylLow = newreg; Physical Surface (latScylLow) = {ps3,ps3+1,ps3+4,ps3+7,ps3+10,ps3+13,ps3+16,ps3+19,ps3+22,ps3+25,ps3+28,ps3+31,ps3+34,ps3+37,ps3+40,ps3+43,ps3+46,ps3+49,ps3+52,ps3+55};


// MID CYLINDER

// Points
     p7 = newp; Point(p7) = {re*Sin(angle1), re*Cos(angle1), -z+height1, lc3};
     p8 = newp; Point(p8) = {0, re, -z+height1, lc3};

// Circles
  c3 = newreg; Circle(c3) = {p7,p4,p8};

// Lines
     l7 = newreg; Line (l7) = {p8,p6};
     l8 = newreg; Line (l8) = {p5,p7};

// Ruled Surface generator
     l9 = newreg; Line Loop(l9) = {-l8,-c2,-l7,-c3}; ps4 = newreg;  Ruled Surface(ps4) = {l9};
     // Printf("Ruled surface has a number %g",ps4);

     For x In {1:end1}
       Rotate {{0,0,1.0},{0,0,0},(x*Pi)/10.0} {Duplicata{Surface{ps4};}}
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor
     lowScylM = newreg; Physical Surface (lowScylM) = {ps4,ps4+1,ps4+6,ps4+11,ps4+16,ps4+21,ps4+26,ps4+31,ps4+36,ps4+41,ps4+46,ps4+51,ps4+56,ps4+61,ps4+66,ps4+71,ps4+76,ps4+81,ps4+86,ps4+91};

// Points
     p9 = newp; Point(p9) = {0, 0, -z+height1+height2, lc1};
     p10 = newp; Point(p10) = {ri*Sin(angle1), ri*Cos(angle1), -z+height1+height2, lc2};
     p11 = newp; Point(p11) = {0, ri, -z+height1+height2, lc2};
     p12 = newp; Point(p12) = {re*Sin(angle1), re*Cos(angle1), -z+height1+height2, lc3};
     p13 = newp; Point(p13) = {0, re, -z+height1+height2, lc3};

// Lines
     l10 = newreg; Line(l10) = {p8,p13};
     l11 = newreg; Line(l11) = {p12,p7};

// Circles
     c4 = newreg; Circle(c4) = {p13,p9,p12};

// Ruled Surface generator
     l12 = newreg; Line Loop(l12) = {l10,c4,l11,c3}; ps5 = newreg;  Ruled Surface(ps5) = {l12};
     // Printf("Ruled surface has a number %g",ps5);

     For x In {1:end1}
       Rotate {{0,0,1.0},{0,0,0},(x*angle1)} { Duplicata{ Surface{ps5}; } }
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     latScylM = newreg; Physical Surface(latScylM) = {ps5,ps5+1,ps5+4,ps5+7,ps5+10,ps5+13,ps5+16,ps5+19,ps5+22,ps5+25,ps5+28,ps5+31,ps5+34,ps5+37,ps5+40,ps5+43,ps5+46,ps5+49,ps5+52,ps5+55};

// Lines
     l13 = newreg; Line(l13) = {p11,p13};
     l14 = newreg; Line(l14) = {p12,p10};

// Circles
     c5 = newreg; Circle(c5) = {p10,p9,p11};

// Ruled Surface generator
     l15 = newreg; Line Loop(l15) = {-l13,-c5,-l14,-c4}; ps6 = newreg;  Ruled Surface(ps6) = {l15};
     // Printf("Ruled surface has a number %g",ps6);

     For x In {1:end1}
       Rotate {{0,0,1.0},{0,0,0},(x*angle1)} { Duplicata{ Surface{ps6}; } }
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     upScylM = newreg; Physical Surface(upScylM) = {ps6,ps6+1,ps6+4,ps6+7,ps6+10,ps6+13,ps6+16,ps6+19,ps6+22,ps6+25,ps6+28,ps6+31,ps6+34,ps6+37,ps6+40,ps6+43,ps6+46,ps6+49,ps6+52,ps6+55};

// UPPER CYLINDER

// Points
     p14 = newp; Point(p14) = {0, 0, -z+height1+height2+height3, lc1};                                                                              p15 = newp; Point(p15) = {ri*Sin(angle1), ri*Cos(angle1), -z+height1+height2+height3, lc2};
     p16 = newp; Point(p16) = {0, ri, -z+height1+height2+height3, lc2};

// Lines
     l16 = newreg; Line(l16) = {p11,p16};
     l17 = newreg; Line(l17) = {p15,p10};

// Circles
     c6 = newreg; Circle(c6) = {p16,p14,p15};

// Ruled Surface generator
     l18 = newreg; Line Loop(l18) = {l16,c6,l17,c5}; ps7 = newreg;  Ruled Surface(ps7) = {l18};
     // Printf("Ruled surface has a number %g",ps7);

     For x In {1:19}
       Rotate {{0,0,1.0},{0,0,0},(x*angle1)} { Duplicata{ Surface{ps7}; } }
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     latScylU = newreg; Physical Surface(latScylU) = {ps7,ps7+1,ps7+4,ps7+7,ps7+10,ps7+13,ps7+16,ps7+19,ps7+22,ps7+25,ps7+28,ps7+31,ps7+34,ps7+37,ps7+40,ps7+43,ps7+46,ps7+49,ps7+52,ps7+55};

// Lines
     l19 = newreg; Line (l19) = {p14,p16};
     l20 = newreg; Line (l20) = {p15,p14};

// Ruled Surface generator
     l21 = newreg; Line Loop(l21) = {-l20,-c6,-l19}; ps8 = newreg;  Ruled Surface(ps8) = {l21};
     // Printf("Ruled surface has a number %g",ps8);

     For x In {1:end1}
       Rotate {{0,0,1.0},{0,0,0},(x*angle1)} {Duplicata{Surface{ps8};}}
       // Printf("Rotated surface %g, has number %g",x,newreg-1);
     EndFor

     l30 = newreg; Line Loop (l30) = {c6,c6+5,(c6+8),(c6+11),(c6+14),(c6+17),(c6+20),(c6+23),(c6+26),(c6+29),(c6+32),(c6+35),(c6+38),(c6+41),(c6+44),(c6+47),(c6+50),(c6+53),(c6+56),(c6+59)};

     ps9 = newreg;  Plane Surface(ps9) = {l30};
     upScylU = newreg; Physical Surface (upScylU) = {ps9};

Return

// Example

// lc1 = 0.025;
// lc2 = 0.015;
// lc3 = 0.005;
// rimin = 0.007317; rimax = 0.161; zcylext = 0.1353; zcylint = 0.2488; heightl = 0.1135; heightm = 0.2706; heightu = 0.1135;
// ri = rimin; re= rimax; angle1=(Pi/10.0); angle2=(Pi/18.0); z = zcylint; height1 = heightl; height2 = heightm; height3 = heightu;
// Call JustOneCylinder;
-------------- next part --------------
// This is a extruder generator for getting a sphere
// N at cho Andres, Delft 2003

Function Sphere

  p1 = newp; Point(p1) = {x,  y,  z,  lc1} ;
  p2 = newp; Point(p2) = {x+r,y,  z,  lc1} ;
  p3 = newp; Point(p3) = {x,  y+r,z,  lc1} ;
  p4 = newp; Point(p4) = {x,  y,  z+r,lc1} ;
  p5 = newp; Point(p5) = {x-r,y,  z,  lc1} ;
  p6 = newp; Point(p6) = {x,  y-r,z,  lc1} ;
  p7 = newp; Point(p7) = {x,  y,  z-r,lc1} ;

  c1 = newreg; Circle(c1) = {p2,p1,p7};
  c2 = newreg; Circle(c2) = {p7,p1,p5};
  c3 = newreg; Circle(c3) = {p5,p1,p4};
  c4 = newreg; Circle(c4) = {p4,p1,p2};
  c5 = newreg; Circle(c5) = {p2,p1,p3};
  c6 = newreg; Circle(c6) = {p3,p1,p5};
  c7 = newreg; Circle(c7) = {p5,p1,p6};
  c8 = newreg; Circle(c8) = {p6,p1,p2};
  c9 = newreg; Circle(c9) = {p7,p1,p3};
  c10 = newreg; Circle(c10) = {p3,p1,p4};
  c11 = newreg; Circle(c11) = {p4,p1,p6};
  c12 = newreg; Circle(c12) = {p6,p1,p7};

  l1 = newreg; Line Loop(l1) = {c5,c10,c4}; ps1 = newreg;  Ruled Surface(ps1) = {l1};
  l2 = newreg; Line Loop(l2) = {c9,-c5,c1}; ps2 = newreg;  Ruled Surface(ps2) = {l2};
  l3 = newreg; Line Loop(l3) = {-c12,c8,c1}; ps3 = newreg; Ruled Surface(ps3) = {l3};
  l4 = newreg; Line Loop(l4) = {c8,-c4,c11}; ps4 = newreg; Ruled Surface(ps4) = {l4};
  l5 = newreg; Line Loop(l5) = {-c10,c6,c3}; ps5 = newreg; Ruled Surface(ps5) = {l5};
  l6 = newreg; Line Loop(l6) = {-c11,-c3,c7}; ps6 = newreg; Ruled Surface(ps6) = {l6};
  l7 = newreg; Line Loop(l7) = {c2,c7,c12}; ps7 = newreg;  Ruled Surface(ps7) = {l7};
  l8 = newreg; Line Loop(l8) = {-c6,-c9,c2}; ps8 = newreg; Ruled Surface(ps8) = {l8};

// Physical Surface
  If (outerS) 
    surfsphere = newreg; Physical Surface (surfsphere) = {-ps1,ps4,ps5,ps6,ps2,ps3,-ps7,ps8};
    Color Red {Surface{ps1,ps4,ps5,ps6,ps2,ps3,ps7,ps8};}
  EndIf

  theloops[t] = newreg ;

  Surface Loop(theloops[t]) = {l8+1, l5+1, l1+1, l2+1, -(l3+1), -(l7+1),
                               l6+1, l4+1};
// Physical Volume
//  thehole = newreg ;
//  Volume(thehole) = theloops[t] ;
//  core = newreg; Physical Volume (core) = {thehole}; 

// Difference between the inner and the outer sphere
  If (innerS)
     intsurfsphereout = newreg; Physical Surface (intsurfsphereout) = {ps1,-ps4,-ps5,-ps6,-ps2,-ps3,ps7,-ps8};
  EndIf

Return

// Example
// innerS= 1.0;
// outerS= 0.0;
// lc1 = .0001;
// rext = 0.001;
// t = 1;
// x= 0.0; y = 0.0; z = 0.0; r = rext;
// Call Sphere;
-------------- next part --------------
// This is a extruder generator for getting a tilted sphere
// N at cho Andres, Delft 2003

Function TiltedSphere

  d1= r*Sin(theta);
  d2= r*Cos(theta);

// SPHERE

  // In the following commands we use the reserved variable name
  // `newp', which automatically selects a new point number. This
  // number is chosen as the highest current point number, plus
  // one. (Note that, analogously to `newp', the variables `newc',
  // `news', `newv' and `newreg' select the highest number amongst
  // currently defined curves, surfaces, volumes and `any entities
  // other than points', respectively.)

  p1 = newp; Point(p1) = {x,  y,  z,  lc1} ;
  p2 = newp; Point(p2) = {x+r,y,  z,  lc1} ;
  p3 = newp; Point(p3) = {x,  y+d2,z-d1,  lc1} ;
  p4 = newp; Point(p4) = {x,  y+d1,  z+d2,lc1} ;
  p5 = newp; Point(p5) = {x-r,y,  z,  lc1} ;
  p6 = newp; Point(p6) = {x,  y-d2,z+d1,  lc1} ;
  p7 = newp; Point(p7) = {x,  y-d1,  z-d2,lc1} ;

  c1 = newreg; Circle(c1) = {p2,p1,p7};
  c2 = newreg; Circle(c2) = {p7,p1,p5};
  c3 = newreg; Circle(c3) = {p5,p1,p4};
  c4 = newreg; Circle(c4) = {p4,p1,p2};
  c5 = newreg; Circle(c5) = {p2,p1,p3};
  c6 = newreg; Circle(c6) = {p3,p1,p5};
  c7 = newreg; Circle(c7) = {p5,p1,p6};
  c8 = newreg; Circle(c8) = {p6,p1,p2};
  c9 = newreg; Circle(c9) = {p7,p1,p3};
  c10 = newreg; Circle(c10) = {p3,p1,p4};
  c11 = newreg; Circle(c11) = {p4,p1,p6};
  c12 = newreg; Circle(c12) = {p6,p1,p7};

  // We need non-plane surfaces to define the spherical cheese
  // holes. Here we use ruled surfaces, which can have 3 or 4
  // sides:

  l1 = newreg; Line Loop(l1) = {c5,c10,c4}; ps1 = newreg;  Ruled Surface(ps1) = {l1};
  l2 = newreg; Line Loop(l2) = {c9,-c5,c1}; ps2 = newreg;  Ruled Surface(ps2) = {l2};
  l3 = newreg; Line Loop(l3) = {-c12,c8,c1}; ps3 = newreg; Ruled Surface(ps3) = {l3};
  l4 = newreg; Line Loop(l4) = {c8,-c4,c11}; ps4 = newreg; Ruled Surface(ps4) = {l4};
  l5 = newreg; Line Loop(l5) = {-c10,c6,c3}; ps5 = newreg; Ruled Surface(ps5) = {l5};
  l6 = newreg; Line Loop(l6) = {-c11,-c3,c7}; ps6 = newreg; Ruled Surface(ps6) = {l6};
  l7 = newreg; Line Loop(l7) = {c2,c7,c12}; ps7 = newreg;  Ruled Surface(ps7) = {l7};
  l8 = newreg; Line Loop(l8) = {-c6,-c9,c2}; ps8 = newreg; Ruled Surface(ps8) = {l8};

// Physical Surface
//  surfsphere = newreg; Physical Surface (surfsphere) = {-ps1,ps4,ps5,ps6,ps2,ps3,-ps7,ps8};

  // Please note that all surface meshes are generated by projecting a
  // 2D planar mesh onto the surface, and that this method gives nice
  // results only if the surface's curvature is small enough. If not,
  // you will have to cut the surface in pieces.

  // We then use an array of variables to store the surface loops
  // identification numbers for later reference (we will need these to
  // define the final volume):

  theloops[t] = newreg ;

  Surface Loop(theloops[t]) = {l8+1, l5+1, l1+1, l2+1, -(l3+1), -(l7+1),
                               l6+1, l4+1};
// Physical Volume
//   thehole = newreg ;
//   Volume(thehole) = theloops[t] ;
//   esfera = newreg; Physical Volume (esfera) = {thehole};

// SURFACE PLANE
// For Data Visualization
// Not possible!!!! GIVES ERRORS!!! => do it in Data Explorer
// l10 = newreg; Line Loop (l10) = {c1,c2,c3,c4};
// ps10 = newreg; Plane Surface (ps10) = {l10};
// plane = newreg; Physical Surface (plane) = {ps10};

// Physical Surface
//   surfilltd = newreg; Physical Surface (surfilltd) = {-ps1,ps4,ps5,ps6};
//   surfdark = newreg; Physical Surface (surfdark) = {ps2,ps3,-ps7,ps8};
//   Color Orange {Surface{ps1,ps4,ps5,ps6};}
//   Color Gray70 {Surface{ps2,ps3,ps7,ps8};}

Return

// Example
//lc1 = .025;
//rext = 0.3;
//t = 1;
//x= 0.0; y = 0.0; z = 0.0; r = rext; theta= Pi/6.0;
//Call TiltedSphere;
-------------- next part --------------
// This program generates the whole LAGEOS geometry using 
// other subprograms
// N at cho Andres, Delft 2004

// Loading the files necessary for the running
      //Include "/home/nacho/3Dmsh/gmsh/functions/just.onecyl.geo";
      //Include "/home/nacho/3Dmsh/gmsh/functions/tiltedsphere.geo";
      //Include "/home/nacho/3Dmsh/gmsh/functions/sphere.geo";
      Include "just.onecyl.geo";
      Include "tiltedsphere.geo";
      Include "sphere.geo";
      Include "ccr.geo";

// General variables
      // lc1 = 0.025;
      // lc2 = 0.020;
      // lc3 = 0.015;
      lc1 = 0.025;
      lc2 = 0.025;
      lc3 = 0.025;

// Geometrical variables
      rext = 0.3;
      rimax = 0.161;
      rimin = 0.007317;
      zcylext = 0.1353;
      zcylint = 0.2488;
      heightl = 0.1135;
      heightm = 0.2706;
      heightu = 0.1135;

// Center of the figure
      innerS=1.0;
      outerS=0.0;
      t=1;
      lc1 = 0.001;        //For the external sphere
      Printf("CENTER");
      x= 0.0; y = 0.0; z = 0.0; r = 0.001;
      Call Sphere;
      ps1000 = newreg; Volume(ps1000) = {theloops[1]};
      center = newreg; Physical Volume (center) = {ps1000};

// Internal + External Cylinders
      lc1 = 0.025;
      Printf("CYLINDERS");
      rimin = 0.007317; rimax = 0.161; zcylext = 0.1353; zcylint = 0.2488; heightl = 0.1135; heightm = 0.2706; heightu = 0.1135;
      ri = rimin; re= rimax; angle1=(Pi/10.0); angle2=(Pi/18.0); z = zcylint; height1 = heightl; height2 = heightm; height3 = heightu;
      Call JustOneCylinder;
      CylSurfLoop = newreg; Surface Loop(CylSurfLoop) = {276,279,282,285,288,291,294,297,300,303,306,309,312,315,318,321,266,267,270,273,336,339,342,345,348,351,354,357,360,363,366,369,372,375,378,381,384,329,330,333,396,399,402,405,408,411,414,417,420,423,426,429,432,435,438,441,444,447,392,393,455,454,527,523,519,515,511,507,503,499,495,491,487,483,479,475,471,467,463,459,176,171,166,165,256,251,246,241,236,231,226,221,216,211,206,201,196,191,186,181,112,115,118,121,124,127,130,133,136,139,142,145,148,151,154,157,102,103,106,109,44,47,50,53,56,59,62,65,68,71,74,77,80,83,86,89,92,37,38,41};
//   BeCu volume
       ps201 = newreg; Volume(ps201) = {CylSurfLoop,theloops[1]};
       BeCu = newreg; Physical Volume (BeCu) = {ps201};
       Color Red {Volume{ps201};}

// Tilted Sphere
      t = 2;
      Printf("TILTED SPHERE");
      x= 0.0; y = 0.0; z = 0.0; r = rext; theta= Pi/6.0;
      // x= 0.0; y = 0.0; z = 0.0; r = rext; theta= 0.0;
      Call TiltedSphere;
      surfilltd = newreg; Physical Surface (surfilltd) = {-ps1,ps4,ps5,ps6};
      surfdark = newreg; Physical Surface (surfdark) = {ps2,ps3,-ps7,ps8};
      Color Orange {Surface{ps1,ps4,ps5,ps6};}
      Color Gray70 {Surface{ps2,ps3,ps7,ps8};}
      bored = newreg; Volume(bored) = {theloops[2],CylSurfLoop};
//      Alluminum = newreg; Physical Volume (Alluminum) = {ps301};
      Color Gray70 {Volume{bored};}

// CCRs;
//  geometrical data (see [Slabinski,1988])
     //Theta_CCR[] = {0.000,10.119,19.849,29.579,39.309,49.039,58.769,67.018,76.748,85.135};
     //Mu_CCR[] = {0.000,55.577,10.577,5.577,12.658,10.339,8.619,2.812,5.625,0.000};
     //n_CCR[] = {1.0,6.0,12.0,18.0,23.0,27.0,31.0,31.0,32.0,32.0};
     Theta_CCR[] = {0.000,10.119};
     Mu_CCR[] = {0.000,55.577};
     n_CCR[] = {1.0,6.0};
//  geometrical variables
     lc1 = 0.005;
     lc10 = 0.025;
     lc20 = 0.015;
     lc30 = 0.005;
     ri = 0.0193; re = 0.0238; height1 = 0.021; height2 = 0.010; angle1=(Pi/2.0); r_LAGEOS = 0.3;
     z = Sqrt(r_LAGEOS*r_LAGEOS-re*re)-height1-height2; height3 = r_LAGEOS-z-height1-height2;
     z_2= z+height1; z_3= z+height1+height2; z_4 = r_LAGEOS; // z_4 = z+height1+height2+height3;
//  Loop
     i = 0; t = 0;
     For i In {0:#Theta_CCR[]-1}
//      Northern Hemisphere
        Theta_c = Theta_CCR[i]*Pi/180.0;
        For j In {0:n_CCR[i]-1}
           Mu_c = Mu_CCR[i]*Pi/180.0+2.0*j*Pi/n_CCR[i];
           Call CCR;
        EndFor
//      Southern Hemisphere
        Theta_c = Pi-Theta_CCR[i]*Pi/180.0;
        For j In {0:n_CCR[i]-1}
           Mu_c = Mu_CCR[i]*Pi/180.0+2.0*j*Pi/n_CCR[i];
           Call CCR;
        EndFor
     EndFor

// Alluminum volume
     u = newreg; Volume(u) = {theloops[2],CylSurfLoop,SCCR[]};
     Al_bored = newreg; Physical Volume (Al_bored) = {u};

// Chorraditas finales
     View "comments" {
 // 10 pixels from the left and 15 pixels from the top of the graphic
 // window:
     T2(10,15,0){"AS Department, DEOS Institute, nacho at deos.tudelft.nl"};

 // 10 pixels from the left and 10 pixels from the bottom of the
 // graphic window:
     T2(10,-10,0){"Copyright (C) N at cho"};
     T2(2000,-10,0){"Hola"};
};