Weird behaviour of 3D Periodic spaces

Hello everyone,

While experimenting with periodic elements on a cubic mesh with an inclusion, I tried to write down a basic loop to display the values of my solution at the interface. The code gives back an array length error telling me that the periodic function has less DOF than the mesh. The following is a minimal script I wrote to reproduce the error.

Sincerly yours,
Salah

load "medit"
load "mmg"
load "msh3"

int nn = 10;
int[int] rt=[1,2,3,4,5,6]; //required triangles, kept unchanged by mmg;

mesh3 Th=cube(nn,nn,nn,[x,y,z]);

func sphere = (x-0.5)^2+(y-0.5)^2+(z-0.5)^2-0.2^2; // level_set initialisation of a sphere

fespace Vh(Th,P1);
Vh phi=sphere;

Th=mmg3d(Th,metric=phi[],iso=1,ls=0,hausd=0.008,hgrad=1,hmin=0.01,hmax=0.1,mem=1000,requiredTriangle=rt); // implicit meshing by mmg

medit("mesh",Th);// display the resulting mesh;

fespace Vhp(Th,P1,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]); //periodic element fespace

Vh F;
Vhp Fp;

cout << " The DOF of the mesh is " << Th.nv << ", The DOF of Vh is " << F[].n << ", finally, the DOF of Vhp (periodic elements space) is " << Fp[].n << endl;

It gives back the following message :

 The DOF of the mesh is 2551, The DOF od Vh is   -- FESpace: Nb of Nodes 2551 Nb of DoF 2551
2551, finally, the DOF of Vhp (periodic elements space) is 2220

Can anyone please tell me what I’m missing about the 331 difference in DOF? Thank you in advance for your help ^^

Best regards,
Salah

Yes, the problem is to day that no periodic mesh generation in mmg3d.

Sorry for that, see we the mmg team.

I don’t know what Frédéric is implying, using requiredTriangle is the proper way to do periodic mesh adaptation with mmg3d. There is a difference because periodic degrees of freedom are not counted multiple times (on each periodic face).

Thank you for your reply,

Does this only occur in 3D? because I don’t recall having this problem with 2D periodic elements, Is there a map that assigns to every mesh nod its adress on a periodic element array? since, for example, u[Th(i)] will not work for all nods.