/***************************************************************************************************
*                            triangulated_rectangles_concatenation.geo                             *
****************************************************************************************************
* AUTHOR: Alejandro Pinio (alepinio@gmail.com)                                                     *
*                                                                                                  * 
* SUMMARY: geometry of a two dimensional domain built from the concatenation of rectangles         *
* divided into 4 triangles. Meshing does not introduce additional nodes nor elements. Rectangles   *
* can have different size and physical name.                                                       *
*                                                                                                  *
* REFERENCES: t5.geo file from the tutorial directory of the Gmsh 2.10.1 distribution.             *
***************************************************************************************************/

Macro GenerateTriangulatedRectangle
  // ~> points
  // points were labelled according to their mappings into natural coordinates and the Figure 3.2.1
  // from [Hughes2000, p. 113]
  _cl = 1.0;
  _p1 = newp; Point (_p1) = {_xMin, _yMin, 0.0, _cl};
  _p2 = newp; Point (_p2) = {_xMax, _yMin, 0.0, _cl};
  _p3 = newp; Point (_p3) = {_xMax, _yMax, 0.0, _cl};
  _p4 = newp; Point (_p4) = {_xMin, _yMax, 0.0, _cl};
  _p5 = newp; Point (_p5) = {_xMin+0.5*(_xMax-_xMin), _yMin+0.5*(_yMax-_yMin), 0.0, _cl};

  // ~> lines
  _l1 = newl; Line (_l1)  = {_p1, _p2}; // line parallel to x direction
  _l2 = newl; Line (_l2)  = {_p2, _p3}; 
  _l3 = newl; Line (_l3)  = {_p3, _p4}; // line parallel to x direction
  _l4 = newl; Line (_l4)  = {_p4, _p1};
  _l5 = newl; Line (_l5)  = {_p1, _p5};
  _l6 = newl; Line (_l6)  = {_p2, _p5};
  _l7 = newl; Line (_l7)  = {_p3, _p5};
  _l8 = newl; Line (_l8)  = {_p4, _p5};

  // ~> line loop
  // outward normals
  _ll1 = news; Line Loop (_ll1) = {_l1, _l6, -_l5};
  _ll2 = news; Line Loop (_ll2) = {_l2, _l7, -_l6};
  _ll3 = news; Line Loop (_ll3) = {_l3, _l8, -_l7};
  _ll4 = news; Line Loop (_ll4) = {_l4, _l5, -_l8};

  // ~> surface
  _s1 = news; Plane Surface (_s1) = {_ll1};
  _s2 = news; Plane Surface (_s2) = {_ll2};
  _s3 = news; Plane Surface (_s3) = {_ll3};
  _s4 = news; Plane Surface (_s4) = {_ll4};

  // ~> append the ids of the surfaces to the list containing the ids of the surfaces of the whole physical surface
  Physical Surface(_physicalSurfaceName) += {_s1};
  Physical Surface(_physicalSurfaceName) += {_s2};
  Physical Surface(_physicalSurfaceName) += {_s3};
  Physical Surface(_physicalSurfaceName) += {_s4};

  // ~> rectangle meshing settings
  // no transfinite indications so mesh exactly follows the geometry
  Transfinite Line{_l1, _l2, _l3, _l4, _l5, _l6, _l7, _l8} = 1;
Return

// ~> create first physical surface
numberOfRectanglesInFirstPhysicalSurface = 5;
firstPhysicalSurfaceLength = 5.0;
rectangleLength = firstPhysicalSurfaceLength/numberOfRectanglesInFirstPhysicalSurface;
_xMin = 0.0;
_xMax = 0.0;
_yMin = 0.0;
_yMax = 1.0;
_physicalSurfaceName = 1;
Physical Surface(_physicalSurfaceName) = {};
For i In {0:numberOfRectanglesInFirstPhysicalSurface-1}
  _xMin = _xMax;
  _xMax = _xMin + rectangleLength;
  Call GenerateTriangulatedRectangle;
EndFor

// ~> create second physical surface next to the first one
numberOfRectanglesInSecondPhysicalSurface = 5;
secondPhysicalSurfaceLength = 5.0;
rectangleLength = secondPhysicalSurfaceLength/numberOfRectanglesInSecondPhysicalSurface;
_physicalSurfaceName = 2;
Physical Surface(_physicalSurfaceName) = {};
For i In {0:numberOfRectanglesInSecondPhysicalSurface-1}
  _xMin = _xMax;
  _xMax = _xMin + rectangleLength;
  Call GenerateTriangulatedRectangle;
EndFor

// ~~> set the mesh elements order
Mesh.ElementOrder = 1;

// ~> optimize mesh after generating it
Mesh.Optimize = 1;

// ~> generate 2D mesh
Mesh 2;

// ~> remove all duplicate mesh nodes
Coherence Mesh;

// ~> save mesh to file
Save "foo.msh";
