[Gmsh] volume mesh
Martin Genet
mgenet at lbl.gov
Thu Mar 29 18:27:31 CEST 2012
Dear GMSH people:
I am having troubles meshing a rather simple geometry: some "bricks" with some "mortar" in between (see .geo files).
If I do not actually mesh the mortar volume, everything seems to be fine: the bricks are correctly meshed (see bricks.png), as well as the mortar surface loop (see mortar_geometry.png) and its orientation (see mortar_orientation.png).
However, when I ask GMSH to actually mesh the mortar volume, when using the frontal algorithm I get a lot of error like
>Error : ERROR: Edge 8 - 2061 multiple times in surface mesh
and then
> Error : ERROR: Surface mesh not consistent
> Error : ERROR: Stop meshing since surface mesh not consistent
and when using the Delaunay algorithm (I use the svn version with netgen) after more than three hours I get a pretty weird mesh: bricks are somehow broken (see bricks_broken.png) and the mortar does not quite fit into its surface mesh (see mortar_broken.png).
I have hard times figuring out what did I do wrong. Does this problem sound familiar to anybody? Many thanks in advance for your help.
Martin
-------------- next part --------------
// ---------
// ------------------------------------------------------- // OPTIONS //
// ---------
General.Terminal = 1;
Mesh.Algorithm = 5; // 2D mesh algorithm (1=MeshAdapt, 5=Delaunay, 6=Frontal)
Mesh.Algorithm3D = 1; // 3D mesh algorithm (1=Delaunay , 4=Frontal)
Mesh.Optimize = 0;
Mesh.OptimizeNetgen = 0;
Include "../TOOLS/Mesh/make_cube.geo";
// ------------
// ---------------------------------------------------- // PARAMETERS //
// ------------
lx = 90.0;
ly = 60.0;
lz = 30.0;
e = 10.0;
n = 10.0;
m = 10.0;
lm = e/n;
lM = (m-1)*lm;
Lx = lx+e;
Ly = ly+e;
Lz = lz+e;
// ------
// ---------------------------------------------------------- // MESH //
// ------
num_F = 1;
num_SI = 0;
// ----------------------------------------------- // brick 1-1 //
xmin = 0;
xmax = lx/2;
ymin = 0;
ymax = ly/2;
zmin = 0;
zmax = lz/2;
cube_type = 1;
Call make_heterogeneous_cube;
Physical Volume ("brick1-1") = {V};
B11_P = P[];
B11_L = L[];
B11_S = S[];
M_SI[num_SI] = S[1]; num_SI++;
M_SI[num_SI] = S[3]; num_SI++;
M_SI[num_SI] = S[5]; num_SI++;
B1_F[0] = num_F-1;
// ----------------------------------------------- // brick 1-2 //
xmin = Lx-lx/2;
xmax = Lx;
ymin = 0;
ymax = ly/2;
zmin = 0;
zmax = lz/2;
cube_type = 2;
Call make_heterogeneous_cube;
Physical Volume ("brick1-2") = {V};
B12_P = P[];
B12_L = L[];
B12_S = S[];
M_SI[num_SI] = S[1]; num_SI++;
M_SI[num_SI] = S[3]; num_SI++;
M_SI[num_SI] = -S[4]; num_SI++;
B1_F[1] = num_F-1;
// ----------------------------------------------- // brick 1-3 //
xmin = 0;
xmax = lx/2;
ymin = Ly-ly/2;
ymax = Ly;
zmin = 0;
zmax = lz/2;
cube_type = 3;
Call make_heterogeneous_cube;
Physical Volume ("brick1-3") = {V};
B13_P = P[];
B13_L = L[];
B13_S = S[];
M_SI[num_SI] = S[1]; num_SI++;
M_SI[num_SI] = -S[2]; num_SI++;
M_SI[num_SI] = S[5]; num_SI++;
B1_F[2] = num_F-1;
// ----------------------------------------------- // brick 1-4 //
xmin = Lx-lx/2;
xmax = Lx;
ymin = Ly-ly/2;
ymax = Ly;
zmin = 0;
zmax = lz/2;
cube_type = 4;
Call make_heterogeneous_cube;
Physical Volume ("brick1-4") = {V};
B14_P = P[];
B14_L = L[];
B14_S = S[];
M_SI[num_SI] = S[1]; num_SI++;
M_SI[num_SI] = -S[2]; num_SI++;
M_SI[num_SI] = -S[4]; num_SI++;
B1_F[3] = num_F-1;
// ----------------------------------------------- // brick 1-5 //
xmin = 0;
xmax = lx/2;
ymin = 0;
ymax = ly/2;
zmin = Lz-lz/2;
zmax = Lz;
cube_type = 5;
Call make_heterogeneous_cube;
Physical Volume ("brick1-5") = {V};
B15_P = P[];
B15_L = L[];
B15_S = S[];
M_SI[num_SI] = -S[0]; num_SI++;
M_SI[num_SI] = S[3]; num_SI++;
M_SI[num_SI] = S[5]; num_SI++;
B1_F[4] = num_F-1;
// ----------------------------------------------- // brick 1-6 //
xmin = Lx-lx/2;
xmax = Lx;
ymin = 0;
ymax = ly/2;
zmin = Lz-lz/2;
zmax = Lz;
cube_type = 6;
Call make_heterogeneous_cube;
Physical Volume ("brick1-6") = {V};
B16_P = P[];
B16_L = L[];
B16_S = S[];
M_SI[num_SI] = -S[0]; num_SI++;
M_SI[num_SI] = S[3]; num_SI++;
M_SI[num_SI] = -S[4]; num_SI++;
B1_F[5] = num_F-1;
// ----------------------------------------------- // brick 1-7 //
xmin = 0;
xmax = lx/2;
ymin = Ly-ly/2;
ymax = Ly;
zmin = Lz-lz/2;
zmax = Lz;
cube_type = 7;
Call make_heterogeneous_cube;
Physical Volume ("brick1-7") = {V};
B17_P = P[];
B17_L = L[];
B17_S = S[];
M_SI[num_SI] = -S[0]; num_SI++;
M_SI[num_SI] = -S[2]; num_SI++;
M_SI[num_SI] = S[5]; num_SI++;
B1_F[6] = num_F-1;
// ----------------------------------------------- // brick 1-8 //
xmin = Lx-lx/2;
xmax = Lx;
ymin = Ly-ly/2;
ymax = Ly;
zmin = Lz-lz/2;
zmax = Lz;
cube_type = 8;
Call make_heterogeneous_cube;
Physical Volume ("brick1-8") = {V};
B18_P = P[];
B18_L = L[];
B18_S = S[];
M_SI[num_SI] = -S[0]; num_SI++;
M_SI[num_SI] = -S[2]; num_SI++;
M_SI[num_SI] = -S[4]; num_SI++;
B1_F[7] = num_F-1;
// -------------------------------------------------- // mortar //
M_L00[0] = B11_L[ 1];
M_L01[0] = -B11_L[ 5];
M_L02[0] = newl; Line (M_L02[0]) = {B11_P[1], B12_P[0]};
M_L03[0] = B12_L[ 4];
M_L04[0] = B12_L[ 1];
M_L05[0] = newl; Line (M_L05[0]) = {B12_P[3], B14_P[1]};
M_L06[0] = -B14_L[ 0];
M_L07[0] = B14_L[ 4];
M_L08[0] = newl; Line (M_L08[0]) = {B14_P[2], B13_P[3]};
M_L09[0] = -B13_L[ 5];
M_L10[0] = -B13_L[ 0];
M_L11[0] = newl; Line (M_L11[0]) = {B13_P[0], B11_P[2]};
M_LL[0] = newll; Line Loop(M_LL[0]) = {M_L00[0], M_L01[0], M_L02[0], M_L03[0], M_L04[0], M_L05[0], M_L06[0], M_L07[0], M_L08[0], M_L09[0], M_L10[0], M_L11[0]};
M_SE[0] = news; Plane Surface (M_SE[0]) = {M_LL[0]};
M_L00[1] = B15_L[ 3];
M_L01[1] = -B15_L[ 7];
M_L02[1] = newl; Line (M_L02[1]) = {B15_P[5], B16_P[4]};
M_L03[1] = B16_L[ 6];
M_L04[1] = B16_L[ 3];
M_L05[1] = newl; Line (M_L05[1]) = {B16_P[7], B18_P[5]};
M_L06[1] = -B18_L[ 2];
M_L07[1] = B18_L[ 6];
M_L08[1] = newl; Line (M_L08[1]) = {B18_P[6], B17_P[7]};
M_L09[1] = -B17_L[ 7];
M_L10[1] = -B17_L[ 2];
M_L11[1] = newl; Line (M_L11[1]) = {B17_P[4], B15_P[6]};
M_LL[1] = newll; Line Loop(M_LL[1]) = {M_L00[1], M_L01[1], M_L02[1], M_L03[1], M_L04[1], M_L05[1], M_L06[1], M_L07[1], M_L08[1], M_L09[1], M_L10[1], M_L11[1]};
M_SE[1] = news; Plane Surface (M_SE[1]) = {M_LL[1]};
M_L00[2] = B11_L[ 9];
M_L01[2] = -B11_L[ 2];
M_L02[2] = newl; Line (M_L02[2]) = {B11_P[4], B15_P[0]};
M_L03[2] = B15_L[ 0];
M_L04[2] = B15_L[ 9];
M_L05[2] = newl; Line (M_L05[2]) = {B15_P[5], B16_P[4]};
M_L06[2] = -B16_L[ 8];
M_L07[2] = B16_L[ 0];
M_L08[2] = newl; Line (M_L08[2]) = {B16_P[1], B12_P[5]};
M_L09[2] = -B12_L[ 2];
M_L10[2] = -B12_L[ 8];
M_L11[2] = newl; Line (M_L11[2]) = {B12_P[0], B11_P[1]};
M_LL[2] = newll; Line Loop(M_LL[2]) = {M_L00[2], M_L01[2], M_L02[2], M_L03[2], M_L04[2], M_L05[2], M_L06[2], M_L07[2], M_L08[2], M_L09[2], M_L10[2], M_L11[2]};
M_SE[2] = news; Plane Surface (M_SE[2]) = {M_LL[2]};
M_L00[3] = B13_L[11];
M_L01[3] = -B13_L[ 3];
M_L02[3] = newl; Line (M_L02[3]) = {B13_P[6], B17_P[2]};
M_L03[3] = B17_L[ 1];
M_L04[3] = B17_L[11];
M_L05[3] = newl; Line (M_L05[3]) = {B17_P[7], B18_P[6]};
M_L06[3] = -B18_L[10];
M_L07[3] = B18_L[ 1];
M_L08[3] = newl; Line (M_L08[3]) = {B18_P[3], B14_P[7]};
M_L09[3] = -B14_L[ 3];
M_L10[3] = -B14_L[10];
M_L11[3] = newl; Line (M_L11[3]) = {B14_P[2], B13_P[3]};
M_LL[3] = newll; Line Loop(M_LL[3]) = {M_L00[3], M_L01[3], M_L02[3], M_L03[3], M_L04[3], M_L05[3], M_L06[3], M_L07[3], M_L08[3], M_L09[3], M_L10[3], M_L11[3]};
M_SE[3] = news; Plane Surface (M_SE[3]) = {M_LL[3]};
M_L00[4] = B11_L[ 6];
M_L01[4] = -B11_L[10];
M_L02[4] = newl; Line (M_L02[4]) = {B11_P[2], B13_P[0]};
M_L03[4] = B13_L[ 8];
M_L04[4] = B13_L[ 6];
M_L05[4] = newl; Line (M_L05[4]) = {B13_P[6], B17_P[2]};
M_L06[4] = -B17_L[ 4];
M_L07[4] = B17_L[ 8];
M_L08[4] = newl; Line (M_L08[4]) = {B17_P[4], B15_P[6]};
M_L09[4] = -B15_L[10];
M_L10[4] = -B15_L[ 4];
M_L11[4] = newl; Line (M_L11[4]) = {B15_P[0], B11_P[4]};
M_LL[4] = newll; Line Loop(M_LL[4]) = {M_L00[4], M_L01[4], M_L02[4], M_L03[4], M_L04[4], M_L05[4], M_L06[4], M_L07[4], M_L08[4], M_L09[4], M_L10[4], M_L11[4]};
M_SE[4] = news; Plane Surface (M_SE[4]) = {M_LL[4]};
M_L00[5] = B12_L[ 7];
M_L01[5] = -B12_L[11];
M_L02[5] = newl; Line (M_L02[5]) = {B12_P[3], B14_P[1]};
M_L03[5] = B14_L[ 9];
M_L04[5] = B14_L[ 7];
M_L05[5] = newl; Line (M_L05[5]) = {B14_P[7], B18_P[3]};
M_L06[5] = -B18_L[ 5];
M_L07[5] = B18_L[ 9];
M_L08[5] = newl; Line (M_L08[5]) = {B18_P[5], B16_P[7]};
M_L09[5] = -B16_L[11];
M_L10[5] = -B16_L[ 5];
M_L11[5] = newl; Line (M_L11[5]) = {B16_P[1], B12_P[5]};
M_LL[5] = newll; Line Loop(M_LL[5]) = {M_L00[5], M_L01[5], M_L02[5], M_L03[5], M_L04[5], M_L05[5], M_L06[5], M_L07[5], M_L08[5], M_L09[5], M_L10[5], M_L11[5]};
M_SE[5] = news; Plane Surface (M_SE[5]) = {M_LL[5]};
// Physical Surface ("mortar_SE") = {M_SE[0], -M_SE[1], M_SE[2], -M_SE[3], M_SE[4], -M_SE[5]};
// Physical Surface ("mortar_SI") = {M_SI[]};
M_SL = newsl; Surface Loop (M_SL) = {M_SE[0], -M_SE[1], M_SE[2], -M_SE[3], M_SE[4], -M_SE[5], M_SI[]};
Physical Surface ("mortar_SL") = {M_SE[0], -M_SE[1], M_SE[2], -M_SE[3], M_SE[4], -M_SE[5], M_SI[]};
M_V = newv; Volume(M_V) = {M_SL};
Physical Volume ("mortar") = {M_V};
// --------------------------------------------- // periodicity //
Periodic Surface M_SE[1] {M_L00[1], M_L01[1], M_L02[1], M_L03[1], M_L04[1], M_L05[1], M_L06[1], M_L07[1], M_L08[1], M_L09[1], M_L10[1], M_L11[1]} = M_SE[0] {M_L00[0], M_L01[0], M_L02[0], M_L03[0], M_L04[0], M_L05[0], M_L06[0], M_L07[0], M_L08[0], M_L09[0], M_L10[0], M_L11[0]};
Periodic Surface B15_S[1] {B15_L[ 2], B15_L[ 7], -B15_L[ 3], -B15_L[ 6]} = B11_S[0] {B11_L[ 0], B11_L[ 5], -B11_L[ 1], -B11_L[ 4]};
Periodic Surface B16_S[1] {B16_L[ 2], B16_L[ 7], -B16_L[ 3], -B16_L[ 6]} = B12_S[0] {B12_L[ 0], B12_L[ 5], -B12_L[ 1], -B12_L[ 4]};
Periodic Surface B17_S[1] {B17_L[ 2], B17_L[ 7], -B17_L[ 3], -B17_L[ 6]} = B13_S[0] {B13_L[ 0], B13_L[ 5], -B13_L[ 1], -B13_L[ 4]};
Periodic Surface B18_S[1] {B18_L[ 2], B18_L[ 7], -B18_L[ 3], -B18_L[ 6]} = B14_S[0] {B14_L[ 0], B14_L[ 5], -B14_L[ 1], -B14_L[ 4]};
Periodic Surface M_SE[3] {M_L00[3], M_L01[3], M_L02[3], M_L03[3], M_L04[3], M_L05[3], M_L06[3], M_L07[3], M_L08[3], M_L09[3], M_L10[3], M_L11[3]} = M_SE[2] {M_L00[2], M_L01[2], M_L02[2], M_L03[2], M_L04[2], M_L05[2], M_L06[2], M_L07[2], M_L08[2], M_L09[2], M_L10[2], M_L11[2]};
Periodic Surface B13_S[3] {B13_L[10], B13_L[ 3], -B13_L[11], -B13_L[ 1]} = B11_S[2] {B11_L[ 8], B11_L[ 2], -B11_L[ 9], -B11_L[ 0]};
Periodic Surface B14_S[3] {B14_L[10], B14_L[ 3], -B14_L[11], -B14_L[ 1]} = B12_S[2] {B12_L[ 8], B12_L[ 2], -B12_L[ 9], -B12_L[ 0]};
Periodic Surface B17_S[3] {B17_L[10], B17_L[ 3], -B17_L[11], -B17_L[ 1]} = B15_S[2] {B15_L[ 8], B15_L[ 2], -B15_L[ 9], -B15_L[ 0]};
Periodic Surface B18_S[3] {B18_L[10], B18_L[ 3], -B18_L[11], -B18_L[ 1]} = B16_S[2] {B16_L[ 8], B16_L[ 2], -B16_L[ 9], -B16_L[ 0]};
Periodic Surface M_SE[5] {M_L00[5], M_L01[5], M_L02[5], M_L03[5], M_L04[5], M_L05[5], M_L06[5], M_L07[5], M_L08[5], M_L09[5], M_L10[5], M_L11[5]} = M_SE[4] {M_L00[4], M_L01[4], M_L02[4], M_L03[4], M_L04[4], M_L05[4], M_L06[4], M_L07[4], M_L08[4], M_L09[4], M_L10[4], M_L11[4]};
Periodic Surface B12_S[5] {B12_L[ 5], B12_L[11], -B12_L[07], -B12_L[ 9]} = B11_S[4] {B11_L[ 4], B11_L[10], -B11_L[ 6], -B11_L[ 8]};
Periodic Surface B14_S[5] {B14_L[ 5], B14_L[11], -B14_L[07], -B14_L[ 9]} = B13_S[4] {B13_L[ 4], B13_L[10], -B13_L[ 6], -B13_L[ 8]};
Periodic Surface B16_S[5] {B16_L[ 5], B16_L[11], -B16_L[07], -B16_L[ 9]} = B15_S[4] {B15_L[ 4], B15_L[10], -B15_L[ 6], -B15_L[ 8]};
Periodic Surface B18_S[5] {B18_L[ 5], B18_L[11], -B18_L[07], -B18_L[ 9]} = B17_S[4] {B17_L[ 4], B17_L[10], -B17_L[ 6], -B17_L[ 8]};
// ------------------------------------------------- // density //
Field[num_F] = MathEval;
Field[num_F].F = Sprintf("%g", lm);
num_F++;
Field[num_F] = Max;
Field[num_F].FieldsList = {num_F-1, B1_F[]};
Background Field = num_F;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
// ------------------------------------------------------- // Scaling //
Mesh.ScalingFactor = 1.0;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: brick_broken.png
Type: image/png
Size: 306305 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20120329/81f45e80/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bricks.png
Type: image/png
Size: 1186177 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20120329/81f45e80/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mortar_broken.png
Type: image/png
Size: 959552 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20120329/81f45e80/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mortar_geometry.png
Type: image/png
Size: 595120 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20120329/81f45e80/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mortar_orientation.png
Type: image/png
Size: 1257019 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20120329/81f45e80/attachment-0004.png>
-------------- next part --------------
Function make_cube
P[0] = newp; Point (P[0]) = {xmin, ymin, zmin, l};
P[1] = newp; Point (P[1]) = {xmax, ymin, zmin, l};
P[2] = newp; Point (P[2]) = {xmin, ymax, zmin, l};
P[3] = newp; Point (P[3]) = {xmax, ymax, zmin, l};
P[4] = newp; Point (P[4]) = {xmin, ymin, zmax, l};
P[5] = newp; Point (P[5]) = {xmax, ymin, zmax, l};
P[6] = newp; Point (P[6]) = {xmin, ymax, zmax, l};
P[7] = newp; Point (P[7]) = {xmax, ymax, zmax, l};
L[ 0] = newl; Line (L[ 0]) = {P[0], P[1]};
L[ 1] = newl; Line (L[ 1]) = {P[2], P[3]};
L[ 2] = newl; Line (L[ 2]) = {P[4], P[5]};
L[ 3] = newl; Line (L[ 3]) = {P[6], P[7]};
L[ 4] = newl; Line (L[ 4]) = {P[0], P[2]};
L[ 5] = newl; Line (L[ 5]) = {P[1], P[3]};
L[ 6] = newl; Line (L[ 6]) = {P[4], P[6]};
L[ 7] = newl; Line (L[ 7]) = {P[5], P[7]};
L[ 8] = newl; Line (L[ 8]) = {P[0], P[4]};
L[ 9] = newl; Line (L[ 9]) = {P[1], P[5]};
L[10] = newl; Line (L[10]) = {P[2], P[6]};
L[11] = newl; Line (L[11]) = {P[3], P[7]};
LL[0] = newll; Line Loop (LL[0]) = {L[ 0], L[ 5], -L[ 1], -L[ 4]};
LL[1] = newll; Line Loop (LL[1]) = {L[ 2], L[ 7], -L[ 3], -L[ 6]};
LL[2] = newll; Line Loop (LL[2]) = {L[ 8], L[ 2], -L[ 9], -L[ 0]};
LL[3] = newll; Line Loop (LL[3]) = {L[10], L[ 3], -L[11], -L[ 1]};
LL[4] = newll; Line Loop (LL[4]) = {L[ 4], L[10], -L[ 6], -L[ 8]};
LL[5] = newll; Line Loop (LL[5]) = {L[ 5], L[11], -L[07], -L[ 9]};
For k In {0:5}
S[k] = news; Plane Surface (S[k]) = {LL[k]};
EndFor
SL = newsl; Surface Loop (SL) = {S[0], -S[1], S[2], -S[3], S[4], -S[5]};
V = newv; Volume (V) = {SL};
Return
Function make_heterogeneous_cube
If (cube_type == 0)
xc = (xmin + xmax)/2;
yc = (ymin + ymax)/2;
zc = (zmin + zmax)/2;
EndIf
If (cube_type == 1)
xc = xmin; yc = ymin; zc = zmin;
EndIf
If (cube_type == 2)
xc = xmax; yc = ymin; zc = zmin;
EndIf
If (cube_type == 3)
xc = xmin; yc = ymax; zc = zmin;
EndIf
If (cube_type == 4)
xc = xmax; yc = ymax; zc = zmin;
EndIf
If (cube_type == 5)
xc = xmin; yc = ymin; zc = zmax;
EndIf
If (cube_type == 6)
xc = xmax; yc = ymin; zc = zmax;
EndIf
If (cube_type == 7)
xc = xmin; yc = ymax; zc = zmax;
EndIf
If (cube_type == 8)
xc = xmax; yc = ymax; zc = zmax;
EndIf
If (cube_type == 0)
dx = (xmax - xmin)/2;
dy = (ymax - ymin)/2;
dz = (zmax - zmin)/2;
EndIf
If (cube_type > 0)
dx = xmax - xmin;
dy = ymax - ymin;
dz = zmax - zmin;
EndIf
Field[num_F] = MathEval;
Field[num_F].F = Sprintf("Fabs(x - %g) / %g", xc, dx);
num_F++;
Field[num_F] = MathEval;
Field[num_F].F = Sprintf("Fabs(y - %g) / %g", yc, dy);
num_F++;
Field[num_F] = MathEval;
Field[num_F].F = Sprintf("Fabs(z - %g) / %g", zc, dz);
num_F++;
Field[num_F] = Max;
Field[num_F].FieldsList = {num_F-3, num_F-2, num_F-1};
num_F++;
Field[num_F] = MathEval;
Field[num_F].F = Sprintf("%g + %g * (1. - F%g)", lm, lM, num_F-1);
num_F++;
P[0] = newp; Point (P[0]) = {xmin, ymin, zmin};
P[1] = newp; Point (P[1]) = {xmax, ymin, zmin};
P[2] = newp; Point (P[2]) = {xmin, ymax, zmin};
P[3] = newp; Point (P[3]) = {xmax, ymax, zmin};
P[4] = newp; Point (P[4]) = {xmin, ymin, zmax};
P[5] = newp; Point (P[5]) = {xmax, ymin, zmax};
P[6] = newp; Point (P[6]) = {xmin, ymax, zmax};
P[7] = newp; Point (P[7]) = {xmax, ymax, zmax};
L[ 0] = newl; Line (L[ 0]) = {P[0], P[1]};
L[ 1] = newl; Line (L[ 1]) = {P[2], P[3]};
L[ 2] = newl; Line (L[ 2]) = {P[4], P[5]};
L[ 3] = newl; Line (L[ 3]) = {P[6], P[7]};
L[ 4] = newl; Line (L[ 4]) = {P[0], P[2]};
L[ 5] = newl; Line (L[ 5]) = {P[1], P[3]};
L[ 6] = newl; Line (L[ 6]) = {P[4], P[6]};
L[ 7] = newl; Line (L[ 7]) = {P[5], P[7]};
L[ 8] = newl; Line (L[ 8]) = {P[0], P[4]};
L[ 9] = newl; Line (L[ 9]) = {P[1], P[5]};
L[10] = newl; Line (L[10]) = {P[2], P[6]};
L[11] = newl; Line (L[11]) = {P[3], P[7]};
LL[0] = newll; Line Loop (LL[0]) = {L[ 0], L[ 5], -L[ 1], -L[ 4]};
LL[1] = newll; Line Loop (LL[1]) = {L[ 2], L[ 7], -L[ 3], -L[ 6]};
LL[2] = newll; Line Loop (LL[2]) = {L[ 8], L[ 2], -L[ 9], -L[ 0]};
LL[3] = newll; Line Loop (LL[3]) = {L[10], L[ 3], -L[11], -L[ 1]};
LL[4] = newll; Line Loop (LL[4]) = {L[ 4], L[10], -L[ 6], -L[ 8]};
LL[5] = newll; Line Loop (LL[5]) = {L[ 5], L[11], -L[07], -L[ 9]};
For k In {0:5}
S[k] = news; Plane Surface (S[k]) = {LL[k]};
EndFor
SL = newsl; Surface Loop (SL) = {S[0], -S[1], S[2], -S[3], S[4], -S[5]};
V = newv; Volume (V) = {SL};
Return
Function set_periodic_conditions
Periodic Surface S[1] {L[ 2], L[ 7], -L[ 3], -L[ 6]} = S[0] {L[ 0], L[ 5], -L[ 1], -L[ 4]};
Periodic Surface S[3] {L[10], L[ 3], -L[11], -L[ 1]} = S[2] {L[ 8], L[ 2], -L[ 9], -L[ 0]};
Periodic Surface S[5] {L[ 5], L[11], -L[07], -L[ 9]} = S[4] {L[ 4], L[10], -L[ 6], -L[ 8]};
Return
Function make_periodic_cube
Call make_cube;
Call set_periodic_conditions;
Return
Function make_periodic_heterogeneous_cube
Call make_heterogeneous_cube;
Call set_periodic_conditions;
Return