[Gmsh] how to make a periodic cubic mesh

[Aleksejs Fomins]

Dear All,

I have written a 2 layer cuboid using Transfinite Lines and Surfaces. The cuboid should have the top and bottom-most surface very coarse, and intermediate surface very fine.

Current problems:
1) Horizontal surfaces are fine. But for vertical surfaces GMSH complains like this

"Surface 37 can not be meshed using transfinite algo"

2) The mesh changes smoothly within the surface, but on the edge it stays constant, thus having too fine mesh towards the top and bottom

Please help
Aleksejs


On 14.12.2016 15:29, David Colignon wrote:
>
> On 14/12/16 14:22, Aleksejs Fomins wrote:
>
>> Dear Gmsh,
>>
>> I would like to create a cubic domain with a few complicated 3D shapes inside, and produce a corresponding tetrahedral mesh. The domain boundaries of the cube should be periodically conformal, e.g. each triangle on the surface with normal -x should be equivalent to a triangle on the surface with normal +x under translation along the x-axis.
>>
>> Using google I have found a solution proposed earlier on this list:
>>
>> http://onelab.info/pipermail/gmsh/2016/010305.html
>>
>> The .geo code suggested in the above link produces a conformal tetrahedral mesh of a cube using Transfinite surfaces and edges. It seems to work, but has a few problems:
>> 1) It completely ignores the value of the characteristic lengthscale variable "lc", always splitting each edge in exactly 4 parts
>> 2) The characteristic lengthscale changes depending on when index number of the Transfinite Line. Currently it is 5. If one would set it to 14, each edge would be split into exactly 25 segments. Why is that?
>
> Hi,
>
> see http://gmsh.info/doc/texinfo/gmsh.html#Structured-grids
>
> Transfinite Line { expression-list } | "*" = expression < Using Progression | Bump expression >;
>
>     Selects the lines in expression-list to be meshed with the 1D transfinite algorithm. The expression on the right hand side gives the number of nodes that will be created on the line (this overrides any other mesh element size prescription—see Specifying mesh element sizes).
>
>
> 5 nodes -> four segments
>
> Regards,
>
> Dave
>

-------------- next part --------------

gridsize = 2.0; // prescribed mesh element size

rCube = {5.0, 7.0};
hCube = {0.0, 4.0, 7.0};

// Points LVL 1
p_lv1[0] = newp;  Point(p_lv1[0]) = {-rCube[0], -rCube[1], hCube[0], gridsize};
p_lv1[1] = newp;  Point(p_lv1[1]) = {rCube[0], -rCube[1], hCube[0], gridsize};
p_lv1[2] = newp;  Point(p_lv1[2]) = {rCube[0], rCube[1], hCube[0], gridsize};
p_lv1[3] = newp;  Point(p_lv1[3]) = {-rCube[0], rCube[1], hCube[0], gridsize};

// Points LVL 1
p_lv2[0] = newp;  Point(p_lv2[0]) = {-rCube[0], -rCube[1], hCube[1], gridsize / 10.0};
p_lv2[1] = newp;  Point(p_lv2[1]) = {rCube[0], -rCube[1], hCube[1], gridsize / 10.0};
p_lv2[2] = newp;  Point(p_lv2[2]) = {rCube[0], rCube[1], hCube[1], gridsize / 10.0};
p_lv2[3] = newp;  Point(p_lv2[3]) = {-rCube[0], rCube[1], hCube[1], gridsize / 10.0};

p_lv3[0] = newp;  Point(p_lv3[0]) = {-rCube[0], -rCube[1], hCube[2], gridsize};
p_lv3[1] = newp;  Point(p_lv3[1]) = {rCube[0], -rCube[1], hCube[2], gridsize};
p_lv3[2] = newp;  Point(p_lv3[2]) = {rCube[0], rCube[1], hCube[2], gridsize};
p_lv3[3] = newp;  Point(p_lv3[3]) = {-rCube[0], rCube[1], hCube[2], gridsize};

