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

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

#include "GeomFigureMesh.hpp"
#include "Hexahedron.hpp"
#include "Tetrahedron.hpp"
#include "Quadrangle.hpp"
#include "Triangle.hpp"
#include "TeXPolygon.hpp"

#include <algorithm>
#include<vector>
using namespace std;

namespace xlifepp {
namespace subdivision {

//-------------------------------------------------------------------------------
//  Public member functions
//-------------------------------------------------------------------------------

/*!
 Returns the list of all the vertices of the element num, with
 minElementNum_ <= num <= numberOfElements+minElementNum_-1.
 Each vertex is given by its number, which is >= minVertexNum_.
 */
template<class T_>
vector<number_t> GeomFigureMesh<T_>::element(const number_t num) const {
   vector<number_t> V = listT_.at(num-minElementNum_).rankOfVertices();
   rankToNum(V);
   return V;
}

/*!
  returns the localization code of the element num, with num such that
  minElementNum_ <= num <= numberOfElements+minElementNum_-1
  The code is computed from those of its main vertices. The function
  returns the localization code of the domain the element belongs to.
*/
template<class T_>
refnum_t GeomFigureMesh<T_>::locCode(const number_t num) const {
   const T_& elem(listT_.at(num-minElementNum_));
   vector<number_t> V = elem.rankOfVertices();
   refnum_t localcod = listV_[V[0]].locCode();
   for (number_t j=1; j<elem.numberOfO1Vertices(); j++) {
      localcod &= listV_[V[j]].locCode();
   }
   return localcod;
}

/*!
  returns a normal vector to the face number i >= 1 of element num,
  with num such that minElementNum_ <= num <= numberOfElements+minElementNum_-1
  The vector is not normalized and is oriented towards the exterior of the element.
 */
template<class T_>
vector<real_t> GeomFigureMesh<T_>::faceExtNormalVec(const number_t num, const number_t i) const {
   return listT_.at(num-minElementNum_).extNormVec(i, listV_);
}

/*!
  returns the orientation of the face number i >= 1 of element num, with num such that
  minElementNum_ <= num <= numberOfElements+minElementNum_-1. The function:
  - returns 1 if the first three vertices (i,j,k) of the face define a vector ij x ik
    oriented towards the exterior of the element,
  - returns -1 otherwise.
*/
template<class T_>
int GeomFigureMesh<T_>::faceOrientation(const number_t num, const number_t i) const {
   return listT_.at(num-minElementNum_).faceOrientation(i);
}

/*!
 Returns the list of elements of the mesh belonging to an area of kind TA
 or belonging to any area of kind TA if num=0. The elements are given by their
 number in the whole set of elements of the mesh.
 */
template<class T_>
vector<number_t> GeomFigureMesh<T_>::elementsIn(const topologicalArea TA, const number_t num) const {
   vector<number_t> V;
   refnum_t sig = lCodeOf(TA,num);

   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      refnum_t localcod = listV_[itT->rankOfVertex(1)].locCode();
      for (number_t j=2; j<=nb_main_vertices_by_element_; j++) {
         localcod &= listV_[itT->rankOfVertex(j)].locCode();
      }
      // If the element belongs to the area, store its number.
      if (localcod & sig) {V.push_back(itT->number());}
   }
   return V;
}

/*!
 Returns the number of edges belonging to an area num of kind TA
 or belonging to any area of this kind if num=0.
 */
template<class T_>
number_t GeomFigureMesh<T_>::numberOfEdgesIn(const topologicalArea TA, const number_t num) const {
   return rk_edgesIn(TA,num).size();
}
/*!
 Returns the list of edges of the elements of the mesh belonging to an area num
 of kind TA or belonging to any area of this kind if num=0.
 Each edge is defined by a pair of the form (nelem[i], nedge[i]) where
   . nelem[i] is the number of the element (>=minElementNum_) the edge belongs to.
   . nedge[i] is the local number of the edge in the element (from 1 to nb_edges_by_element_)
 These numbers are stored in two separate vectors nelem and nedge, with the same
 length, grouped as a pair returned by the function.
 -> The core of the algorithm is the same in rk_edgesIn.
 Nota: If TA is an interface, the result consists in fact in two lists, since each
       segment of the interface belongs to 2 elements. The two lists are intermixed,
       except when there are only two subdomains, in which case, due to the construction
       algorithm, the first half of the elements belong to one subdomain and the second
       half of the elements belong to the other subdomain.
       If needed, the two lists should be separated using geometric computations.
 */
template<class T_>
pair< vector<number_t>,vector<number_t> >
GeomFigureMesh<T_>::edgeElementsIn(const topologicalArea TA, const number_t num) const {
   vector<number_t> numEdge;
   vector<number_t> numElem;
   refnum_t sig = lCodeOf(TA,num);

   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      for (number_t i=1; i<=nb_edges_by_element_; i++) {
         pair_nn prV = itT->rkOfO1VeOnEdge(i);
         refnum_t localcod = listV_[prV.first].locCode() & listV_[prV.second].locCode();
         // If the edge belongs to the area, store it.
         if (localcod & sig) { numEdge.push_back(i); numElem.push_back(itT->number()); }
      }
   }
   return make_pair(numElem, numEdge);
}
/*!
 Returns the list of edges of a triangle mesh belonging to an area num of kind TA
 or belonging to any area of this kind if num=0.
 Edges are defined by the number of their vertices (starting from minVertexNum_).
 */
