Need a border cube but it looks like flags in msh3.cpp commented out

I could not ifind a way to make cube “holes” in a cubes mesh. I thought bcube may
do that but looking here, plugin/seq/msh3.cpp, it looks like there may be
flag bits for that but none seem to be coded in. Is there some simple way to do this?

Before getting to the Dirac Eqn, I thought I would play with this, lol,

This code sort of runs except it seems to add the BC to each cube and solve inside
each one. And the adaptive thing does not work but I can probably figure that out.

load “msh3”
load “medit”
real h=.25;
real s1=.25;
real s2=.25;
real cen=1;
real first=1.5;
real d=.6;
mesh3 Th3 = cube(20,20,20,[2x,2y,2*z]);
mesh3 Th4 = bcube(10,10,10, [s1*x+first,s1*y+cen,s1*z+h]);
mesh3 Th5 = bcube(10,10,10, [s2*x+first-d,s2*y+cen,s2*z+h]);
mesh3 Thx=Th3+Th4+Th5;
mesh3 Tsb=buildBdMesh(Th3);
//meshS surface=Tsb.Gamma;
int[int] labs = [1];
//meshS Thl = extract(Th3.Gamma, label=labs);
meshS Thl = extract(Th3, label=labs);


real tol=.5;
real cut=.01;

for(int i=0; i<1; ++i)

fespace Xh(Thx,P13d);
Xh psi, wx;
macro grad(u) [dx(u),dy(u),dz(u)]// def of grad operator
// Solve
solve potential(psi, wx)
= int3d(Thx)(
grad(psi)'*grad(wx)) // scalar product
plot(psi, value=true, wait=1);

//Thx = adaptmesh(Thx, [dx(psi), dy(psi),dz(psi)], splitpbedge=1, nbvx=1e5, abserror=1, cutoff=cut, err=tol, inquire=0, ratio=1.5, hmin=1./100);

} // i

in include file MeshSurface.idp

like in file example/3d/Poisson-cube-ballon.edp or pluton

Thanks, I will take another look through the examples but that is difficult without
some uncommon key words :slight_smile: I did get code to run using, SurfaceHex that I found
in the manual but the orientation flag does not seem to do anything and I have
yet to look at the cpp source code.
The code I have is below but let me look at those examples before doing anything
else. Thanks.

real h=-.5;
real s1=.5;
real s2=.5;
real cen=1;
real first=1.;
real d=.6;
int[int,int] label1=[[101,102],[103,104],[105,106]];
int[int] NX1=[20,20,20];
int[int,int] label2=[[201,202],[203,204],[205,206]];
int[int] NX2=[20,20,20];
real [int,int] BX1=[[first,s1+first],[cen,cen+s1],[h,h+s1]];
real [int,int] BX2=[[first-d,s2+first-d],[cen,cen+s2],[h,h+s2]];
int[int] NX3=[20,20,20];
real [int,int] BX3=[[0,2],[0,2],[-1,1 ]];
int[int,int] label3=[[1,2],[3,4],[5,6]];
meshS Th5=SurfaceHex(NX1,BX1,label1,1);
meshS Th3=SurfaceHex(NX2,BX2,label2,1);
meshS Th4=SurfaceHex(NX3,BX3,label3,-1);
meshS Thxs=Th3+Th4+Th5;
mesh3 Thx = tetg(Thxs, switch=“pqaAAYYQ”); // , nbofregions=1);