subdivision/subUtil/TetrahedronMesh.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file TetrahedronMesh.cpp
  \author Y. Lafranche
  \since 07 Dec 2007
  \date 08 Mar 2014

  \brief Implementation of xlifepp::subdivision::TetrahedronMesh class members and related functions
*/

#include "TetrahedronMesh.hpp"

using namespace std;

namespace xlifepp {
namespace subdivision {

//-------------------------------------------------------------------------------
//  I/O utilities
//-------------------------------------------------------------------------------

/*!
 Prints, on stream ftex, Fig4TeX instructions to draw the numbering convention
 for a tetrahedron mesh of order k.
 To be called as TetrahedronMesh(1,0,k,0).printTeXNumberConvention(ftex);
 */
void TetrahedronMesh::printTeXNumberConvention(ostream& ftex) const {
   if (subdiv_level_ > 0 || listT_.size() > 1) {
      cerr << "*** Error in printTeXNumberConvention:" << endl;
      cerr << "*** Numbering convention output function expects only one element.";
      cerr << endl;
   }
   int dim = max(5,int(order_+1));
   ftex << "\\input fig4tex.tex" << endl;
   ftex << "\\def\\drawFace#1#2#3{" << endl;
   ftex << "\\psset(color=\\FaceColor, fill=yes)\\psline[#1,#2,#3]" << endl;
   ftex << "\\psset(color=\\defaultcolor, fill=no)\\psline[#1,#2,#3,#1]}"<<endl;
   ftex << "%" << endl;
   ftex << "\\def\\defpts{" << endl;
   printTeXInArea(ftex,0,boundaryArea);
   ftex << "}" << endl;
   ftex << "% 1. Definition of characteristic points" << endl;
   ftex << "\\figinit{" << dim << "cm,orthogonal}" << endl;
   ftex << "\\defpts" << endl;
   ftex << "\\figset proj(psi=40, theta=45)" << endl;
   ftex << "%" << endl;
   ftex << "% 2. Creation of the graphical file" << endl;
   ftex << "\\psbeginfig{}" << endl;
   ftex << "\\figpt 0:(0,-1,0)\\psaxes 0(0.2)" << endl;
   for (number_t numBoundary=1; numBoundary<=3; numBoundary++) {
      printTeXFacesInArea(ftex,boundaryArea,numBoundary);
   }
   ftex << "\\psendfig" << endl;
   ftex << "%" << endl;
   ftex << "% 3. Writing text on the figure" << endl;
   ftex << "\\figvisu{\\figBoxA}{Points on the background faces.}{" << endl;
   ftex << "\\figsetmark{$\\figBullet$}\\figsetptname{{\\bf #1}}" << endl;
   for (number_t numBoundary=1; numBoundary<=3; numBoundary++) {
      printTeXInArea(ftex,1,boundaryArea,numBoundary);
   }
   ftex << "}" << endl;
// Last face of the tetrahedron
   ftex << "% 1. Definition of characteristic points" << endl;
   ftex << "\\figinit{4cm,orthogonal}" << endl;
   ftex << "\\defpts" << endl;
   ftex << "\\figset proj(psi=40, theta=45)" << endl;
   ftex << "%" << endl;
   number_t numBoundary4 = 4;
   ftex << "% 2. Creation of the graphical file" << endl;
   ftex << "\\psbeginfig{}" << endl;
   printTeXFacesInArea(ftex,boundaryArea,numBoundary4);
   ftex << "\\psendfig" << endl;
   ftex << "%" << endl;
   ftex << "% 3. Writing text on the figure" << endl;
   ftex << "\\figvisu{\\figBoxB}{Points on the foreground face (smaller scale).}{" << endl;
   ftex << "\\figsetmark{$\\figBullet$}\\figsetptname{{\\bf #1}}" << endl;
   printTeXInArea(ftex,1,boundaryArea,numBoundary4);
   ftex << "}" << endl;
   ftex << "\\centerline{\\box\\figBoxA\\hfil\\box\\figBoxB}" << endl;
// Internal vertices
   vector<number_t> VI = rk_verticesInside();
   if (VI.size() > 0)
   {
   ftex << "% 1. Definition of characteristic points" << endl;
   ftex << "\\figinit{" << dim << "cm,orthogonal}" << endl;
   ftex << "\\defpts" << endl;
   ftex << "\\figset proj(psi=40, theta=45)" << endl;
   ftex << "%" << endl;
   ftex << "% 2. Creation of the graphical file" << endl;
   ftex << "\\psbeginfig{}" << endl;
   for (number_t numBoundary=1; numBoundary<=3; numBoundary++) {
      printTeXFacesInArea(ftex,boundaryArea,numBoundary);
   }
   ftex << "\\psendfig" << endl;
   ftex << "%" << endl;
   ftex << "% 3. Writing text on the figure" << endl;
   ftex << "\\figvisu{\\figBoxA}{Internal points.}{" << endl;
   ftex << "\\figsetmark{$\\figBullet$}\\figsetptname{{\\bf #1}}" << endl;
   ftex << "% Vertices inside the boundaries" << endl;
   printTeXfigpt(ftex,VI);
   printTeXfigwrite(ftex,VI);
   ftex << "}" << endl;
   ftex << "\\smallskip\\centerline{\\box\\figBoxA}" << endl;
   }
   ftex << "\\medskip" << endl;
   ftex << "\\centerline{\\bf Nodes numbering convention for the tetrahedron of order "<< order_ << ".}" << endl;
   ftex << "%-------------------------------- End of figure --------------------------------" << endl;
   ftex << "\\vfill\\eject" << endl;
}


//-------------------------------------------------------------------------------
//  Protected member functions
//-------------------------------------------------------------------------------
/*!
 Subdivision of a tetrahedron T into 8 tetrahedrons.

 Input arguments :
 \param T         : tetrahedron to be subdivided

 Input and output arguments (updated by this function) :
 \param ElementNum: number of the last tetrahedron created in the mesh
 \param VertexNum : number of the last vertex created in the mesh
 \param listT     : list of tetrahedrons in the mesh
 \param SeenEdges : temporary map used to build the mesh (contains the edges already seen,
                    i.e. taken into account, along with the associated vertex number,
                    in order to avoid vertex duplication)

 Numbering and orientation convention
 - - - - - - - - - - - - - - - - - - -
 The vertices, faces and edges are numbered according to a convention given in Tetrahedron.hpp.
 The tetrahedron T=(1,2,3,4) is subdivided into 8 "sub-"tetrahedrons: 4 at each
 vertex plus 4 in the "kernel". Thus, the subdivision algorithm involves 6 more
 vertices. They are the midpoints of the edges of the tetrahedron T for every
 internal edge. For boundary edges, the additional points are also the midpoints
 if a flat subdivision is requested ; otherwise, the midpoints are projected onto
 the boundary surface. Their numbers follow the numbering convention given below:
 \verbatim
                   3                                         3
                  /.\                                       /.\
                 / . \                                     / . \
                /  .  \                                   /  .  \
               /   .   \                                 /   .   \
              /    .    \         SUBDIVISION           /    7    \
             /     .     \                             /     .     \
            /      .      \          ---->            /      .      \
           /       .       \                         9       .       8
          /      . 4 .      \                       /      . 4 .      \
         /     .       .     \                     /     .       .     \
        /    .           .    \                   /    .           .    \
       /   .               .   \                 /   .5             6.   \
      /  .                   .  \               /  .                   .  \
     / .                       . \             / .                       . \
    /.                           .\           /.                           .\
   /_______________________________\         /_______________10______________\
  1                                2        1                                2
 \endverbatim
 If, for instance, the points 1, 2 and 3 lie on the boundary and if a curved mesh (non
 flat) mesh is requested, then the points 8, 9 and 10 are projected onto the boundary
 surface.

 Each "sub-"tetrahedron is defined by a sequence of vertices that matches the above
 numbering convention.
 Nota: the orientation of the edges is not used in this subdivision algorithm.
 */
void TetrahedronMesh::algoSubdiv(const Tetrahedron& T, number_t& ElementNum,
                                 number_t& VertexNum, vector<Tetrahedron>& listT,
                                 map_pair_num& SeenEdges) {
// Creation of at most nb_edges_by_element_ new vertices: their rank in the list listV_ is retained in rV
   vector<number_t> rV(nb_edges_by_element_);
   for (number_t i=0; i<nb_edges_by_element_; i++) {
      rV[i]=createVertex(VertexNum,T.rkOfO1VeOnEdge(i+1),SeenEdges);}

// To help transmission of the curved boundary face number to the correct first 3 sub-tetrahedrons
// (which are defined with the same vertex ordering as T, so this boundary face number is the same):
   number_t bdSideNo[]={0,0,0,0}, bdfanum=T.bdSideOnCP();
   if (bdfanum > 0) {
      for (number_t i=0; i<3; i++) {
         bdSideNo[T.getrkFace(bdfanum-1,i)] = bdfanum;
      }
   }

// 1. Tetrahedrons at each vertex
// Vi is T.rankOfVertex(i) for i=1,...4 and is rV[i-5] for i=5,...10.
   listT.push_back(Tetrahedron(++ElementNum,T.rankOfVertex(1),rV[5],rV[4],rV[0],bdSideNo[0])); // ( 1,10,9,5)
   listT.push_back(Tetrahedron(++ElementNum,rV[5],T.rankOfVertex(2),rV[3],rV[1],bdSideNo[1])); // (10, 2,8,6)
   listT.push_back(Tetrahedron(++ElementNum,rV[4],rV[3],T.rankOfVertex(3),rV[2],bdSideNo[2])); // ( 9, 8,3,7)
   listT.push_back(Tetrahedron(++ElementNum,rV[0],rV[1],rV[2],T.rankOfVertex(4),bdSideNo[3])); // ( 5, 6,7,4)

// 2. Tetrahedrons inside the two pyramids
// Index arrays to create the following 4 tetrahedrons, according to 3 possible subdivisions:
//  s[0][t[i][j]] creates ( 5,8,10,9), ( 5,8,7, 6), (8, 5,7,9), (8, 5,10, 6), separation plane (5,10,8,7) ;
//  s[1][t[i][j]] creates ( 6,9, 5,7), ( 6,9,8,10), (9, 6,8,7), (9, 6, 5,10), separation plane (5, 6,8,9) ;
//  s[2][t[i][j]] creates (10,7, 6,8), (10,7,9, 5), (7,10,9,8), (7,10, 6, 5), separation plane (10,6,7,9).
   int t[4][4]={{0,1,3,2}, {0,1,5,4}, {1,0,5,2}, {1,0,3,4}};
   int s[3][6]={{0,3,4,5,1,2}, {1,4,2,0,5,3}, {5,2,3,1,0,4}};
// In order to obtain the 4 tetrahedrons, we choose the internal edge to be
// the shortest of the 3 diagonals 5-8, 6-9 or 7-10.
   real_t d[2]={listV_[rV[0]].squareDistance(listV_[rV[3]]),
          listV_[rV[1]].squareDistance(listV_[rV[4]])};
   int k = d[0] < d[1] ? 0 : 1;
   if (listV_[rV[2]].squareDistance(listV_[rV[5]]) < d[k]) k = 2;
// To help transmission of the curved boundary face number to the last sub-tetrahedron: we must identify
// which sub-tetrahedron is concerned and which face of this sub-tetrahedron is concerned.
// If bdfanum = 1 (resp. 2, 3, 4), the curved boundary subface is (6,7,8) (resp. (5,7,9), (5,6,10), (8,9,10))
// which belongs to one of the 4 tetrahedrons to be created, according to the subdivision k used. The
// correspondence is stored in the following index array, so that 0 <= u[k][bdfanum-1] <= 3 defines the right
// sub-tetrahedron in one of the 3 lists above:
//        subface (6,7,8)  (5,7,9), (5,6,10), (8,9,10)
   int u[3][4] = {{1,      2,       3,        0},
                  {2,      0,       3,        1},
                  {0,      1,       3,        2}};
// Due to the numbering used, the corresponding curved face number is 1 for all of them.
// Use this correspondence array to store the curved boundary face number for each tetrahedron to be created:
   number_t bdFaceNo[3][4] = {{0,0,0,0},{0,0,0,0},{0,0,0,0}};
   if (bdfanum > 0) {
      for (int ik=0; ik<3; ik++) {
         bdFaceNo[ik][u[ik][bdfanum-1]] = 1; // curved boundary face number is 1 for all
      }
   }
   for (int i=0; i<4; i++){
      listT.push_back(Tetrahedron(++ElementNum,rV[s[k][t[i][0]]],rV[s[k][t[i][1]]],
                                               rV[s[k][t[i][2]]],rV[s[k][t[i][3]]],bdFaceNo[k][i]));
   }
}

/*!
 Create high order vertices inside the tetrahedron T.
 On input, VertexNum is the previous number used. On output, it is the last one.
 */
void TetrahedronMesh::createHOiV(Tetrahedron& T, const number_t order, number_t& VertexNum){
   vector<Point> VP(nb_main_vertices_by_element_);
   vector<refnum_t> lc(nb_main_vertices_by_element_);
   for (number_t i=0; i<nb_main_vertices_by_element_; i++) {
      number_t rV = T.vertices_[i];
      VP[i] = listV_[rV].geomPt();
      lc[i] = listV_[rV].locCode();
   }
   refnum_t localcod = lc[0];
   for (number_t i=1; i<nb_main_vertices_by_element_; i++) { localcod &= lc[i]; }// localcod should be zero

   Point P;
   vector<real_t> coef(nb_main_vertices_by_element_);
   for (number_t i1=1; i1<order; i1++) {
      coef[0] = i1;
      for (number_t i2=1; i2<order-i1; i2++) {
         coef[1] = i2;
         for (number_t i3=1; i3<order-i1-i2; i3++) {
            coef[2] = i3;
            coef[3] = order-i1-i2-i3;
            P = barycenter(coef,VP);
           // add rank of the new vertex in the tetrahedron list:
            T.vertices_.push_back(VertexNum);
           // store the new vertex in the general list listV_:
            listV_.push_back(Vertex(++VertexNum,localcod,P));
         }
      }
   }
}


/*!
 Create the initial mesh (or "seed" of the mesh) to be subdivided for a geometry of revolution,
 i.e. defined by an axis of revolution and circular sections. The general shape is a truncated
 cone, that can be a cylinder if the radius is constant.
 In the following, the word "object" is used to designate the geometric shape.

 nbslices    : number of slices of elements orthogonal to the axis of the object
               If this input value is equal to 0, a number of slices is computed from
               the geometrical data.
 radius1, radius2 : radii of the object
 CharacPts   : the two end points of the axis of the object
 VertexNum   : number of the last vertex created (modified on output)
 ElementNum  : number of the last element created (modified on output)
 vSI         : vector for shape information at both ends of the object
 */
void TetrahedronMesh::initMesh(const number_t nbslices, const real_t radius1, const real_t radius2,
                               const vector<Point>& CharacPts, number_t& VertexNum,
                               number_t& ElementNum, const vector<ShapeInfo>& vSI){
   real_t R1, R2, slice_height;
   Point P1, P2;
   vector<ShapeInfo> cvSI;
   string bottomPt, topPt, shape;
   bool iscone;
   DefaultGeometry *DG;
   SurfPlane *SP;
   SurfCone *SC;
   vector <PatchGeometry *> EBS;
   int nb_sl;
   initRevMesh(nbslices, radius1,radius2, CharacPts, vSI,
               R1,R2, P1,P2, cvSI, bottomPt, topPt, iscone, shape, DG, SP, SC, EBS, nb_sl, slice_height);

   number_t iPh = 0, iPhSD;
   title_ = shape + " - Tetrahedron mesh";
   { // Definition of 3 boundaries...
      //vector<number_t> B_patch{1,2,3}, B_dimen{2,2,2}; Ready for C++11
      vector<number_t> B_patch(3, 1), B_dimen(3, 2);
      B_patch[1] = 2; B_patch[2] = 3;
      number_t nbBound = B_patch.size();
     // ... nbIntrf interfaces...
     // (two orthogonal internal planes containing the axis, a plane between
     //  two successive slices, plus eventually one plane between each end slice
     //  and the corresponding "hat" end shape)
      number_t nbIntrf = 2 + nb_sl-1, nbSbDom = nb_sl;
      if (EBS[0]->curvedShape()) {nbIntrf++; nbSbDom++;} // if non Flat
      if (EBS[1]->curvedShape()) {nbIntrf++; nbSbDom++;}
      vector<number_t> I_patch(nbIntrf), I_dimen(nbIntrf);
      number_t numPa = nbBound;
      for (size_t i=0; i<nbIntrf; i++) { I_patch[i] = ++numPa; I_dimen[i] = 2; }
     // ...and nbSbDom subdomains
     // (each slice plus eventually each "hat" end shape)
      vector<number_t> D_patch(nbSbDom), D_dimen(nbSbDom);
      for (size_t i=0; i<nbSbDom; i++) { D_patch[i] = ++numPa; D_dimen[i] = 3; }

      number_t nbBI =  nbBound + nbIntrf;
      number_t nbPatches =  nbBI + nbSbDom;
      vector<PatchGeometry *> P_geom(nbPatches);
      P_geom[0] = EBS[0], P_geom[1] = EBS[1];
      for (size_t i=2; i<nbBI; i++) { P_geom[i] = SP; }
      for (size_t i=nbBI; i<nbPatches; i++) { P_geom[i] = DG; }

      if (type_ == 0) {
         // P_geom = {EBS[0],EBS[1],SP, SP,...   SP, DG,...   DG}
         //           \__ nbBound  __/  \_nbIntrf_/  \_nbSbDom_/
         delete SC; // not needed
      }
      else {
         // P_geom = {EBS[0],EBS[1],SC, SP,...   SP, DG,...   DG}
         //           \__ nbBound  __/  \_nbIntrf_/  \_nbSbDom_/
         P_geom[2] = SC;
      }
      TG_ = TopoGeom(B_patch,B_dimen, I_patch,I_dimen, D_patch,D_dimen, P_geom);
      iPhSD = nbBound+nbIntrf;
   }
   refnum_t sBEs1, sBEs2, sBS, sIIp1, sIIp2; // localization code of the following patches
//       Description of the boundary patches
   TG_.setDescription(++iPh) = "Boundary: End surface on the side of end point " + bottomPt; sBEs1 = TG_.sigma(iPh); // patch 1
   TG_.setDescription(++iPh) = "Boundary: End surface on the side of end point " + topPt;    sBEs2 = TG_.sigma(iPh); // patch 2
   TG_.setDescription(++iPh) = "Boundary: Surface of the " + shape;                          sBS   = TG_.sigma(iPh); // patch 3
//       Description of the first two interface patches
   TG_.setDescription(++iPh) = "Interface: Internal plane 1 containing axis";                sIIp1 = TG_.sigma(iPh); // patch 4
   number_t iPhIIp1 = iPh;
   TG_.setDescription(++iPh) = "Interface: Internal plane 2 containing axis";                sIIp2 = TG_.sigma(iPh); // patch 5

// Steps 1 and 2 of the construction
   vector<Point> Pt;
   real_t prevRitrf;
   number_t rVaFirst, rVc;
   ElementNum = minElementNum_ - 1;
   VertexNum = minVertexNum_ - 1;
   startConeTrunk(R1, P1, EBS, SC, slice_height, nb_sl, iscone, sBEs1, sBS, sIIp1, sIIp2, iPhIIp1,
                  iPh, iPhSD, VertexNum, ElementNum, Pt, prevRitrf, rVaFirst, rVc);

// 3. Set of 4 vertices on top of last slice (on patch 2) and associated tetrahedra
//    The central point Pt[4] is inserted later.
   Vect AxV = SC -> AxisVector();
   for (size_t np=0; np<5; np++) { Pt[np] = translate(Pt[np],slice_height,AxV); }
      if (iscone) {
        Point& homcen(Pt[4]);
        // Compute the homothety ratio between 2 consecutive slices. homcen is the homothety center.
        real_t ratio = SC -> radiusAt(homcen) / prevRitrf;
        // Apply homothety to all the points in the interface, except the center.
        for (size_t np=0; np<4; np++) { Pt[np] = translate(homcen,ratio,Vect(homcen,Pt[np])); }
      }
   refnum_t sITF, sITFslice, sDOM;
   ++iPhSD, sDOM = TG_.sigma(iPhSD);    // Last slice
   if (EBS[1]->curvedShape()) { // if non Flat
      sITFslice = TG_.sigma(++iPh);     // Last interface between last slice and top end subdomain
      sDOM = sDOM | TG_.sigma(iPhSD+1); // Last slice and top end subdomain sharing this interface
   }
   else {
      sITFslice = 0;                    // No more transversal interface
   }
   for (size_t np=0; np<4; np++) {
      sITF = TG_.sigma(iPhIIp1 + np%2)|sITFslice; // one internal plane + transversal interface if any
      listV_.push_back(Vertex(++VertexNum,sBEs2|sBS|sITF|sDOM,Pt[np]));
   }
   // Associated tetrahedra
   number_t rVa = rVaFirst, rVb = rVa+1;
   subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum); rVa++; rVb++;
   subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum); rVa++; rVb++;
   subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum); rVa++; rVb=rVaFirst;
   subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum);

   if (EBS[1]->curvedShape()) { // if non Flat
// 4.a. Add tetrahedrons on the top
//      The central point on top of the last slice is internal.
         sITF = sIIp1|sIIp2|sITFslice; // center belongs to 3 interfaces
         listV_.push_back(Vertex(++VertexNum,           sITF|sDOM,Pt[4])); // rank rVc
      Point TopPt(EBS[1]->EndPt2()); // get Apex computed in endBoundaryShapes

      // Insert last vertex associated to top point in the global list
         sDOM = TG_.sigma(iPhSD+1); // Top end subdomain
         sITF = sIIp1|sIIp2; // two internal planes
         listV_.push_back(Vertex(++VertexNum,sBEs2     |sITF|sDOM,TopPt)); // rank rVd
      // Insert last 4 tetrahedra in the global list
      rVaFirst+=5; rVc+=5;
      number_t rVa = rVaFirst, rVb = rVa+1, rVd = rVc+1;
/*!
 \verbatim
                 d
                 |
             a+2 |                  d is the top point
              \  |                  a, b, c are in the top circular section of the object
               \ |                  c is the center
     a+3________\|_________b=a+1
                 \c
                  \
                   \
                    a
 \endverbatim
*/
      listT_.push_back(Tetrahedron(++ElementNum,rVa,rVb,rVd,rVc,4)); rVa++; rVb++;
      listT_.push_back(Tetrahedron(++ElementNum,rVa,rVb,rVd,rVc,4)); rVa++; rVb++;
      listT_.push_back(Tetrahedron(++ElementNum,rVa,rVb,rVd,rVc,4)); rVa++; rVb=rVaFirst;
      listT_.push_back(Tetrahedron(++ElementNum,rVa,rVb,rVd,rVc,4));
//      Description of the interface patch
      TG_.setDescription(iPh) = "Interface: Top section (containing end point 2)";
   }
   else {
// 4.b The central point of the top face lies on patch 2
         sITF = sIIp1|sIIp2; // two internal planes
         listV_.push_back(Vertex(++VertexNum,sBEs2      |sITF|sDOM,Pt[4]));
   }

