//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mesh_construction/mesh_construction.edp
//
//  Discussion:
//
//    Construct a mesh for a contraction flow problem.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 May 2020
//
//  Reference:
//
//    John Chrispell, Jason Howell,
//    Finite Element Approximation of Partial Differential Equations
//    Using FreeFem++, or, How I Learned to Stop Worrying and Love
//    Numerical Analysis.
//
//    Frederic Hecht,
//    New development in FreeFem++,
//    Journal of Numerical Mathematics,
//    Volume 20, Number 3-4, 2012, pages 251-265.
//
cout << "\n";
cout << "mesh_construction\n";
cout << "  FreeFem++ version\n";
cout << "  Demonstrate how to define a region with line segments,\n";
cout << "  and then mesh it.\n";
//
//  Use N to control the relative fineness of the mesh.
//  N = 1 will be the coarsest possible mesh.
//
int n = 3;
//
//  Define line segments that form the boundary.
//  We break up each line so that we can force parts to have a finer mesh.
//
// Bottom edge.
//
border l00 ( t = 0, 1 ) { x = 0.0 + 2.0 * t; y = 0.0;               label = 1; }; 
border l01 ( t = 0, 1 ) { x = 2.0 + 1.5 * t; y = 0.0;               label = 1; };
border l02 ( t = 0, 1 ) { x = 3.5 + 1.0 * t; y = 0.0;               label = 1; };
border l03 ( t = 0, 1 ) { x = 4.5 + 3.5 * t; y = 0.0;               label = 1; };
//
// Outflow edge
//
border l04 ( t = 0, 1 ) { x = 8.0;           y = 0.0   + 0.125 * t; label = 2; };
border l05 ( t = 0, 1 ) { x = 8.0;           y = 0.125 + 0.125 * t; label = 2; };
//
// Top of contraction channel
//
border l06 ( t = 0, 1 ) { x = 8.0 - 3.5 * t; y = 0.250;             label = 3; }; 
border l07 ( t = 0, 1 ) { x = 4.5 - 0.5 * t; y = 0.250;             label = 3; };
//
// Contraction wall
//
border l08 ( t = 0, 1 ) { x = 4.0;           y = 0.250 + 0.125 * t; label = 4; }; 
border l09 ( t = 0, 1 ) { x = 4.0;           y = 0.375 + 0.250 * t; label = 4; };
border l10 ( t = 0, 1 ) { x = 4.0;           y = 0.625 + 0.375 * t; label = 4; };
//
// Top of inflow channel
//
border l11 ( t = 0, 1 ) { x = 4.0 - 0.5 * t; y = 1.0;               label = 5; }; 
border l12 ( t = 0, 1 ) { x = 3.5 - 1.5 * t; y = 1.0;               label = 5; };
border l13 ( t = 0, 1 ) { x = 2.0 - 2.0 * t; y = 1.0;               label = 5; };
//
// Inflow walls.
//
border l14 ( t = 0, 1 ) { x = 0.0;           y = 1.000 - 0.375 * t; label = 6; };
border l15 ( t = 0, 1 ) { x = 0.0;           y = 0.625 - 0.250 * t; label = 6; };
border l16 ( t = 0, 1 ) { x = 0.0;           y = 0.375 - 0.250 * t; label = 6; };
border l17 ( t = 0, 1 ) { x = 0.0;           y = 0.125 - 0.125 * t; label = 6; };
//
//  Build the mesh.
//  We particularly want a fine mesh near the transition region at x = 4.0.
//
mesh Th = buildmesh 
(
    l00 ( 4 * n ) // Bottom edge
  + l01 ( 4 * n )
  + l02 ( 8 * n )
  + l03 ( 10 * n )
  + l04 ( 4 * n ) // Outflow edge
  + l05 ( 4 * n )
  + l06 ( 10 * n ) // Top of contraction channel
  + l07 ( 4 * n )
  + l08 ( 4 * n ) // Contraction wall
  + l09 ( 4 * n ) 
  + l10 ( 4 * n )
  + l11 ( 4 * n ) // Top of inflow channel
  + l12 ( 4 * n )
  + l13 ( 4 * n )
  + l14 ( 4 * n ) // Inflow walls.
  + l15 ( 4 * n )
  + l16 ( 8 * n )
  + l17 ( 2 * n )
);
//
//  Display the mesh.
//
plot ( Th, wait = true, ps = "mesh_construction.ps" );
//
//  Terminate.
//
cout << "\n";
cout << "mesh_construction\n";
cout << "  Normal end of execution.\n";

exit ( 0 );