template<class T_>
vector< pair_nn > GeomFigureMesh<T_>::edgesIn(const topologicalArea TA, const number_t num) const {
   vector< pair_nn > EonB = rk_edgesIn(TA,num);
   vector< pair_nn >::iterator itEonB;

   for (itEonB=EonB.begin(); itEonB != EonB.end(); itEonB++) { rankToNum(*itEonB); }
   return EonB;
}

//  functions for 3D geometries
/*!
 Returns the number of faces belonging to an area num of kind TA
 or belonging to any area of this kind if num=0.
 */
template<class T_>
number_t GeomFigureMesh<T_>::numberOfFacesIn(const topologicalArea TA, const number_t num) const {
   return rk_facesIn(TA,num).size();
}
/*!
 Returns the number of the internal faces of the mesh
 */
template<class T_>
number_t GeomFigureMesh<T_>::numberOfFacesInside() const {
// lazy algorithm: it would be more efficient not to buid the list in memory,
//                 but just count the internal faces as they appear in the loop
//                  (see function facesInside).
   return facesInside().size();
}

/*!
 Returns the list of faces of the elements of the mesh belonging to an area num
 of kind TA or belonging to any area of this kind if num=0.
 Each face is defined by a pair of the form (nelem[i], nface[i]) where
   . nelem[i] is the number of the element (>=minElementNum_) the face belongs to,
   . nface[i] is the local number of the face in the element (from 1 to nb_faces_by_element_).
 These numbers are stored in two separate vectors nelem and nface, with the same
 length, grouped as a pair returned by the function.
 -> The core of the algorithm is the same in rk_facesIn.
 Nota: If TA is an interface, the result consists in fact in two lists, since each 2D
       element of the interface belongs to two 3D elements. The two lists are intermixed,
       except when there are only two subdomains, in which case, due to the construction
       algorithm, the first half of the elements belong to one subdomain and the second
       half of the elements belong to the other subdomain.
       If needed, the two lists should be separated using geometric computations.
 */
template<class T_>
pair< vector<number_t>,vector<number_t> >
GeomFigureMesh<T_>::faceElementsIn(const topologicalArea TA, const number_t num) const {
   vector<number_t> numFace;
   vector<number_t> numElem;
   refnum_t sig = lCodeOf(TA,num);

   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      for (number_t i=1; i<=nb_faces_by_element_; i++) {
         vector<number_t> vrV = itT->rkOfO1VeOnFace(i);
         vector<number_t>::const_iterator itvrV(vrV.begin());
         refnum_t localcod = listV_[*itvrV].locCode();
         while (++itvrV < vrV.end()) { localcod &= listV_[*itvrV].locCode(); }
         // If the face belongs to the area, store it.
         if (localcod & sig) { numFace.push_back(i); numElem.push_back(itT->number()); }
      }
   }
   return make_pair(numElem, numFace);
}

