xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/geometry/subdivision/subUtil/SubdivisionMesh.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 SubdivisionMesh.cpp
  \author Y. Lafranche
  \since 22 May 2008
  \date 08 Mar 2014

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

#include "SubdivisionMesh.hpp"
#include "TeXutil.hpp"
#include "PointUtils.hpp"

using namespace std;

namespace xlifepp {
namespace subdivision {
//-------------------------------------------------------------------------------
//  Static variables
//-------------------------------------------------------------------------------
//HOV_error SubdivisionMesh::PB;

//-------------------------------------------------------------------------------
//  Constructors, Destructor
//-------------------------------------------------------------------------------
/*!
 Main constructor.
 */
SubdivisionMesh::SubdivisionMesh(const number_t nbsubdiv, const number_t order,
         const number_t type, const number_t minVertexNum, const number_t minElementNum)
: subdiv_level_(nbsubdiv), order_(order), type_(type),
  minVertexNum_(minVertexNum), minElementNum_(minElementNum),
  newVertexPt_(&SubdivisionMesh::newVertexPtDef) {
   if (type != 0 ) { newVertexPt_ = &SubdivisionMesh::newVertexPtGen; }
}

/*!
 Destructor.
 Free memory allocated for patch shapes taking care of duplicated pointers
 (see IMPORTANT NOTE in TopoGeom.hpp).
 */
SubdivisionMesh::~SubdivisionMesh(){
   for (number_t numPatch=1; numPatch<=TG_.numberOfPatches(); numPatch++) {
      bool notAlreadySeen = true;
      for (number_t j=1; j<numPatch; j++) {
         if (TG_.shape(numPatch) == TG_.shape(j)) {
            notAlreadySeen = false;
            break;
         }
      }
      if (notAlreadySeen) {delete TG_.shape(numPatch);}
   }
}
//-------------------------------------------------------------------------------
//  Public Access functions
//-------------------------------------------------------------------------------
/*!
 Returns the coordinates of a vertex, given its number num, num >= minVertexNum_
 */
vector<real_t> SubdivisionMesh::vertexCoord(const number_t num) const{
   return listV_.at(num-minVertexNum_).geomPt();
}
/*!
 Returns the coordinates of a vertex, given its rank rk>=0 in the global list of vertices
 */
vector<real_t> SubdivisionMesh::rkvertexCoord(const number_t rk) const{
   return listV_.at(rk).geomPt();
}

int DECHOL(const vector<vector<real_t> >& A, int n, vector<vector<real_t> >& L, real_t eps);
void DRCHOL(const vector<vector<real_t> >&L, int n, vector<vector<real_t> >&B, int m);
void IMP(const vector< vector<number_t> >& comesh,
         const vector<pair<short,short> >& rkEV,
         const vector<number_t>& rkUnk, const vector<number_t>& rkData) {
   cout << " MAILLAGE : " << endl;
   for (vector< vector<number_t> >::const_iterator itcm=comesh.begin(); itcm<comesh.end(); itcm++) {
      for (vector<number_t>::const_iterator itt=itcm->begin(); itt<itcm->end(); itt++) {
         cout << *itt << " ";
      }
      cout << endl;
   }
   cout << " rkEV :" << endl;
   for (vector<pair<short,short> >::const_iterator itEV=rkEV.begin(); itEV<rkEV.end(); itEV++) {
      cout << itEV->first << " " << itEV->second << endl;
   }
   cout << " rkUnk :" << endl;
   for (vector<number_t>::const_iterator itt=rkUnk.begin(); itt<rkUnk.end(); itt++) { cout << *itt << " "; }
   cout << endl;
   cout << " rkData :" << endl;
   for (vector<number_t>::const_iterator itt=rkData.begin(); itt<rkData.end(); itt++) { cout << *itt << " "; }
   cout << endl;
}
// connectivity mesh, ranks of edge vertices, ranks of unknowns, ranks of data
vector<Point> SubdivisionMesh::unifMesh(const vector< vector<number_t> >& comesh,
                                        const vector<pair<short,short> >& rkEV,
                                        const vector<number_t>& rkUnk, const vector<number_t>& rkData) {

//IMP(comesh, rkEV, rkUnk, rkData);
   number_t dimM=rkUnk.size(), tN=2*dimM, fN=4*dimM;
// Associate a row index to each unknown number and tN to all data number
   map<number_t, number_t> row;
   vector<number_t>::const_iterator itv;
   number_t i=0;
   for (itv=rkUnk.begin(); itv<rkUnk.end(); itv++) { row[*itv] = i++; }
   for (itv=rkData.begin(); itv<rkData.end(); itv++) { row[*itv] = tN; }

// Build the matrix M (lower triangular and initialized to 0)
   vector<vector<real_t> > M(dimM);
   for (number_t p=0; p<dimM; p++) { M[p].assign(p+1, 0.); }
   vector<set<number_t> > S(dimM);
   // for each element in the mesh
   for (vector< vector<number_t> >::const_iterator itcm=comesh.begin(); itcm<comesh.end(); itcm++) {
      // for each edge AiAj
      for (vector<pair<short,short> >::const_iterator itEV=rkEV.begin(); itEV<rkEV.end(); itEV++) {
         number_t i=(*itcm)[itEV->first], j=(*itcm)[itEV->second];
         number_t p=row[i], q=row[j], qJ=j, ppq=p+q;
         if (ppq < tN) { // Ai and Aj are both unknown points
            if (p < q) { number_t r=p; p=q; q=r; qJ=i; }// swap indices
            // Now, we have p > q and p cannot be equal to q here.
            // Update line p of the matrix
            if (M[p][q] >= 0) { M[p][q] = -1.; M[p][p] += 1.; M[q][q] += 1.; }
            }
         else {
            if (ppq < fN) { // One point is unknown and the other data.
               if (q < p)  { number_t r=p; p=q; q=r; qJ=i; }// swap indices
               // Now, AqJ is a data point. Store it to compute the rhs later.
               pair<set<number_t>::iterator,bool> res = S[p].insert(qJ); // qJ is inserted only once in the set
               if (res.second) { M[p][p] += 1.; } // update line p of the matrix only once
            }
            else { continue; }// ppq == fN : Ai and Aj are both data points (AiAj boundary edge)
         }
      }
   }
// Build the rhs
   dimen_t dimpt=rkvertexCoord(comesh[0][0]).size(); // space dimension
   vector<vector<real_t> > B(dimM,vector<real_t>(dimpt, 0.));
   // Update the rhs according to S'contents. B is initialized with 0 which is important
   // for values of p such that S[p] is empty.
   for (number_t p=0; p<dimM; p++) {
      for (set<number_t>::const_iterator itSp=S[p].begin(); itSp!=S[p].end(); itSp++) {
         vector<real_t> Coor = rkvertexCoord(*itSp);
         for (number_t k=0; k<dimpt; k++) { B[p][k] += Coor[k]; }
      }
   }
// Solve the system and return the result
   if ( DECHOL(M,dimM,M,theTolerance) ) { error("mat_noinvert"); }
   else { DRCHOL(M,dimM,B,dimpt); }
   vector<Point> vPts(dimM, Point(std::vector<real_t>(dimpt)));
   for (number_t p=0; p<dimM; p++) { vPts[p] = B[p]; }
   return vPts;
}

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

/*!
 Prints, on stream os, the general informations about the mesh.
 */
void SubdivisionMesh::printInfo(ostream& os, const bool forTeX) const {
   string prefix = forTeX ? "\\" : "" ;
   os << prefix << "   "; for (size_t i=0; i<title_.size(); i++) {os << "=";} os << endl;
   string title = forTeX ? fmtTeX(title_) : title_;
   os << prefix << "   " << title << endl;
   os << prefix << "   "; for (size_t i=0; i<title_.size(); i++) {os << "=";} os << endl;
   os << prefix << " Size: " << numberOfElements() << " elements, "
      << listV_.size() << " vertices" << endl;
   os << prefix << " Order " << order_ << ", subdivision level " << subdiv_level_ << endl;
   if (type_ > 0) {
      os << prefix << " Curved mesh";
      if (order_ > 2) {
         if (type_ == 1) {os << " with radial projection algorithm";}
         else if (type_ == 2) {os << " with rotation algorithm";}
      }
      os << endl;
   }
   TG_.print(os,forTeX);
}

//-------------------------------------------------------------------------------
//  Other public member functions
//-------------------------------------------------------------------------------
/*!
 Initializes the attribute of the areas of the mesh. The areas are considered in
 the order boundary, then interface and at last subdomain.
 If the input data vector holds less attributes than the number of areas, the
 data are reused cyclically starting again from the beginning.
 */
void SubdivisionMesh::initUserAttribute(const vector<string>& attribute) {
   size_t  k = 0;
   initAttribute_(boundaryArea,attribute,k);
   initAttribute_(interfaceArea,attribute,k);
   initAttribute_(subdomainArea,attribute,k);
}
/*!
 Initializes the attribute of the areas of kind TA of the mesh from vector attribute.
 Called by initUserAttribute. k is the current index in vector attribute.
 The vector attribute is scanned cyclically until all the areas have been assigned
 an attribute.
 */
void SubdivisionMesh::initAttribute_(const topologicalArea TA, const vector<string>& attribute, size_t& k){
   number_t nbAreas = TG_.numberOf(TA);
   number_t nbAttrib = attribute.size();
   number_t ica=0;

   while (ica<nbAreas) {
      number_t nb1 = nbAttrib-k, nb2 = nbAreas-ica;
      if (nb1 >= nb2) {
         for (number_t i=0; i<nb2; i++,ica++){
            TG_.setAttribute(TA,ica+1) = attribute[k++];
         }
      }
      else {
         for (number_t i=0; i<nb1; i++,ica++){
            TG_.setAttribute(TA,ica+1) = attribute[k++];
         }
         k=0;
      }
   }
}
/*!
 Returns the number of the internal vertices of any order in the mesh
 */
number_t SubdivisionMesh::numberOfVerticesInside() const {
   return rk_verticesInside().size();
}
/*!
 Returns the number of vertices of any order in the mesh, belonging to the area
 num (num>=1) of kind TA or belonging to any area of kind TA if num=0.
 */
number_t SubdivisionMesh::numberOfVerticesIn(const topologicalArea TA, const number_t num) const {
   return rk_verticesIn(TA,num).size();
}
/*!
 Returns the list of vertices of any order in a tetrahedron mesh, which are
 not on the boundaries. Thus, they are inside the domain.
 Vertices are defined by their number (starting from 1).
 */
vector<number_t> SubdivisionMesh::verticesInside() const {
   vector<number_t> VI = rk_verticesInside();
   rankToNum(VI);
   return VI;
}
/*!
 Returns the list of vertices of any order in a tetrahedron mesh, belonging to
 the area num (num>=1) of kind TA or belonging to any area of kind TA if num=0.
 Vertices are defined by their number (starting from 1).
 */
vector<number_t> SubdivisionMesh::verticesIn(const topologicalArea TA, const number_t num) const {
   vector<number_t> VonB = rk_verticesIn(TA,num);
   rankToNum(VonB);
   return VonB;
}

/*!
 Returns the list of vertices of order 1 in a mesh of simplices.
 Vertices are defined by their number (starting from 1).
 */
vector<number_t> SubdivisionMesh::verticesOfOrder1() const {
   vector<number_t> V(nb_main_vertices_);

   for (size_t i=0; i<nb_main_vertices_; i++) {
    V[i] = listV_[i].number();
   }
   return V;
}

//-------------------------------------------------------------------------------
//  Protected member functions
//-------------------------------------------------------------------------------
/*!
 Create a new vertex starting from a list of vertices and insert it in the general vertex list.
 The function returns the rank (in the general vertex list) of the new vertex created.
 The vertices are those of an element or a face or an edge.
 The new vertex is computed from the isobarycenter of the given vertices (it is this
 isobarycenter or its projection onto a boundary surface).
 Input arguments :
 \param VertexNum   : number of the last vertex created
 \param rkVert      : list of vertices being processed (rank of the vertices
                       in the general vertex list)
 Output arguments :
 \param VertexNum   : number of the new vertex created
 */
number_t SubdivisionMesh::createVertex(number_t& VertexNum, const vector<number_t>& rkVert){
   vector<Point> VP(rkVert.size());
   real_t* coef = new real_t[rkVert.size()];
   vector<number_t>::const_iterator itrV = rkVert.begin();
   number_t ip(0);
   Vertex& V = listV_[*itrV];
   refnum_t localcod = V.locCode();
   VP[ip] = V.geomPt();
   coef[ip++] = 1.;
   while (++itrV < rkVert.end()) {
      Vertex& V = listV_[*itrV];
      localcod &= V.locCode();
      VP[ip] = V.geomPt();
      coef[ip++] = 1.;
   }
   Point P = (this->*newVertexPt_)(localcod,coef,VP);
   number_t num = listV_.size();
   listV_.push_back(Vertex(++VertexNum,localcod,P));

   delete [] coef;
   return num;
}

/*!
 Create a new vertex starting from the two endpoints V1 and V2 of an edge.
 The function returns the rank (in the general vertex list) of the new vertex
 or the rank of the corresponding vertex if it has already been created.
 Input arguments :
 \param VertexNum   : number of the last vertex created
 \param Edge        : edge being processed
 \param SeenEdges   : list of the edges already processed
 Output arguments :
 \param VertexNum   : if a new vertex was created, number of this vertex,
                      unchanged otherwise
 \param SeenEdges   : list of the edges already processed, updated with the current
                      one if a new vertex was created
 */
number_t SubdivisionMesh::createVertex(number_t& VertexNum,
                                      const pair_nn Edge, map_pair_num& SeenEdges){
   const number_t rV1=Edge.first, rV2=Edge.second;
   map_pair_num::iterator itSE;
   number_t num;

   if ((itSE=SeenEdges.find(Edge)) == SeenEdges.end()) {
// This edge has not yet been seen: compute the corresponding vertex
      Vertex V1=listV_[rV1];
      Vertex V2=listV_[rV2];
      refnum_t localcod = V1.locCode() & V2.locCode();
      vector<Point> VP(2);
      VP[0] = V1.geomPt();
      VP[1] = V2.geomPt();
      real_t coef[2] = {1., 1.};
      Point P;
      P = (this->*newVertexPt_)(localcod,coef,VP);
      num = listV_.size();
// A new vertex is being created on this edge: the edge is stored along with
// the rank of the new vertex. The edge is stored in the 2 forms (V1,V2) and
// (V2,V1) for simplicity reason, to make its future retrieval easier.
      SeenEdges.insert(make_pair(Edge,num));
      SeenEdges.insert(make_pair(make_pair(rV2,rV1),num));
      listV_.push_back(Vertex(++VertexNum,localcod,P));
   }
   else {
// Edge found: get back the rank of the corresponding vertex in the general list
      num = itSE->second;
   }
   return num;
}
/*!
 Create a new vertex starting from the vertices of a face.
 The function returns the rank (in the general vertex list) of the new vertex
 or the rank of the corresponding vertex if it has already been created.
 Input arguments :
 \param VertexNum   : number of the last vertex created
 \param rkVFace     : list of vertices defining the face being processed (rank of the vertices
                      in the local list of the element)
 \param SeenFaces   : list of the faces already processed
 Output arguments :
 \param VertexNum   : if a new vertex was created, number of this vertex,
                      unchanged otherwise
 \param SeenFaces   : list of the faces already processed, updated with the current
                      one if a new vertex was created
 */
number_t SubdivisionMesh::createVertex(number_t& VertexNum,
                                       const vector<number_t>& rkVFace, map_set_num& SeenFaces){
   set_n Face(rkVFace.begin(),rkVFace.end());
   map_set_num::iterator itSF;
   number_t num;

   if ((itSF=SeenFaces.find(Face)) == SeenFaces.end()) {
// This face has not yet been seen: compute the corresponding new vertex
      num = createVertex(VertexNum,rkVFace);
// A new vertex is being created on this face: the face is stored along with
// the rank of the new vertex.
      SeenFaces.insert(make_pair(Face,num));
   }
   else {
// Face found: get back the rank of the corresponding vertex in the general list
      num = itSF->second;
   }
   return num;
}
   //--- Utilitary functions

/*!
 Returns the reference number (localization code) of the area num or of the union
 of the areas if num=0. The kind of area is given by TA ; it may be boundary,
 interface or subdomain.
 */
refnum_t SubdivisionMesh::lCodeOf(const topologicalArea TA, const number_t num) const {
   if (num > 0) {
      return TG_.localCodeOf(TA,num);
   }
   else {
      return TG_.maskOf(TA);
   }
}
/*!
 Returns color associated to an area of kind TA or default color.
 The number of the area num must be >= 1.
 */
string SubdivisionMesh::colorOf(const topologicalArea TA, const number_t num) const {
   string defaultColor("0.5"); // medium gray
   string Color = defaultColor;
   if (num <= TG_.numberOf(TA)) {
      Color = TG_.getAttribute(TA,num);
      if (Color.empty()) Color = defaultColor;
   }
   return Color;
}

/*!
 Initializes the user attribute of the geometrical areas of the mesh,
 boundaries, interfaces and subdomains, given in this order, by a color.
 */
void SubdivisionMesh::initDefaultUserAttribute(void){
   vector <string> color;

   color.push_back("\\DarkOrangergb");
   color.push_back("\\ForestGreenrgb");
   color.push_back("\\Cyanrgb");
   color.push_back("\\Yellowrgb");
   color.push_back("\\Magentargb");
   color.push_back("\\Bluergb");
   color.push_back("\\Redrgb");
   color.push_back("\\Firebrickrgb");
   color.push_back("\\Pinkrgb");
   color.push_back("\\DarkGoldenrodrgb");
   color.push_back("\\Goldrgb");
   color.push_back("\\HotPinkrgb");
   color.push_back("\\Maroonrgb");
   color.push_back("\\RoyalBluergb");

   initUserAttribute(color);
}
/*!
   Apply rotations defined in rots and the translation of vector U to each point of Pt
*/
// Returns the rotation matrix in R3 of angle theta around the absolute axis number 1, 2 or 3.
// The angle is to be given in degrees.
Matrix<real_t> rotationMatrix(real_t theta, dimen_t axis) {
  dimen_t nbrows(3);
  Matrix<real_t> rotmat(nbrows,nbrows, 0.);
  dimen_t i(axis%nbrows), j((i+1)%nbrows), k(axis-1);
  real_t trad(theta*pi_/180.), st(std::sin(trad)), ct(std::cos(trad));
  rotmat[i*nbrows+i] = ct;
  rotmat[j*nbrows+i] = st;
  rotmat[i*nbrows+j] = -st;
  rotmat[j*nbrows+j] = ct;
  rotmat[k*nbrows+k] = 1.;
  return rotmat;
}
void SubdivisionMesh::rotNtrans(const vector<pair<real_t, dimen_t> >& rots, const Vect& U, vector<Point>& Pt) {
  vector<pair<real_t, dimen_t> >::const_iterator it_rots;
  for (it_rots=rots.begin(); it_rots!=rots.end(); it_rots++) {
     Matrix<real_t> R(rotationMatrix(it_rots->first, it_rots->second));
     for (size_t np=0; np<Pt.size(); np++) {
        Vector<real_t> v(Pt[np].toVect());
        Pt[np] = Point(R*v);
     }
  }
  for (size_t np=0; np<Pt.size(); np++) { Pt[np] = translate(Pt[np],1.,U); }
}

/*!
 Return the boundary shapes to be used at both end of an object of revolution (a cylinder, a cone
 or a truncated cone), according to the data provided by the input arguments:
 \param vSI : information needed to define the boundary shapes and compute their
              intersection with the axis of the object (apexes).
 \param SC  : definition of the object of revolution
 */
vector<PatchGeometry *> SubdivisionMesh::endBoundaryShapes(const vector<ShapeInfo>& vSI,
                                                           const SurfCone& SC){
   vector <PatchGeometry *> EBS;
   Vect AxV = SC.AxisVector();
   vector<real_t> R(SC.radii());
   vector<real_t>::const_iterator itR=R.begin();
   vector<Point> EndPts;
   EndPts.push_back(SC.EndPt1());
   EndPts.push_back(SC.EndPt2());
   vector<Point>::const_iterator itEP=EndPts.begin();
   vector<ShapeInfo>::const_iterator itvSI;
   int direction = 1, shape;
   for (itvSI=vSI.begin(); itvSI != vSI.end(); itvSI++, itEP++, itR++) {
      direction = -direction; // first translate P1 in the direction of P2P1, then P2 in the direction of P1P2
      Point Apex=translate(*itEP,direction * itvSI->dist_, AxV);
      if (itvSI->dist_ > theTolerance) { shape = itvSI->shapeCode_; }
      else { shape = Flat; }// discard the shape if degenerated
      switch (shape) {
         case Cone:
            EBS.push_back(new SurfCone(*itEP,Apex,*itR));
            break;
         case Ellipsoid:
            EBS.push_back(new SurfEllipsoid(*itEP,Apex,*itR));
            break;
         case Sphere: {
            // Compute PtOnAxis, intersection of the sphere and the axis of the object
            Point PtOnAxis=translate(*itEP,direction * *itR, AxV);
            EBS.push_back(new SurfSphere(*itEP,PtOnAxis,*itR));
            }
            break;
         default: // Flat
            EBS.push_back(new SurfPlane());
            break;
      }
   }
   return EBS;
}

/*!
  Defines a set of points used to build the initial mesh in a cube.
*/
vector<Point> SubdivisionMesh::cubePoints(const real_t edLen, const vector<pair<real_t, dimen_t> >& rots, const Point& Center) {
// Define the 27 characteristic points of the cube near the origin Pt[4] (corners,
// middle of the edges and faces) in order to define the selected octants.
   vector<Point> Pt;
   real_t halfEL = edLen / 2;
   Pt.push_back(Point(      0,      0,      0)); // unused                      7
   Pt.push_back(Point( halfEL,      0,      0)); // Pt[1]     Z>0     11+-------+-------+6
   Pt.push_back(Point( halfEL, halfEL,      0)); // Pt[2]              /.      /.      /|
   Pt.push_back(Point(      0, halfEL,      0)); // Pt[3]             / .    8/ .     / |
   Pt.push_back(Point(      0,      0,      0)); // Pt[4]          12+-------+-------+5 |
   Pt.push_back(Point( halfEL,      0, halfEL)); // Pt[5]           /.10+.../..3+.../|..+2
   Pt.push_back(Point( halfEL, halfEL, halfEL)); // Pt[6]          / .     / .     / | /
   Pt.push_back(Point(      0, halfEL, halfEL)); // Pt[7]       13+-------+16-----+17|/
   Pt.push_back(Point(      0,      0, halfEL)); // Pt[8]         | 9+....|..+4...|..+1
   Pt.push_back(Point(-halfEL,      0,      0)); // Pt[9]         | .     | .     | /
   Pt.push_back(Point(-halfEL, halfEL,      0)); // Pt[10]        |.      |.      |/
   Pt.push_back(Point(-halfEL, halfEL, halfEL)); // Pt[11]      14+-------+-------+18
   Pt.push_back(Point(-halfEL,      0, halfEL)); // Pt[12]               15
   Pt.push_back(Point(-halfEL,-halfEL, halfEL)); // Pt[13]
   Pt.push_back(Point(-halfEL,-halfEL,      0)); // Pt[14]                      3
   Pt.push_back(Point(      0,-halfEL,      0)); // Pt[15]    Z<0     10+-------+-------+2
   Pt.push_back(Point(      0,-halfEL, halfEL)); // Pt[16]             /.      /.      /|
   Pt.push_back(Point( halfEL,-halfEL, halfEL)); // Pt[17]            / .    4/ .     / |
   Pt.push_back(Point( halfEL,-halfEL,      0)); // Pt[18]          9+-------+-------+1 |
   Pt.push_back(Point( halfEL,      0,-halfEL)); // Pt[19]          /.24+.../.21+.../|..+20
   Pt.push_back(Point( halfEL, halfEL,-halfEL)); // Pt[20]         / .     / .     / | /
   Pt.push_back(Point(      0, halfEL,-halfEL)); // Pt[21]      14+-------+15-----+18|/
   Pt.push_back(Point(      0,      0,-halfEL)); // Pt[22]        |23+....|..+22..|..+19
   Pt.push_back(Point(-halfEL,      0,-halfEL)); // Pt[23]        | .     | .     | /
   Pt.push_back(Point(-halfEL, halfEL,-halfEL)); // Pt[24]        |.      |.      |/
   Pt.push_back(Point(-halfEL,-halfEL,-halfEL)); // Pt[25]      25+-------+-------+27
   Pt.push_back(Point(      0,-halfEL,-halfEL)); // Pt[26]               26
   Pt.push_back(Point( halfEL,-halfEL,-halfEL)); // Pt[27]
// Rotate the points starting from their initial position, then translate them to their final position.
   Vect OC(Point(0,0,0),Center);
   rotNtrans(rots, OC, Pt);
   return Pt;
}

/*!
   Start mesh building process for revolution objects.
*/
void SubdivisionMesh::initRevMesh(const number_t nbsubdom, const real_t radius1, const real_t radius2,
                                  const std::vector<Point>& CharacPts, const std::vector<ShapeInfo>& vSI,
                                  real_t& R1, real_t& R2, Point& P1, Point& P2, std::vector<ShapeInfo>& cvSI,
                                  std::string& bottomPt, std::string& topPt, bool& iscone, std::string& shape,
                                  DefaultGeometry *& DG, SurfPlane *& SP, SurfCone *& SC,
                                  std::vector <PatchGeometry *>& EBS, int& nb_sl, real_t& slice_height) {
   cvSI = vSI;
   bottomPt=string("1"), topPt=string("2");
   // Since one of the radii may be close to 0, we identify what we call the "first" basis of the
   // object, defined by (P1,R1) with R1 the largest radius, because we need geometrically separated
   // points to start the mesh.
   if (radius1 < radius2) {
      R1 = radius2; P1 = CharacPts[1];
      R2 = radius1; P2 = CharacPts[0];
      cvSI[0] = vSI[1]; cvSI[1] = vSI[0];
      swap(bottomPt, topPt);
   }
   else {
      R1 = radius1; P1 = CharacPts[0];
      R2 = radius2; P2 = CharacPts[1];
   }
   iscone = R1 != R2;
   shape = iscone ? "cone" : "cylinder";
   DG = new DefaultGeometry();
   SP = new SurfPlane();
   SC = new SurfCone(P1,P2,R1,R2);
   EBS = endBoundaryShapes(cvSI,*SC);

//  Compute effective number of slices (nb_sl) and height of a slice (slice_height)
//  -> The number of slices is now computed to be a multiple of the number of requested
//     number of subdomains (nbsubdom) along the axis.
   nb_sl = 0;// = nbslices formerly, imposed by input argument nbslices
   slice_height = SC -> Height();
   if (nb_sl == 0) {// force computation of number of slices of elements, which is now the default behavior
      nb_sl = int((slice_height + 0.5*R1) / R1);
      if (nb_sl == 0) { nb_sl = 1; }
   }
   number_t nbdom = nbsubdom;
   if (nbdom == 0) { nbdom = 1; }
   nb_sl = std::max(number_t(1), nb_sl/nbdom) * nbdom;
   slice_height /= nb_sl;
}