//       Description of the subdomain patches
   if (EBS[0]->curvedShape()) { TG_.setDescription(++iPh) = "End subdomain on the side of end point " + bottomPt; }
   for (int isl=1; isl<=nb_sl; isl++) {
      ostringstream ss;
      ss << "Slice " << isl;
      TG_.setDescription(++iPh) = ss.str();
   }
   if (EBS[1]->curvedShape()) { TG_.setDescription(++iPh) = "End subdomain on the side of end point " + topPt; }
}

/*!
 This function makes the first part of the initial mesh (or "seed" of the mesh) to be subdivided
 for a geometry of revolution,
 i.e. defined by an axis of revolution and circular sections. The general shape is a truncated
 cone, that can be a cylinder if the radius is constant.
 In the following, the word "object" is used to designate the geometric shape.
 */
void TetrahedronMesh::startConeTrunk(real_t R1, const Point& P1, const vector <PatchGeometry *>& EBS,
                                     const SurfCone *SC, real_t slice_height, int nb_sl, bool iscone,
                                     refnum_t sBEs1, refnum_t sBS, refnum_t sIIp1, refnum_t sIIp2, number_t iPhIIp1,
                                     number_t& iPh, number_t& iPhSD, number_t& VertexNum, number_t& ElementNum,
                                     vector<Point>& Pt, real_t& prevRitrf, number_t& rVaFirst, number_t& rVc) {
// Define 5 "reference" points on the circle centered at the origin.
// Then, rotate and translate them so that they lie on the bottom face of the object.
    Pt.push_back(Point(R1,0,0));  // Pt[0]
    Pt.push_back(Point(0,R1,0));  // Pt[1]
    Pt.push_back(Point(-R1,0,0)); // Pt[2]
    Pt.push_back(Point(0,-R1,0)); // Pt[3]
    Pt.push_back(Point(0,0,0));   // Pt[4]
   Vect Z(0,0,1);
   Vect AxV = SC -> AxisVector();
   Vect Nor = crossProduct(Z,AxV); // rotation axis (Z and AxV are both unitary vectors)
   real_t st = norm(Nor);
   if (st > theTolerance){
   // the axis of the object is not Z: we need to rotate the points around the axis (Pt[4],Nor)
      Nor *= (1./st); // normalize Nor
      real_t ct = dot(Z,AxV);
      for (size_t np=0; np<4; np++) { Pt[np] = rotInPlane(Pt[np],ct,st,Pt[4],Nor); }
   }
   Vect OCP1(Pt[4],P1);
   for (size_t np=0; np<5; np++) { Pt[np] = translate(Pt[np],1.,OCP1); }

   refnum_t sITF, sITFslice;
   refnum_t sDOM = TG_.sigma(++iPhSD); // first subdomain : bottom end subdomain or first slice
   rVaFirst = 0, rVc = 4;
   if (EBS[0]->curvedShape()) { // if non Flat
// 1.a. Add tetrahedrons at bottom.
//      Set of 5 vertices at the bottom of the first slice whose center is internal, plus the bottom point.
      Point BottomPt(EBS[0]->EndPt2()); // get Apex computed in endBoundaryShapes
      // Insert first 6 vertices in the global list
         sITF = sIIp1|sIIp2; // two internal planes
         listV_.push_back(Vertex(++VertexNum,sBEs1    |sITF|sDOM,BottomPt)); // rank 0
      sITFslice = TG_.sigma(++iPh);     // first transversal interface
      sDOM = sDOM | TG_.sigma(iPhSD+1); // Bottom end subdomain and first slice
      for (size_t np=0; np<4; np++) {
         sITF = TG_.sigma(iPhIIp1 + np%2)|sITFslice; // one internal plane + transversal interface
         listV_.push_back(Vertex(++VertexNum,sBEs1|sBS|sITF|sDOM,Pt[np]));   // rank 1, 2, 3, 4
      }
         sITF = sIIp1|sIIp2|sITFslice; // center belongs to 3 interfaces
         listV_.push_back(Vertex(++VertexNum,          sITF|sDOM,Pt[4]));    // rank 5 (rVc)
         rVc++;
      // Insert first 4 tetrahedra in the global list
/*!
 \verbatim
             3
              \
               \                    0 is the bottom point
       4________\c_________2        1, 2, 3, 4, c are in the bottom circular section of the object
                 \                  c=5 is the center
                 |\
                 | \                Internal plane 1 is (0,1,3)
                 |  1               Internal plane 2 is (0,2,4)
                 |
                 0
 \endverbatim
*/
      listT_.push_back(Tetrahedron(++ElementNum,2,1,0,rVc,4));
      listT_.push_back(Tetrahedron(++ElementNum,3,2,0,rVc,4));
      listT_.push_back(Tetrahedron(++ElementNum,4,3,0,rVc,4));
      listT_.push_back(Tetrahedron(++ElementNum,1,4,0,rVc,4));
      rVaFirst++;
//      Description of the interface patch
      TG_.setDescription(iPh) = "Interface: Bottom section (containing end point 1)";
   }
   else {
// 1.b Set of 5 vertices at the bottom of the first slice (on patch 1)
      for (size_t np=0; np<4; np++) {
         sITF = TG_.sigma(iPhIIp1 + np%2); // one internal plane
         listV_.push_back(Vertex(++VertexNum,sBEs1|sBS|sITF|sDOM,Pt[np])); // rank 0, 1, 2, 3
      }
         sITF = sIIp1|sIIp2; // two internal planes
         listV_.push_back(Vertex(++VertexNum,sBEs1    |sITF|sDOM,Pt[4]));  // rank 4 (rVc)
      --iPhSD;
   }
// 2. Set of 5 vertices on top of each slice and associated tetrahedra
   // We consider the center point of the larger basis of the object (P1).
   // The images of this point by translations along the axis will be the centers of successive
   // homotheties in each interface plane.
   Point& homcen(Pt[4]);// warning : homcen is a reference to a point whose coordinates change
   prevRitrf = R1;
   for (int isl=1; isl<nb_sl; isl++,rVaFirst+=5,rVc+=5) {
      for (size_t np=0; np<5; np++) { Pt[np] = translate(Pt[np],slice_height,AxV); }
      if (iscone) {
        // Compute the homothety ratio between 2 consecutive slices. homcen is the homothety center.
        real_t Ritrf = SC -> radiusAt(homcen), ratio = Ritrf / prevRitrf;
        prevRitrf = Ritrf;
        // Apply homothety to all the points in the interface, except the center.
        for (size_t np=0; np<4; np++) { Pt[np] = translate(homcen,ratio,Vect(homcen,Pt[np])); }
      }
      sITFslice = TG_.sigma(++iPh);                          // Next transversal interface
      ++iPhSD, sDOM = TG_.sigma(iPhSD) | TG_.sigma(iPhSD+1); // Two slices sharing this interface
      for (size_t np=0; np<4; np++) {
         sITF = TG_.sigma(iPhIIp1 + np%2)|sITFslice; // one internal plane + transversal interface
         listV_.push_back(Vertex(++VertexNum,sBS|sITF|sDOM,Pt[np]));
      }
         sITF = sIIp1|sIIp2|sITFslice; // center belongs to 3 interfaces
         listV_.push_back(Vertex(++VertexNum,    sITF|sDOM,Pt[4]));
      // Associated tetrahedra
      number_t rVa = rVaFirst, rVb = rVa+1;
      subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum); rVa++; rVb++;
      subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum); rVa++; rVb++;
      subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum); rVa++; rVb=rVaFirst;
      subdivPrism(rVa,rVb,rVc,rVa+5,rVb+5,rVc+5, ElementNum);
