Assigning labels to surface mesh


I am trying to create a mesh on the parallelepiped (Xmin,Xmax)-by-(Ymin,Ymax)-by-(Zmin,Zmax). The mesh must be uniform in the x- and z-directions, but non-uniform in the y-direction, so I cannot use cube. My approach is to create meshes on the six faces of the parallelepiped using triangulate() followed by movemesh23; then glue these meshS objects together and feed the resulting mesh to tetg. It works and the resulting mesh looks good, but somehow the labels of the faces (which I will need for imposing boundary conditions later) get lost. First, when I glue the six faces together, only the labels of the first four get remembered. Then, after the call to tetg, all the faces of the resulting parallelepiped have the same label zero. Here are the pertinent parts of the code. Any help will be greatly appreciated.

int[int] rinit1 = [1, 31], rinit2=[1, 32], rinit3=[1, 33], r1=[33,10000],r2=[33,20000],r3=[32,30000],r4=[32,40000],r5=[31,50000],r6=[31,60000];

mesh Thyz=triangulate(“yzf”);
mesh Thyzlab = change(Thyz,label=rinit1);

mesh Thxz=triangulate(“xzf”);
mesh Thxzlab = change(Thxz,label=rinit2);

mesh Thxy=triangulate(“xyf”);
mesh Thxylab = change(Thxy,label=rinit3);

meshS Thxy1 = movemesh23(Thxylab, transfo=[XX1, YY1, ZZ1max], label=r1, orientation=1);
meshS Thxy2 = movemesh23(Thxylab, transfo=[XX1, YY1, ZZ1min], label=r2, orientation=-1);
meshS Thxz1 = movemesh23(Thxzlab, transfo=[XX2, YY2max, ZZ2], label=r3, orientation=-1);
meshS Thxz2 = movemesh23(Thxzlab, transfo=[XX2, YY2min, ZZ2], label=r4, orientation=1);
meshS Thyz1 = movemesh23(Thyzlab, transfo=[XX3max, YY3, ZZ3], label=r5, orientation=1);
meshS Thyz2 = movemesh23(Thyzlab, transfo=[XX3min, YY3, ZZ3], label=r6, orientation=-1);

meshS ThSurf = Thxy1 + Thxy2 + Thxz1 + Thxz2 + Thyz1 + Thyz2;

real voltet = ((XmaxVal-XminVal)/N)((YmaxVal-YminVal)/dofY)((ZmaxVal-ZminVal)/N)/6.0;
real[int] domaintetg =[Xmid, Ymid, Zmid, 1, voltet];
mesh3 Th3 = tetg(ThSurf, switch=“pqaAAYYQ”, nbofregions=1, regionlist=domaintetg);

Here is how I check the labels of the resulting meshes, e.g., he surface mesh ThSurf:

int[int] lab = labels(ThSurf);
cout<< "lab.n = “<< lab.n<<endl;
for(int i=0;i<lab.n;i++)
cout<< “lab(”<<i<<”) = "<<lab(i)<<endl;

Here is the result:
lab.n = 4
lab(0) = 10000
lab(1) = 20000
lab(2) = 30000
lab(3) = 40000

If I change the order in which I glue meshes (keeping correct orientation), I still see only the labels of the first four meshes being glued. And, as I said earlier, all the six faces of the resulting 3D mesh Th3 have label zero.


  1. you want to change region number in 2d or in meshS case not label/

mesh Thyzlab = change(Thyz,fregion=31);/ use region => function to set region

meshS Thxy1 = movemesh23(Thxylab, transfo=[XX1, YY1, ZZ1max], region=r1, orientation=1);

int[int] lab = regions(ThSurf);

This worked, thank you very much!