   //--- Fig4TeX printing functions

/*!
 Prints, on stream ftex, fig4tex instructions to define vertices (action = 0) or
 to display vertices (action = 1) of any order in a tetrahedron mesh.
 The vertices belong to area num of kind TA, or to any area of kind TA if num=0.
 */
void SubdivisionMesh::printTeXInArea(ostream& ftex, const number_t action,
                                    const topologicalArea TA, const number_t num) const {
   vector<number_t> VonB = rk_verticesIn(TA,num);

   if (num > 0) {
      ftex << "% Vertices on " << TG_.kindOf(TA) << " " << num << endl; }
   else {
      ftex << "% " << TG_.kindOf(TA) << " vertices" << endl; }

   switch (action) {
   default:
   case 0: // defines points
      printTeXfigpt(ftex,VonB);
      break;
   case 1: // display points
      printTeXfigwrite(ftex,VonB);
      break;
   }
}
/*!
 Prints, on stream ftex, fig4tex instructions to define vertices of a tetrahedron
 mesh stored in vector V.
 */
void SubdivisionMesh::printTeXfigpt(ostream& ftex, const vector<number_t>& V) const {
   vector<number_t>::const_iterator itV;
   for (itV=V.begin(); itV != V.end(); itV++) {
      ftex << "\\figpt " << listV_[*itV].number() << ":";
      listV_[*itV].printTeX(ftex);
      ftex << endl;
   }
}
/*!
 Prints, on stream ftex, fig4tex instructions to display vertices of a tetrahedron
 mesh stored in vector V.
 */
void SubdivisionMesh::printTeXfigwrite(ostream& ftex, const vector<number_t>& V) const {
   vector<number_t>::const_iterator itV = V.begin();
   ftex << "\\def\\dist{4pt}\\figwriten " << listV_[*itV].number();
   for (itV++; itV != V.end(); itV++) {
      ftex << "," << listV_[*itV].number();
   }
   ftex << ":(\\dist)" << endl;
}
/*!
 Prints, on stream ftex, Fig4TeX instructions to draw faces of a mesh
 belonging to area num of kind TA.
 */
void SubdivisionMesh::printTeXFacesInArea(ostream& ftex, const topologicalArea TA,
                                         const number_t num) const {
   vector< vector<number_t> > FonB = rk_facesIn(TA,num);
   vector< vector<number_t> >::iterator itFonB;
   vector<number_t>::iterator itF;

   ftex << "% Faces on " << TG_.kindOf(TA) << " " << num << endl;
   ftex << "\\def\\FaceColor{"<< colorOf(TA,num) <<"}" << endl;

   for (itFonB=FonB.begin(); itFonB != FonB.end(); itFonB++) {
      ftex << "\\drawFace";
      for (itF=itFonB->begin(); itF != itFonB->end(); itF++) {
         ftex << "{" << listV_[*itF].number() << "}";
      }
      ftex << endl;
   }
}