/*!
 Returns the list of faces of a mesh belonging to an area num of kind TA
 or belonging to any area of this kind if num=0.
 Faces are defined by the number of their vertices (starting from minVertexNum_).
 */
template<class T_>
vector< vector<number_t> > GeomFigureMesh<T_>::facesIn(const topologicalArea TA, const number_t num) const {
   vector< vector<number_t> > FonB = rk_facesIn(TA,num);
   vector< vector<number_t> >::iterator itFonB;

   for (itFonB=FonB.begin(); itFonB != FonB.end(); itFonB++) { rankToNum(*itFonB); }
   return FonB;
}

/*!
 Returns the list of faces of a mesh which are inside the boundaries.
 Each face is defined by a pair of the form ( vertices, VnTnF ) where
   . vertices is a set of integers, corresponding to the numbers of the vertices
     defining the face,
   . VnTnF is a vector of two pairs of the form (nelem, nface) corresponding to the
     two elements sharing the face, where
     . nelem is the number (>=1) of one of the two elements the face belongs to,
     . nface is the local number of the face in this element (from 1 to nb_faces_by_element_).
 */
template<class T_>
map_set_vec GeomFigureMesh<T_>::facesInside() const {
   map_set_vec FI;
   pair< map_set_vec::iterator, bool > retval;
   refnum_t BMask = TG_.maskOf(boundaryArea);

   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      for (number_t i=1; i<=nb_faces_by_element_; i++) {
         vector<number_t> vrV = itT->rkOfO1VeOnFace(i);
         vector<number_t>::const_iterator itvrV(vrV.begin());
         refnum_t localcod = BMask & listV_[*itvrV].locCode();
         while (++itvrV < vrV.end()) { localcod &= listV_[*itvrV].locCode(); }
     // If the vertices of the face do not belong all to the same boundary patch,
     // the face is internal and is stored.
     // Since each such face belongs to two elements, in order to ensure the
     // face is stored only once, the vector vrV is converted into a set (ordered) used
     // as the key value of the map, whose occurrence is thus unique.
     if (! localcod) {
        rankToNum(vrV);
        set_n Face(vrV.begin(),vrV.end());
        pair_nn nTnF(make_pair(itT->number(),i));
        vector< pair_nn > VnTnF; VnTnF.push_back(nTnF);
        retval = FI.insert(make_pair(Face, VnTnF));
        if (! retval.second) {// add second element sharing this face
           (retval.first->second).push_back(nTnF);
           }
        }
      }
   }
   return FI;
}

//-------------------------------------------------------------------------------
//  I/O utilities
//-------------------------------------------------------------------------------
/*!
 Prints, on stream os, the definition of a mesh with vertices coordinates.
 */
template<class T_>
void GeomFigureMesh<T_>::printall(ostream& os) const {
   printInfo(os);

   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      os << "Element " << itT->number() << endl;
      for (number_t j=1; j<=nb_main_vertices_by_element_; j++) {
         listV_[itT->rankOfVertex(j)].print(os,TG_);
      }
   }
   os << endl;

   os << endl << "List of vertices :" << endl;
   vector<Vertex>::const_iterator itV;
   for (itV=listV_.begin(); itV != listV_.end(); itV++) {itV->print(os,TG_);}
   os << endl;
}

/*!
 Prints, on stream os, the definition of a mesh.
 */
