Using mshmet and tetgreconstruction with periodic domains

General Discussion
1

Hi all, I am trying to implement an adaptive meshing scheme using mshmet and tetgreconstruction for my 3D, periodic system.

I am getting an error stating that

Exec error : Periodic 3d:  the both number of face is not the same

I am running it with isotropic remeshing, so I don’t think it should be causing any changes. Do you know of any workaround to enforce continuity on the periodic faces?

I know that there is the YY switch in tetgreconstruction that keep the boundary values, but that kind of defeats my point if there are errors at the boundaries.

I have included my code below, but the key part is probably somewhere in
Th3 = tetgreconstruction(Th3, switch="raAQ", sizeofvolume=h*h*h/6.)


load "msh3"
load "medit"
load "tetgen"
load "mshmet"
load "MUMPS_seq"
searchMethod=1; // More safe seach algo

/// Rules about macros
/// MACROS MUST END WITH THE TWO FORWARDSLASHES "//"
  macro Grad(u) [dx(u),dy(u)]	    // End of Macro
  macro Lap(u,v) (Grad(u)'*Grad(v)) //  End of Macro
  macro Cdot(fL,fR) (fL'*fR)  // End of Macro
  macro RotVec(u) [cos(u),sin(u)] // End of Macro






real zmin = 0;
real zmax = 2*pi;
real errm = 1e-2; //level of error

// Mesh 2D
real ax = 5.;
real ay = 5.;

int nnouter = 20 , nz = 10;


border rightborder(t=-0.5, 0.5){x=ax/2.; y=ay*t; label=2;}
border leftborder(t=-0.5, 0.5){x=-ax/2.; y=ay*t; label=4;}
border bordero1(t=0.5, -0.5){x=ax*t; y=-ay/2; label=1;}//Bottom
border bordero3(t=0.5, -0.5){x=ax*t; y=ay/2; label=3;}//Top
mesh Th = buildmesh(bordero1(-nnouter) + bordero3(nnouter) + rightborder(nnouter) + leftborder(-nnouter) );





{ // Cleaning the memory correctly
    int[int] old2new(0:Th.nv-1);
    fespace Vh2(Th, P1);
    Vh2 sorder = x + y;
    sort(sorder[], old2new);
    int[int] new2old = old2new^-1; // Inverse permutation
    Th = change(Th, renumv=new2old);
    sorder[] = 0:Th.nv-1;
}
{
    fespace Vh2(Th, P1);
    Vh2 nu;
    nu[] = 0:Th.nv-1;
    //plot(nu, cmm="nu=", wait=true);
}

// Mesh 3D
int[int] rup = [0, 5], rlow = [0, 6], rmid = [1, 1, 2, 2, 3, 3, 4, 4], rtet = [0, 41]; 
mesh3 Th3 = buildlayers(Th, nz, zbound=[zmin, zmax],
    reftet=rtet, reffacemid=rmid, reffaceup=rup, reffacelow=rlow);
for(int i = 1; i <= 6; ++i)
    cout << " int " << i << " : " << int2d(Th3,i)(1.) << " " << int2d(Th3,i)(1./area) << endl;



fespace Vh(Th3, P1, periodic=[ [5, x, y], [6, x, y]]);
int n = Vh.ndof;
int n1 = n+1;
Vh uh, vh, dxu, dyu , usol, h; // probabiility and test function.




// Problem
problem va(uh,vh,solver=sparsesolver) = // definion of the problem
  int3d(Th3)(dx(vh)*dx(uh) + dy(vh)*dy(uh) + dz(uh)*dz(vh))
  + on(1,uh=1) // bottom hand wall
  + on(3,uh=0.9) // top hand wall
;


// Loop
for (int ii = 0; ii < 5; ii++){
    // Solve
    va;
    dxu = dx(uh);
    dyu = dy(uh);

    cout << "u min, max = " << uh[].min << " "<< uh[].max << endl;

    h=0.; //for resizing h[] because the mesh change
    h[] = mshmet(Th3, uh, normalization=1, aniso=0, nbregul=1, hmin=1e-3, hmax=0.3, err=errm);
    cout << "h min, max = " << h[].min << " "<< h[].max << " " << h[].n << " " << Th3.nv << endl;
    plot(uh, wait=true);

    errm *= 0.8; //change the level of error
    cout << "Th3 " << Th3.nv < " " << Th3.nt << endl;
    Th3 = tetgreconstruction(Th3, switch="raAQ", sizeofvolume=h*h*h/6.); //rebuild mesh
    medit("U-adap-iso-"+ii, Th3, uh, wait=true);
}