   //--- Low level utilitary functions

/*!
 Utilitary function to replace the rank in the internal list by the vertex number
 of each vertex stored in the vector V.
 */
void SubdivisionMesh::rankToNum(vector<number_t>& V) const {
   vector<number_t>::iterator itV;

   for (itV=V.begin(); itV != V.end(); itV++) {
      *itV = listV_[*itV].number(); // should be equivalent to (*itV)+minVertexNum_
   }
}
/*!
 Utilitary function to replace the rank in the internal list by the vertex number
 of each vertex stored in the pair V.
 */
void SubdivisionMesh::rankToNum(pair_nn& V) const {
   V.first  = listV_[V.first].number();  // should be equivalent to (V.first)+minVertexNum_
   V.second = listV_[V.second].number(); // should be equivalent to (V.second)+minVertexNum_
}

/*!
 Returns the list of vertices of any order in a tetrahedron mesh, which are
 not on the boundaries. Thus, they are inside the domain.
 Vertices are defined by their rank in the general internal list of vertices.
 */
vector<number_t> SubdivisionMesh::rk_verticesInside() const {
   vector<number_t> VI;
   refnum_t BMask = TG_.maskOf(boundaryArea);
   number_t k;

   vector<Vertex>::const_iterator itV;
   for (itV=listV_.begin(), k=0; itV != listV_.end(); itV++,k++) {
      refnum_t localcod = itV->locCode() & BMask;
      // If the vertex belongs to none of the boundaries, store its rank.
      if (! localcod) {VI.push_back(k);}
   }
   return VI;
}
/*!
 Returns the list of vertices of any order in a tetrahedron mesh, belonging to
 the area num of kind TA or belonging to any area of kind TA if num=0.
 Vertices are defined by their rank in the general internal list of vertices.
 */
vector<number_t> SubdivisionMesh::rk_verticesIn(const topologicalArea TA, const number_t num) const {
   vector<number_t> VonB;
   number_t k;
   refnum_t sig = lCodeOf(TA,num);

   vector<Vertex>::const_iterator itV;
   for (itV=listV_.begin(), k=0; itV != listV_.end(); itV++,k++) {
      refnum_t localcod = itV->locCode();
      // If the vertex belongs to the area, store its rank.
      if (localcod & sig) {VonB.push_back(k);}
   }
   return VonB;
}