template<class T_>
void GeomFigureMesh<T_>::print(ostream& os) const {
   printInfo(os);

   os << endl << "List of elements (vertices are given by their number):" << endl;
   for (number_t nel=1; nel<=numberOfElements(); nel++) {
      os << "Element " << nel << " : ";
      vector<number_t> nV(element(nel));
      for (size_t i=0; i<nV.size(); i++) {os << nV[i] << " ";}
      os << endl;
   }

   os << endl << "List of elements (vertices are given by their rank in the vertex list):" << endl;
   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      itT->print(os);
   }

   os << endl << "List of vertices :" << endl;
   vector<Vertex>::const_iterator itV;
   number_t k;
   for (itV=listV_.begin(), k=0; itV != listV_.end(); itV++,k++) {
      os << "Rank " << setw(4) << k << " ";
      itV->print(os,TG_);}
   os << endl;
}

/*!
 Prints, on stream ftex, Fig4TeX instructions to draw a mesh.
 The observation direction is defined by the longitude Psi and the latitude Theta,
 given in degrees. nbviews figures are drawn, all with the same latitude, the
 longitude varying by a 360/nbviews increment.
*/
template<class T_>
void GeomFigureMesh<T_>::printTeX(ostream& ftex, const float Psi, const float Theta,
                               const number_t nbviews, const string& DimProj,
                               const bool withInterface, const bool withElems) const {
   ftex << "\\let\\showfigOne\\centerline" << endl;
   ftex << "\\def\\showfigTwo#1#2{\\centerline{#1}\\nobreak\\medskip\\centerline{#2}}" << endl;
   ftex << "\\input fig4tex.tex" << endl;
   printTeXHeader(ftex);
   ftex << "%" << endl;
   ftex << "% 1. Definition of characteristic points" << endl;
   ftex << "\\figinit{"<< DimProj << "}" << endl;

   bool drawInterface = withInterface && numberOfInterfaces() > 0;
   printTeXPoints(ftex,drawInterface);
   ftex << "%" << endl;

   float psi = Psi, deltapsi = 360./nbviews;
   for (number_t numv=1; numv <= nbviews; numv++) {
      ostringstream ss;
      ss << "Subdiv. level "<< subdiv_level_
            << ", long. "<< psi << "$^\\circ$, lat. " << Theta << "$^\\circ$";
      if (drawInterface) {
         createFileNBox(ftex,boundaryArea,psi,Theta,"A",string("Domain"));
         createFileNBox(ftex,interfaceArea,psi,Theta,"B",string("Interfaces"));
         ftex << "\\showfigTwo{\\box\\figBoxA\\hfil\\quad\\box\\figBoxB}{"
              << ss.str() << "}" << endl;
         }
      else {
         createFileNBox(ftex,boundaryArea,psi,Theta,"A",ss.str());
         ftex << "\\showfigOne{\\box\\figBoxA}" << endl;
         }
      ftex << "%-------------------------------- End of figure --------------------------------" << endl;
      ftex << "\\bigskip\\vfill" << endl;
      psi += deltapsi;
   }

// Draw the elements if requested (essentially usefull in 3D and to check high degree vertices)
   if (withElems) {
      ftex << "\\bigskip\\vfill\\eject" << endl;
      ftex << "% Draw all the elements of the mesh" << endl;
      printTeXInArea(ftex,0,subdomainArea);
      ftex << "\\figdrawbegin{}"<< endl;
      vector<number_t>::iterator itN;
      number_t ind, indx(numberOfMainVerticesByElement());
      // Draw all the elements of the mesh
      number_t nbElts = numberOfElements();
      for (number_t k=1; k<=nbElts; k++) {
        vector<number_t> N = element(k);
        ftex << "\\drawElem";
        for (itN=N.begin(), ind=0; ind<indx; ind++,itN++) {
          ftex << "{" << *itN << "}";
          }
        ftex << endl;
      }
      ftex << "\\figdrawend" << endl;
      ftex << "\\figvisu{\\figBoxA}{" << nbElts << " elements of order " << order_ << "}{" << endl;
      ftex << "% Write all the vertices as a whole" << endl;
      ftex << "%\\figshowpts[1," << numberOfVertices() << "]" << endl;
      ftex << "% Write all the vertices, element by element, including high order vertices if any" << endl;
      ftex << "\\figset write(ptname={\\bf{#1}})" << endl;
      // Write all the vertices, element by element, including high order vertices if any
      for (number_t k=1; k<=nbElts; k++) {
        vector<number_t> N = element(k);
        itN=N.begin();
        ftex << "\\figwritec[" << *itN++ ;
        for (; itN<N.end(); itN++) {
          ftex << "," << *itN ;
          }
        ftex << "]{}" << endl;
      }
      ftex << "}" << endl << "\\centerline{\\box\\figBoxA}" << endl;
   }
   ftex << "\\bye" << endl;
}

