Periodic BCs in 3D


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. Next, I need to introduce FE spaces with periodic boundary conditions in the x and z directions - and I get the error message. I have seen other messages regarding this type of error in the General Discussion, but none have been answered.

Here are the pertinent parts of the code, along with the error message. Any help will be greatly appreciated.

int[int] 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,fregion=31);

mesh Thxz=triangulate(“xzf”);
mesh Thxzlab = change(Thxz,fregion=32);

mesh Thxy=triangulate(“xyf”);
mesh Thxylab = change(Thxy,fregion=33);

meshS Thxy1 = movemesh23(Thxylab, transfo=[XX1, YY1, ZZ1max], region=r1, orientation=1);
meshS Thxy2 = movemesh23(Thxylab, transfo=[XX1, YY1, ZZ1min], region=r2, orientation=-1);
meshS Thxz1 = movemesh23(Thxzlab, transfo=[XX2, YY2max, ZZ2], region=r3, orientation=-1);
meshS Thxz2 = movemesh23(Thxzlab, transfo=[XX2, YY2min, ZZ2], region=r4, orientation=1);
meshS Thyz1 = movemesh23(Thyzlab, transfo=[XX3max, YY3, ZZ3], region=r5, orientation=1);
meshS Thyz2 = movemesh23(Thyzlab, transfo=[XX3min, YY3, ZZ3], region=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);

fespace Vh2(Th3,P2, periodic=[[60000,y,z],[50000,y,z],[20000,x,y],[10000,x,y]]);

fespace Vh(Th3,P1, periodic=[[60000,y,z],[50000,y,z],[20000,x,y],[10000,x,y]]);

problem ADC1 ([u1,u2,u3,mixDzeta1,mixDzeta2,mixDzeta3,p,lambda1],[v1,v2,v3,mixv1,mixv2,mixv3,q,qlambda1], solver=UMFPACK) =
(1.0/dt)( u1v1 + u2v2 + u3v3 ));

– Build Nodes/DF on mesh : n.v. 609, n. elmt. 2625, n b. elmt. 768
nb of Nodes 4226 nb of DoF 26574 DFon=8600
Problem build of GFESpace (3d) (may be : due to periodic Boundary condition missing ) FH
The number of DF must be 22902 and it is 26574
current line = 1163
Assertion fail : (snbdf == NbOfDF)
line :560, in file …/femlib/FESpacen.cpp
Assertion fail : (snbdf == NbOfDF)
line :560, in file …/femlib/FESpacen.cpp
err code 6 , mpirank 0

I was able to figure this out on my own. Apparently, the FreeFEM++ bug (from at least 2012) persists, which doesn’t allow to “mix” periodic FE spaces. Thus, I had to introduce the vector valued FE space via

fespace Vh2(Th3,[P2,P2,P2,P1,P2,P2,P2,P1], periodic=[[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);

This resolved the issue.