load "gmsh"

int nn=20;
mesh Th2=square(nn,nn,region=0);
fespace Vh2(Th2,P2);
Vh2 ux,uz,p2;
int[int] rup=[0,2],  rdown=[0,1], rmid=[1,1,2,1,3,1,4,1];
real zmin=0,zmax=1;
mesh3 Th=buildlayers(Th2,nn,
  zbound=[zmin,zmax],
  labelmid=rmid,
  reffaceup = rup,
  reffacelow = rdown);

Th = buildBdMesh(Th); // build the surface mesh
meshS ThS= Th.Gamma;

fespace Vh(ThS,P2);
Vh u,v;

func ue =   2*x*x + 3*y*y + 4*z*z + 5*x*y+6*x*z+1;
func uex=   4*x+  5*y+6*z;
func uey=   6*y + 5*x;
func uez=   8*z +6*x;
func f= -18. ;

macro Grad3(u) [dx(u),dy(u),dz(u)]  // EOM
// 
problem Lap3dS(u,v,solver=CG)  =
  int2d(ThS)(Grad3(v)' *Grad3(u)) //') for emacs
  + int1d(ThS,2)(u*v)
  - int2d(ThS)(f*v)
  - int1d(ThS,2) ( ue*v + (uex*N.x +uey*N.y +uez*N.z)*v )
  + on(1,u=ue);
Lap3dS;

plot(u,wait=1);


savegmsh(ThS,"ThS");

 