//-------------------------------------------------------------------------------
//  Protected member functions
//-------------------------------------------------------------------------------
// defined here because of instantiation, but never used in practice
template<class T_>
void GeomFigureMesh<T_>::algoSubdiv(const T_& T, number_t& ElementNum, number_t& VertexNum,
                                    vector<T_>& listT, map_pair_num& SeenEdges) {
   error("nofunc","GeomFigureMesh" ,"algoSubdiv");
}

   //! build the mesh, compute high order vertices and check consistency
template<class T_>
void GeomFigureMesh<T_>::buildNcheck(number_t& VertexNum) {
   if (order_ < 1) {order_ = 1;}
//   Subdivisions
   buildMesh(VertexNum);
//   Store the number of vertices of order 1
   nb_main_vertices_ = listV_.size();

//   End with computing high order vertices, if required, and some ckecks
   try {
      // Take order into account: compute additionnal high order vertices
      if (order_ > 1) {
         createHOV(VertexNum); // throw (HOV_error)
      }
      // Check number of vertices for each element and the number assigned to each of them
      checkElements(); // throw (number_t *)
      // Check the number assigned to each vertex
      vector<Vertex>::const_iterator itV;
      number_t k;
      for (itV=listV_.begin(), k=minVertexNum_; itV != listV_.end(); itV++, k++) {
         if (itV->number() != k) {
            throw std::make_pair(itV->number(), k);
         }
      }
   }
   catch (HOV_error pb) {
      cout << "*******************************************************" << endl;
      cout << "From subdivision mesh constructor" << endl;
      cout << " Error during high order vertices generation while processing" << endl;
      cout << " element " << pb.elemNum << " and vertex " << pb.vertexNum << endl;
      cout << pb.mes << endl;
      cout << " Rank of vertex 1 = " << pb.rV1 << endl;
      cout << " Rank of vertex 2 = " << pb.rV2 << endl;
      cout << " Rank of vertex 3 = " << pb.rV3 << endl;
      cout << "*******************************************************" << endl;
   }
   catch (number_t *values) {
      cout << "*******************************************************" << endl;
      cout << "From subdivision mesh constructor" << endl;
      cout << " Error detected during post construction check:" << endl;
      cout << "Computed element number : " << values[0] << endl;
      cout << "Expected element number : " << values[1] << endl;
      cout << values[2] << " vertices computed on this element " << endl;
      cout << values[3] << " vertices expected." << endl;
      cout << "*******************************************************" << endl;
   }
   catch (pair_nn values) {
      cout << "*******************************************************" << endl;
      cout << "From subdivision mesh constructor" << endl;
      cout << " Error detected during post construction check:" << endl;
      cout << "Computed vertex number : " << values.first << endl;
      cout << "Expected vertex number : " << values.second << endl;
      cout << "*******************************************************" << endl;
   }
}

/*!
 Build the mesh by successive subdivisions of the initial mesh.
 Input arguments :
   VertexNum   : number of the last vertex created in the initial mesh
 Output arguments :
   VertexNum   : number of the last vertex created in the mesh
  */
