[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"};
};