   //--- For high order vertices generation

/*!
 Comparison between two vector<number_t> in lexicographical order according to the
 first two components.
 */
bool SubdivisionMesh::cmpvect(const vector<number_t>& U, const vector<number_t>& V){
   if (U[0] < V[0]) {return true;}
   else if (U[0] == V[0]) {return U[1] < V[1];}
   else return false;
}
//-------------------------------------------------------------------------------
//  Private member functions
//-------------------------------------------------------------------------------
/*!
 Default function that computes the geometrical point associated to a new vertex.
 The new point is the barycenter of points in VP with coefficients coef,
 whatever the localization code localcod is.
 */
Point SubdivisionMesh::newVertexPtDef(const refnum_t localcod, const real_t *coef,
                                      const vector<Point>& VP) const{
   return barycenter(coef,VP);
}
/*!
 Computational routine of the geometrical point associated to a new vertex,
 taking into account the definition of the areas of the domain by the
 mean of the localization code localcod of the set of points in VP.
 Let P be the barycenter of the points in VP with coefficients coef.
 The new point is:
  - the projection of P on a curved area, if points in VP all belong to this area,
  - the point P otherwise.
 If the points in VP belong to several areas, the first curved area checked is
 considered, which is consistent since the areas are assumed to have continuous
 connections.
 */
Point SubdivisionMesh::newVertexPtGen(const refnum_t localcod, const real_t *coef,
                                      const vector<Point>& VP) const{
// Check whether points in VP belong to a curved patch whose list has been prepared in TG_
   for (number_t i=0; i < TG_.numberOfCurvedPatches(); i++) {
      pair<PatchGeometry*,refnum_t> CP = TG_.shapeNlcode(i);
      if (localcod  &  CP.second) {
         // points in VP all belong to this curved patch
         return (CP.first) -> projOnBound(coef,VP);
      }
   }

// Default behavior when points in VP do not belong to any curved patch
   return barycenter(coef,VP);
}


int DECHOL(const vector<vector<real_t> >& A, int n, vector<vector<real_t> >& L, real_t eps) {
      int i, j, k;
      real_t s;
/*
       ----------------------------------------------------------------
       Decomposition de Cholesky de la matrice A sous la forme A = LLt.
       Le calcul de la matrice L est fait colonne par colonne.
         Arguments d'entree  : A,n,eps
         Arguments de sortie : L
       A      : matrice donnee d'ordre n symetrique
       n      : ordre de la matrice
       L      : matrice triangulaire inferieure resultat de la
                decomposition de A (si A n'est pas reutilisee, A et L
                peuvent correspondre au meme emplacement memoire)
       eps    : reel servant a determiner la positivite de la matrice A
       La fonction renvoie un code de retour valant 0 si le calcul a abouti,
       1 si la matrice A n'est pas definie positive.

       --> seuls sont utilises les elements des matrices
           triangulaires inferieures
       --> utiliser ensuite DRCHOL pour resoudre un systeme lineaire
       ----------------------------------------------------------------
*/
      for (i=0; i<n; i++) {
    /* Pre-calcul de la ieme colonne de L */
         for (j=i; j<n; j++) {
            s = A[j][i];
            for (k=0; k<i; k++) { s -= L[i][k]*L[j][k]; }
            L[j][i] = s;
         }
    /* Calcul du terme diagonal et test de positivite */
         if (L[i][i] < 0.) return 1;
         L[i][i] = std::sqrt(L[i][i]);
         if (L[i][i] < eps) return 1;
    /* Division de la colonne par l'element diagonal */
         s = L[i][i];
         for (j=i+1; j<n; j++) { L[j][i] /= s; }
      }
      return 0;
}
// *********************************************************************** //
void DRCHOL(const vector<vector<real_t> >&L, int n, vector<vector<real_t> >&B, int m) {
      int i, j, l;
      real_t s;
/*
      -----------------------------------------------------------------
      Resolution des m systemes lineaires Mxj=Bj, j=1,m.
      L est le resultat de la decomposition de Cholesky (cf. module
      DECHOL) d'une matrice M = LLt, L triangulaire inferieure.
      La resolution se fait par descente-remontee.
         Arguments d'entree  : L,n,m
         Arguments de sortie : B
       L      : matrice triangulaire inferieure de Cholesky
       n      : ordre de la matrice
       B      : matrice (n,m) contenant les seconds membres et les solutions
       m      : nombre de systemes a resoudre
      --> En sortie, les solutions sont memorises a l'emplacement des
          seconds membres.
      -----------------------------------------------------------------
*/
/*        Pour chaque second membre : */
      for (l=0; l<m; l++) {
/*        Descente      */
         for (i=0; i<n; i++) {
            s = B[i][l];
            for (j=0; j<i; j++) { s -= L[i][j]*B[j][l]; }
            B[i][l] = s / L[i][i];
          }
/*        Remontee      */
         for (i=n-1; i>=0; i--) {
            s = B[i][l];
            for (j=i+1; j<n; j++) { s -= L[j][i]*B[j][l]; }
            B[i][l] = s / L[i][i];
         }
      }
}
} // end of namespace subdivision
} // end of namespace xlifepp