/**
 * @brief    Minimal Yet Complete example a problem with the OpenCASCADE
 *           factory in GMSH when performing a full 360 degree extrusion.
 * @author   Sigve Karolius (SK)
 * @date     2018-08-20
 * @details
 *
 * The purpose if to generate a 3D structured hexahedral mesh of a hollow pipe.
 * Please note that this is just an illustration of the problem. This specific
 * mesh can be generated in a number of other ways but the actual geometry
 * requires an analogous approach.
 *
 * The file consists of four parts:
 *
 * 1. Defines relevant geometry/mesh parameters and displays them in the GUI
 * 2. Defines a macro for generating a rectangle that is to be extruded about
 *    about an axis effectively forming a hollow pipe.
 * 3. Defines a macro for extruding a "surface" four times about and axis
 * 4. Calls the aforementioned macros and the Mesh command
 *
 * By opening the file in the GMSH (v. 3.x or greater) and changing the tag
 * "Use OpenCASCADE" in the menu on the right hand side from "No" to "Yes" the
 * error appears.
*/

// ========================================================================= //
// Define Parameters Governing Geometry/Mesh Generation
// ========================================================================= //
R = 1;                                                                // Radius
dR = 0.5;                                                    // Shell thickness
cx = 0; cy = 0; cz = 0;                      // Coordinates of the center point
ex = 1; ey = 0; ez = 0;                                 // Rotate around x axis
L    = 2;                                                   // Length of "pipe"
Nr0 = 5;                                     // Volumes in the radial direction
Na0 = 3;                                    // Volumes in the angular direction
Nl0 = 7;                                     // Volumes in the length direction
m = 1;            // Scaling factor to change (SI) units (not implemented here)

DefineConstant [
  radius    = {   R*m, Name "Input/Characteristic Sizes/Geometry/Radius (m)"}
  thickness = {  dR*m, Name "Input/Characteristic Sizes/Geometry/Thickness (m)"}
  len       = {   L*m, Name "Input/Characteristic Sizes/Geometry/Length (m)"}
  Nlen      = {   Nl0, Name "Input/Characteristic Sizes/Mesh/Layers"}
  Nrad      = {   Nr0, Name "Input/Characteristic Sizes/Mesh/Radial Layers"}
  Nang      = {   Na0, Name "Input/Characteristic Sizes/Mesh/Angular Layers"}
  CASCADE   = {0, Choices{0="No", 1="Yes"},
                Name "Input/0Use OpenCASCADE", Highlight "Red"}
];

If( CASCADE == 1)
  SetFactory("OpenCASCADE");
EndIf

// ========================================================================= //
// Macro Generating a Rectangle
// ========================================================================= //

Function DeclareRectangularSurface
  /**
   *
   *       p2 +-------------+ p3
   *          |  "surface"  |
   *       p1 +-------------+ p4
   *
   * (center) *
   *
  */

  /* Define coordinates for the points relative to the center of the pipe */
  P1x = cx; P1y = cy;  P1z = cz;
  P2x = cx; P2y = cy;  P2z = cz;
  P3x = cx; P3y = cy;  P3z = cz;
  P4x = cx; P4y = cy;  P4z = cz;

              P1y += thickness;
              P2y += radius;
  P3x += len; P3y += radius;
  P4x += len; P4y += thickness;

  /* Define Points */
  p1 = newp; Point(p1) = {P1x, P1y, P1z, 1.0};
  p2 = newp; Point(p2) = {P2x, P2y, P2z, 1.0};
  p3 = newp; Point(p3) = {P3x, P3y, P3z, 1.0};
  p4 = newp; Point(p4) = {P4x, P4y, P4z, 1.0};

  /* Define Topology */
  l1 = newl; Line(l1) = {p1, p2};                                  // p1 --> p2
  l2 = newl; Line(l2) = {p2, p3};                                  // p2 --> p3
  l3 = newl; Line(l3) = {p3, p4};                                  // p3 --> p4
  l4 = newl; Line(l4) = {p4, p1};                                  // p4 --> p1

  ll1 = newll; Line Loop(ll1) = {l1, l2, l3, l4}; // p1 -> p2 -> p3 -> p4 -> p1

  surface = news; Plane Surface(surface) = {ll1};

  /* Mesh Generation */
  Transfinite Line{l1, l3} = Nrad + 1 Using Bump 0.1;
  Transfinite Line{l2, -l4} = Nlen + 1 Using Progression 0.9;

  Transfinite Surface{surface};
  Recombine Surface{surface};

Return

// ========================================================================= //
// Macros Performing Extrusions
// ========================================================================= //

Function ExtrudeSurfaceFourTimesAboutAxis
  /*
   * Rotates the variable "surface" (see macro above) four times to form a
   * circle.
   *
   *          0
   *          +
   *          |
   *          +
   * 3  +--+  *  +--+ 1
   *          +
   *          |
   *          +
   *          2
  */
  // Define lists for labeling purposes
  volumes     = {};
  innerfaces  = {};
  outerfaces  = {};
  inletfaces  = {};
  outletfaces = {};

  // Perform Extrusion
  For i In {0:3}
    out[] = Extrude{ {ex, ey, ez}, {cx, cy, cz}, Pi/2 } {
      Surface{surface};
      Layers{Nang};
      Recombine;
    };
    surface      = out[0];
    volumes     += out[1];
    inletfaces  += out[2];
    outerfaces  += out[3];
    outletfaces += out[4];
    innerfaces  += out[5];
  EndFor

Return

// ========================================================================= //
// Generate and Mesh Geometry (for completeness we also add labels)
// ========================================================================= //
Call DeclareRectangularSurface;                       // Create patch "surface"
Call ExtrudeSurfaceFourTimesAboutAxis;     // Rotate patch "surface" four times

Physical Surface("inlets")  = inletfaces[];
Physical Surface("outlets") = outletfaces[];
Physical Surface("walls")   = {outerfaces[], innerfaces[]};
Physical Volume("volume")   = volumes[];

Mesh 3;

// Vim 'modelines', standardise settings for the 'vim' text editor...
// vim: filetype=gmsh fileencoding=ascii syntax=on colorcolumn=79
// vim: ff=unix tabstop=2 softtabstop=0 expandtab shiftwidth=2 smarttab