template<class T_>
void GeomFigureMesh<T_>::buildMesh(number_t& VertexNum){
//  At each subdivision, the elements are numbered starting from minElementNum_ (via
//  the statement ElementNum=minElementNum_-1), since theses are all new elements.
//  On the contrary, the vertices are all kept from the beginning and the vertex counter
//  VertexNum must not be reinitialized.
//  Nota: By construction, the elements stored in listT_ are grouped by subdomains.
   for (number_t i=0; i<subdiv_level_; i++) {
      number_t ElementNum = minElementNum_ - 1;
      vector<T_> tmplistT;
      // listT_.size() is the number of elements built during previous subdivision
      tmplistT.reserve(subdivisionFactor_*listT_.size());
      map_pair_num SeenEdges;
      typename vector<T_>::iterator itT;
      for (itT=listT_.begin(); itT != listT_.end(); itT++) {
         algoSubdiv(*itT,ElementNum,VertexNum,tmplistT,SeenEdges);
      }
      listT_ = tmplistT;
   }
   initDefaultUserAttribute();
}
/*!
 Check number of vertices for each element and the number assigned to each of them
 */
template<class T_>
void GeomFigureMesh<T_>::checkElements() {
   number_t k;
   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(), k=minElementNum_; itT != listT_.end(); itT++, k++) {
      if (itT->numberOfVertices() != total_nb_vertices_by_element_ || itT->number() != k) {
         number_t values[] = {itT->number(), k, itT->numberOfVertices(), total_nb_vertices_by_element_};
         throw values;
      }
   }
}

/*!
 Create high order vertices.
 Input arguments :
 \param VertexNum : number of the last vertex created so far in the mesh
 Output arguments :
 \param VertexNum : number of the last vertex created in the mesh by this procedure
 */
template<class T_>
void GeomFigureMesh<T_>::createHOV(number_t& VertexNum){
   map_set_pair SeenEdges;
   map_set_pair_in SeenFaces;
   typename vector<T_>::iterator itT;

   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
// Vertices inside each edge
      for (number_t numEdge=1; numEdge<=nb_edges_by_element_; numEdge++) {
         createHOeV(*itT,order_,VertexNum,numEdge,SeenEdges);
      }
// Vertices inside each face
      for (number_t numFace=1; numFace<=nb_faces_by_element_; numFace++) {
         createHOfV(*itT,order_,VertexNum,numFace,SeenFaces);
      }
// Internal vertices
      createHOiV(*itT,order_,VertexNum);
   }
}
/*!
 Create or retrieve high order vertices on the edge numEdge of element Elem.
 Their numbers are consecutive on the edge.
 The creation process of the new vertices and the numbering defined in
 T_::numberingOfVertices are coherent.
 */
template<class T_>
void GeomFigureMesh<T_>::createHOeV(T_& Elem, const number_t order, number_t& VertexNum,
                                    const number_t numEdge, map_set_pair& SeenEdges){
   pair_nn prV = Elem.rkOfO1VeOnEdge(numEdge);
   // Map key is a set in order to create a unique entry in the map, regardless of
   // the way the edge is defined.
   number_t rV1 = prV.first;
   number_t rV2 = prV.second;
   set_n Edge; Edge.insert(rV1); Edge.insert(rV2);
   map_set_pair::iterator itSE;

   if ((itSE=SeenEdges.find(Edge)) == SeenEdges.end()) {
      // This edge has not yet been seen: compute the corresponding vertices
      // Store the edge along with the corresponding couple (rank of the
      // first new vertex created in the edge, rank of the first vertex
      // defining the edge)
      pair_nn NV (VertexNum, rV1);
      SeenEdges.insert(make_pair(Edge,NV));

      // Compute the new vertices along this edge
      Vertex V1=listV_[rV1];
      Vertex V2=listV_[rV2];
      refnum_t localcod = V1.locCode() & V2.locCode();
      // We compute the position of each point as a barycenter of the two end-points of the edge.
      // The numbering is decreasing from V1 towards V2.
      //           V1                                                                           V2
      //           |----------|----------|----------|----------|----------|----------|----------|
      // coef[0] =            1          2             . . .             k-2        k-1
      // coef[1] =           k-1        k-2            . . .              2          1
      vector<Point> VP(2);
      VP[0] = V1.geomPt();
      VP[1] = V2.geomPt();
      real_t coef[2];
      Point P;
      for (number_t i1=1; i1<order; i1++) {
         coef[0] = i1;
         coef[1] = order-i1;
         P = (this->*newVertexPt_)(localcod,coef,VP);
         Elem.vertices_.push_back(VertexNum);// add rank of the new vertex
         listV_.push_back(Vertex(++VertexNum,localcod,P));// store the new vertex
      }
   }
   else {
// Edge found: get back the corresponding vertices whose numbers are consecutive
// along this edge, starting from num.
      number_t num = (itSE->second).first;
      // Find the edge orientation
      if ( (itSE->second).second == rV1){
         for (number_t i1=1; i1<order; i1++) {
            Elem.vertices_.push_back(num++);// add rank of the vertex
         }
      }
      else {
         num += order-2;
         for (number_t i1=1; i1<order; i1++) {
            Elem.vertices_.push_back(num--);// add rank of the vertex
         }
      }
   }
}
// defined here because of instantiation, but never used in practice
template<class T_>
void GeomFigureMesh<T_>::printTeXHeader(ostream& ftex) const {
   error("nofunc","GeomFigureMesh","printTeXHeader");
}
/*!
 Prints, on stream ftex, Fig4TeX instructions to define points to be used for
 the display of the mesh. These points are located on the boundaries of the domain
 or on the interfaces.
 */
