//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mshmet/mshmet.edp
//
//  Discussion:
//
//    The domain is a 3D version of the L shaped region, a cube with
//    a corner removed.
//
//    Morice's mshmet plugin computes an isotropic adapted mesh.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 May 2020
//
//  Reference:
//
//    Frederic Hecht,
//    New development in FreeFem++,
//    Journal of Numerical Mathematics,
//    Volume 20, Number 3-4, 2012, pages 251-265.
//
cout << "\n";
cout << "mshmet:\n";
cout << "  FreeFem++ version\n";

load "msh3" 
load "tetgen" 
load "mshmet" 
load "medit"

int n = 6;
int[int] l1 = [1,1,1,1],l01=[0,1],l11=[1,1];
//
//  Label numbering for boundary condition.
//
mesh3 Th3 = buildlayers ( square ( n, n, region = 0, label = l1 ),
  n, zbound = [0,1], labelmid = l11, labelup = l01,
  labeldown = l01 );
//
//  Remove the [0.5,1.0]^3 cube.
//
Th3 = trunc ( Th3, ( x < 0.5 ) | ( y < 0.5 ) | ( z < 0.5 ),
  label = 1 ); 

fespace Vh ( Th3, P1 );
Vh u, v, usol;

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

solve Poisson ( u, v, solver = CG ) =
    int3d ( Th3 ) ( 1.0e-6*u*v + Grad(u)' * Grad(v) )
  - int3d ( Th3 ) ( ( z - y / 2 - x / 2 ) * v );

real errm = 1.0e-02;

for ( int ii = 0; ii < 5; ii++ )
{
  Vh h;
  h[] = mshmet ( Th3, u, normalization = 1, aniso = 0, nbregul = 1,
    hmin = 1.0E-03, hmax = 0.3, err = errm );

  errm = errm * 0.6;

  Th3 = tetgreconstruction ( Th3, switch = "raAQ",
    sizeofvolume = h*h*h / 6.0 );
  Poisson;
//
//  A 3D plot can't be saved to a PostScript file.
//
  plot ( u, wait = true, nbiso = 15 );
}
//
//  Invoke medit to display the solution.
//
medit ( "U3", Th3, u );
//
//  Terminate.
//
cout << "\n";
cout << "mshmet\n";
cout << "  Normal end of execution.\n";

exit ( 0 );