// Lines - Horizontal - LVL 1
linh_lv1[0] = newl;  Line(linh_lv1[0]) = {p_lv1[0], p_lv1[1]}; 
linh_lv1[1] = newl;  Line(linh_lv1[1]) = {p_lv1[1], p_lv1[2]};
linh_lv1[2] = newl;  Line(linh_lv1[2]) = {p_lv1[2], p_lv1[3]};
linh_lv1[3] = newl;  Line(linh_lv1[3]) = {p_lv1[3], p_lv1[0]};

// Lines - Horizontal - LVL 2
linh_lv2[0] = newl;  Line(linh_lv2[0]) = {p_lv2[0], p_lv2[1]}; 
linh_lv2[1] = newl;  Line(linh_lv2[1]) = {p_lv2[1], p_lv2[2]};
linh_lv2[2] = newl;  Line(linh_lv2[2]) = {p_lv2[2], p_lv2[3]};
linh_lv2[3] = newl;  Line(linh_lv2[3]) = {p_lv2[3], p_lv2[0]};

// Lines - Horizontal - LVL 3
linh_lv3[0] = newl;  Line(linh_lv3[0]) = {p_lv3[0], p_lv3[1]}; 
linh_lv3[1] = newl;  Line(linh_lv3[1]) = {p_lv3[1], p_lv3[2]};
linh_lv3[2] = newl;  Line(linh_lv3[2]) = {p_lv3[2], p_lv3[3]};
linh_lv3[3] = newl;  Line(linh_lv3[3]) = {p_lv3[3], p_lv3[0]};

// Lines - Vertical - LVL 1
linv_lv1[0] = newl;  Line(linv_lv1[0]) = {p_lv1[0], p_lv2[0]}; 
linv_lv1[1] = newl;  Line(linv_lv1[1]) = {p_lv1[1], p_lv2[1]};
linv_lv1[2] = newl;  Line(linv_lv1[2]) = {p_lv1[2], p_lv2[2]};
linv_lv1[3] = newl;  Line(linv_lv1[3]) = {p_lv1[3], p_lv2[3]};

// Lines - Vertical - LVL 2
linv_lv2[0] = newl;  Line(linv_lv2[0]) = {p_lv2[0], p_lv3[0]}; 
linv_lv2[1] = newl;  Line(linv_lv2[1]) = {p_lv2[1], p_lv3[1]};
linv_lv2[2] = newl;  Line(linv_lv2[2]) = {p_lv2[2], p_lv3[2]};
linv_lv2[3] = newl;  Line(linv_lv2[3]) = {p_lv2[3], p_lv3[3]};

// Line Loops
llh_lv1[0] = newreg;  Line Loop(llh_lv1[0]) = {linh_lv1[]};
llh_lv2[0] = newreg;  Line Loop(llh_lv2[0]) = {linh_lv2[]};
llh_lv3[0] = newreg;  Line Loop(llh_lv3[0]) = {linh_lv3[]};

llv_lv1[0] = newreg;  Line Loop(llv_lv1[0]) = {linh_lv1[0], linv_lv1[1], -linh_lv2[0], -linv_lv1[0]};
llv_lv1[1] = newreg;  Line Loop(llv_lv1[1]) = {linh_lv1[1], linv_lv1[2], -linh_lv2[1], -linv_lv1[1]};
llv_lv1[2] = newreg;  Line Loop(llv_lv1[2]) = {linh_lv1[2], linv_lv1[3], -linh_lv2[2], -linv_lv1[2]};
llv_lv1[3] = newreg;  Line Loop(llv_lv1[3]) = {linh_lv1[3], linv_lv1[0], -linh_lv2[3], -linv_lv1[3]};