template<class T_>
void GeomFigureMesh<T_>::printTeXPoints(ostream& ftex, const bool withInterface) const {
   printTeXInArea(ftex,0,boundaryArea);
   if (withInterface) { printTeXInArea(ftex,0,interfaceArea); }
}

/*!
 Prints, on stream ftex, Fig4TeX instructions to draw faces of a mesh belonging
 to any area of kind TA. The faces are ordered from the farthest to the nearest
 along the observation direction defined by the longitude psi and the latitude
 theta, given in degrees.
 */
template<class T_>
void GeomFigureMesh<T_>::printTeXSortedAreaFoE(ostream& ftex, const topologicalArea TA,
                                               const float psi, const float theta) const {
   vector<TeXPolygon> BF;
   number_t number_of_area = TG_.numberOf(TA);

//  Collect all the boundary faces and keep their associated boundary number
   for (number_t numArea=1; numArea<=number_of_area; numArea++) {
      vector< vector<number_t> > FonB = rk_facesIn(TA,numArea);
      vector< vector<number_t> >::iterator itFonB;
      for (itFonB=FonB.begin(); itFonB != FonB.end(); itFonB++) {
         BF.push_back(TeXPolygon(*itFonB,numArea,listV_));
      }
      // Define boundary colors
      ftex << "\\def\\Color" << string(1,'@'+numArea)  // ColorA, ColorB...
           << "{" << colorOf(TA,numArea) << "}% " << TG_.nameOf(TA,numArea) << endl;
   }
//  Sort the faces from the farthest to the nearest
//  1. Compute the vector defining the observation direction
   real_t psird = psi*pi_/180.;
   real_t thetard = theta*pi_/180.;
   real_t ct = std::cos(thetard);
//   TeXPolygon::OD = Vect(ct*cos(psird), ct*sin(psird), sin(thetard));
   TeXPolygon::initObsDir(Vect(ct*std::cos(psird), ct*std::sin(psird), std::sin(thetard)));
//  2. Apply sort algorithm
   sort(BF.begin(),BF.end());
//  3. Print the sorted faces
   if (number_of_area > 0) {
      ftex << "% " << BF.size() << " faces on " << TG_.kindOf(TA) << " " << TG_.nameOf(TA,1);
   }
   for (number_t numArea=2; numArea<=number_of_area; numArea++) {
      ftex << ", " << TG_.nameOf(TA,numArea);
   }
   ftex << endl;
   vector<TeXPolygon>::iterator itBF;
   vector<number_t>::const_iterator itF;
   for (itBF=BF.begin(); itBF != BF.end(); itBF++) {
      ftex << "\\drawFace";
      for (itF=itBF->Vrank().begin(); itF != itBF->Vrank().end(); itF++) {
         ftex << "{" << listV_[*itF].number() << "}";
      }
      ftex << "{\\Color" << string(1,'@'+itBF->attrib()) << "}";
      ftex << endl;
   }
}