//      Description of the interface patch
      ostringstream ss;
      ss << "Interface: Section " << isl;
      TG_.setDescription(iPh) = ss.str();
   }
}

/*!
 Create the tetrahedrons that subdivide a prism.
 A prism can be subdivided in 3 tetrahedrons in six ways. We must take care to
 have compatible subdivisions in prisms put together side by side. Here, since
 we restrict to the case where neighbour prisms share the edge (3,6), the compa-
 tibility is achieved with 4 subdivisions among the six possible. This function
 selects one of those.

 Ui, i=1,...6 are the ranks of the vertices in the general vertex list.
 The vertices must be ordered as shown on the following figure.
 ElementNum is the number of the last tetrahedron created and is updated here.

 The face (1,2,5,4) is on the boundary of the object and thus corresponds to
 a curved patch. This is used to set the face number of the tetrahedrons lying
 on this patch.
 \verbatim
                            5 +
                            / | \
                          /   |  \
                        /     |   \             Subdivision :
                      /       |    \               1 2 6 3
                    /         |     \              2 4 6 5
                6 +-----------|------+ 4           2 6 4 1
  Z               |           |      |
  ^     Y         |         2 +      |
  |   /           |         /   \    |
  | /             |       /      \   |
  +-----> X       |     /         \  |
                  |   /            \ |
                  | /               \|
                3 +------------------+ 1
 \endverbatim
 */