llv_lv2[0] = newreg;  Line Loop(llv_lv2[0]) = {linh_lv2[0], linv_lv2[1], -linh_lv3[0], -linv_lv2[0]};
llv_lv2[1] = newreg;  Line Loop(llv_lv2[1]) = {linh_lv2[1], linv_lv2[2], -linh_lv3[1], -linv_lv2[1]};
llv_lv2[2] = newreg;  Line Loop(llv_lv2[2]) = {linh_lv2[2], linv_lv2[3], -linh_lv3[2], -linv_lv2[2]};
llv_lv2[3] = newreg;  Line Loop(llv_lv2[3]) = {linh_lv2[3], linv_lv2[0], -linh_lv3[3], -linv_lv2[3]};

// Surfaces
sh_lv1[0] = news;  Plane Surface(sh_lv1[0]) = {llh_lv1[0]};
sh_lv2[0] = news;  Plane Surface(sh_lv2[0]) = {llh_lv2[0]};
sh_lv3[0] = news;  Plane Surface(sh_lv3[0]) = {llh_lv3[0]};
sv_lv1[0] = news;  Plane Surface(sv_lv1[0]) = {llv_lv1[0]};
sv_lv1[1] = news;  Plane Surface(sv_lv1[1]) = {llv_lv1[1]};
sv_lv1[2] = news;  Plane Surface(sv_lv1[2]) = {llv_lv1[2]};
sv_lv1[3] = news;  Plane Surface(sv_lv1[3]) = {llv_lv1[3]};
sv_lv2[0] = news;  Plane Surface(sv_lv2[0]) = {llv_lv2[0]};
sv_lv2[1] = news;  Plane Surface(sv_lv2[1]) = {llv_lv2[1]};
sv_lv2[2] = news;  Plane Surface(sv_lv2[2]) = {llv_lv2[2]};
sv_lv2[3] = news;  Plane Surface(sv_lv2[3]) = {llv_lv2[3]};

// Surface Loops and Volumes
sl[0] = newreg;  Surface Loop(sl[0]) = {sh_lv1[], sh_lv2[], sv_lv1[]};
sl[1] = newreg;  Surface Loop(sl[1]) = {sh_lv2[], sh_lv3[], sv_lv2[]};
vol[0] = newv;   Volume(vol[0]) = {sl[0]};
vol[1] = newv;   Volume(vol[1]) = {sl[1]};

// Transfinite lines
Transfinite Line{linh_lv1[]} = 5;
Transfinite Line{linh_lv2[]} = 20;
Transfinite Line{linh_lv3[]} = 5;
Transfinite Line{linv_lv1[]} = 5;
Transfinite Line{linv_lv2[]} = 7;

// Transfinite Surfaces
// need to specify orientation (through vertex list) by hand to have correctly
// matching volume/surface edges
Transfinite Surface {sh_lv1[0]} = {p_lv1[]};
Transfinite Surface {sh_lv3[0]} = {p_lv3[]};
Transfinite Surface {sv_lv1[0]} = {p_lv1[0], p_lv1[1], p_lv2[1], p_lv2[0]};
Transfinite Surface {sv_lv1[1]} = {p_lv1[1], p_lv1[2], p_lv2[2], p_lv2[1]};
Transfinite Surface {sv_lv1[2]} = {p_lv1[2], p_lv1[3], p_lv2[3], p_lv2[2]};
Transfinite Surface {sv_lv1[3]} = {p_lv1[3], p_lv1[0], p_lv2[0], p_lv2[3]};
Transfinite Surface {sv_lv2[0]} = {p_lv2[0], p_lv2[1], p_lv3[1], p_lv3[0]};
Transfinite Surface {sv_lv2[1]} = {p_lv2[1], p_lv2[2], p_lv3[2], p_lv3[1]};
Transfinite Surface {sv_lv2[2]} = {p_lv2[2], p_lv2[3], p_lv3[3], p_lv3[2]};
Transfinite Surface {sv_lv2[3]} = {p_lv2[3], p_lv2[0], p_lv3[0], p_lv3[3]};

// Transfinite Volumes
Transfinite Volume{vol[0]} = {p_lv1[], p_lv2[]};
Transfinite Volume{vol[1]} = {p_lv2[], p_lv3[]};