/*!
 Returns the list of edges of a mesh belonging to area num of kind TA
 or belonging to any area of this kind if num=0.
 Edges are defined by the rank of their vertices in the internal list of vertices.
 */
template<class T_>
vector< pair_nn > GeomFigureMesh<T_>::rk_edgesIn(const topologicalArea TA, const number_t num) const {
   vector< pair_nn > EonB;
   refnum_t sig = lCodeOf(TA,num);

   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      for (number_t i=1; i<=nb_edges_by_element_; i++) {
         pair_nn prV = itT->rkOfO1VeOnEdge(i);
         refnum_t localcod = listV_[prV.first].locCode() & listV_[prV.second].locCode();
         // If the edge belongs to the area, store it.
         if (localcod & sig) {EonB.push_back(prV);}
      }
   }
   return EonB;
}

/*!
 Returns the list of faces of a mesh belonging to area num of kind TA
 or belonging to any area of this kind if num=0.
 Faces are defined by the rank of their vertices in the internal list of vertices.
 */
template<class T_>
vector< vector<number_t> > GeomFigureMesh<T_>::rk_facesIn(const topologicalArea TA, const number_t num) const {
   vector< vector<number_t> > FonB;
   refnum_t sig = lCodeOf(TA,num);

   typename vector<T_>::const_iterator itT;
   for (itT=listT_.begin(); itT != listT_.end(); itT++) {
      for (number_t i=1; i<=nb_faces_by_element_; i++) {
         vector<number_t> vrV = itT->rkOfO1VeOnFace(i);
         vector<number_t>::const_iterator itvrV(vrV.begin());
         refnum_t localcod = listV_[*itvrV].locCode();
         while (++itvrV < vrV.end()) { localcod &= listV_[*itvrV].locCode(); }
         // If the face belongs to the area, store it.
         if (localcod & sig) {
            if (itT->faceOrientation(i) > 0) { FonB.push_back(vrV); }
            else {// orientation is negative : store vertices in the reverse order
               FonB.push_back(vector<number_t> (vrV.rbegin(),vrV.rend())); }
         }
      }
   }
   return FonB;
}
//-------------------------------------------------------------------------------
//  Private member functions
//-------------------------------------------------------------------------------
/*!
 Write on stream ftex Fig4TeX instructions to create the graphical file and
 the associated TeX box. Selected faces and edges belong to areas of kind TA.
 Longitude psi and latitude theta define the observation direction.
 boxName is the name of the TeX box created.
 */
template<class T_>
void GeomFigureMesh<T_>::createFileNBox(ostream& ftex, const topologicalArea TA,
                                        const float psi, const float theta, const char* boxName,
                                        const string& caption) const {
   ftex << "\\figset proj(psi=" << psi << ", theta=" << theta << ")" << endl;
   ftex << "% 2. Creation of the graphical file" << endl;
   ftex << "\\figdrawbegin{}" << endl;
   printTeXSortedAreaFoE(ftex,TA,psi,theta);
   ftex << "\\figdrawend" << endl;
   ftex << "%" << endl;
   ftex << "% 3. Writing text on the figure" << endl;
   ftex << "\\figvisu{\\figBox" << boxName << "}{" << caption << "}{" << endl;
   // Display first points:
   ftex << "\\figshowpts[1," << nb_main_vertices_by_element_ << "]" << endl;
   ftex << "}" << endl;
}

// class instantiations:
template class GeomFigureMesh<Hexahedron>;
template class GeomFigureMesh<Tetrahedron>;
template class GeomFigureMesh<Quadrangle>;
template class GeomFigureMesh<Triangle>;

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