
// Parameters
int nn = 30; // Mesh quality

// Mesh
int[int] labs = [1, 2, 2, 1, 1, 2]; // Label numbering
mesh3 Th = cube(nn, nn, nn, label=labs);
// extract the surface (boundary) of the cube
int[int] llabs = [1, 2]; //  Th = buildSurface(Th); 
meshS ThS = extract(Th,label=llabs);


// Remove the ]0.5,1[^3 domain of the surface cube
ThS = trunc(ThS, (x < 0.5) | (y < 0.5) | (z < 0.5), label=1);
plot(ThS);
// Fespace
fespace Vh(ThS, P1);
Vh u, v;

// Macro
macro Grad(u) [dx(u), dy(u), dz(u)] //

// Define the weak form and solve
solve Poisson(u, v, solver=CG)
    = int2d(ThS)(
          Grad(u)' * Grad(v)
    )
    -int2d(ThS)(
          1 * v
    )
    + on(1, u=0)
    ;

// Plot
plot(u, nbiso=15);
