// test of mmg3d for move objets in a mesh ...
load "msh3" 
load "tetgen" 
load "medit" 
load "mmg3d-v4.0"                                                               
include "MeshSurface.idp"

// build mesh of a box (311)  wit 2 holes  (300,310)

real hs = 0.8; 
int[int]  N=[4/hs,8/hs,11.5/hs];
real [int,int]  B=[[-2,2],[-2,6],[-10,1.5]];
int [int,int]  L=[[311,311],[311,311],[311,311]];
mesh3 ThH = SurfaceHex(N,B,L,1);
mesh3 ThSg =Sphere(1,hs,300,-1); // "gluing" surface meshs to tolat boundary meshes
mesh3 ThSd =Sphere(1,hs,310,-1); 

ThSd=movemesh3(ThSd,transfo=[x,4+y,z]);

mesh3 ThHS=ThH+ThSg+ThSd;
medit("ThHS", ThHS);

real voltet=(hs^3)/6.;
cout << " voltet = " << voltet << endl;
real[int] domaine = [0,0,-4,1,voltet];
real [int] holes=[0,0,0,0,4,0];
mesh3 Th = tetg(ThHS,switch="pqaAAYYQ",regionlist=domaine,holelist=holes);    
medit("Box-With-two-Ball",Th);
// End build mesh 

int[int] opt=[9,0,64,0,0,3];   // options  of mmg3d see freeem++ doc 

real[int] vit=[0,0,-0.3];
func zero = 0.;
func dep  = vit[2];

fespace Vh(Th,P1); 
macro Grad(u) [dx(u),dy(u),dz(u)] //

Vh uh,vh; //  to compute the displacemnt field 
problem Lap(uh,vh,solver=CG) = int3d(Th)(Grad(uh)'*Grad(vh))  //') for emacs
				  + on(310,300,uh=dep) +on(311,uh=0.); 

for(int it=0; it<29; it++){ 
  cout<<"  ITERATION       "<<it<<endl;
  Lap;
medit("Box-With-two-Ball",Th,uh);
  Th=mmg3d(Th,opt="-O 9",displacement=[zero,zero,uh]); 
 }