void TetrahedronMesh::subdivPrism(const number_t U1, const number_t U2, const number_t U3,
                                  const number_t U4, const number_t U5, const number_t U6,
                                  number_t& ElementNum){
   listT_.push_back(Tetrahedron(++ElementNum,U1,U2,U6,U3));// internal
   listT_.push_back(Tetrahedron(++ElementNum,U2,U4,U6,U5,3));// face 3 of this tetrahedron on the boundary
   listT_.push_back(Tetrahedron(++ElementNum,U2,U6,U4,U1,2));// face 2 of this tetrahedron on the boundary
}

/*!
 Define some macros in the file associated to the stream ftex in order to display
 the mesh using Fig4TeX.
 */
void TetrahedronMesh::printTeXHeader(ostream& ftex) const {
   ftex << "\\def\\drawFace#1#2#3#4{" << endl;
   ftex << "\\figset(color=#4, fill=yes)\\figdrawline[#1,#2,#3]" << endl;
   ftex << "\\figset(color=default, fill=no)\\figdrawline[#1,#2,#3,#1]}" << endl;
   ftex << "\\def\\drawElem#1#2#3#4{" << endl;
   ftex << "\\figdrawline[#1,#2,#3,#1,#4,#2]" << endl;
   ftex << "\\figdrawline[#4,#3]}" << endl;
}

} // end of namespace subdivision
} // end of namespace xlifepp