Skip to content
Snippets Groups Projects
Select Git revision
  • 83a1a852b7405612e563efb02171b0ceb3bc2f5e
  • master default protected
  • hierarchical-basis
  • alphashapes
  • bl
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • fixedMeshIF
  • save_edges
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

robustPredicates.h

Blame
  • HighOrder.cpp 60.03 KiB
    // Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
    //
    // Contributor(s):
    //   Koen Hillewaert
    //
    
    #include <sstream>
    #include <vector>
    #include "GmshConfig.h"
    #include "HighOrder.h"
    #include "MLine.h"
    #include "MTriangle.h"
    #include "MQuadrangle.h"
    #include "MTetrahedron.h"
    #include "MHexahedron.h"
    #include "MPrism.h"
    #include "MPyramid.h"
    #include "GmshMessage.h"
    #include "OS.h"
    #include "fullMatrix.h"
    #include "BasisFactory.h"
    #include "InnerVertexPlacement.h"
    #include "Context.h"
    
    #if defined(HAVE_OPTHOM)
    #include "HighOrderMeshFastCurving.h"
    #include "HighOrderMeshPeriodicity.h"
    #endif
    
    // Functions that help optimizing placement of points on geometry
    
    // The aim here is to build a polynomial representation that consist
    // in polynomial segments of equal length
    
    static double mylength(GEdge *ge, int i, double *u)
    {
      return ge->length(u[i], u[i + 1], 10);
    }
    
    static void myresid(int N, GEdge *ge, double *u, fullVector<double> &r,
                        double *weight = NULL)
    {
      double L[100];
      for(int i = 0; i < N - 1; i++) L[i] = mylength(ge, i, u);
      if(weight)
        for(int i = 0; i < N - 2; i++)
          r(i) = L[i + 1] / weight[i + 1] - L[i] / weight[i];
      else
        for(int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
    }
    
    static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
                                             double *u, double underRelax)
    {
      const double PRECISION = 1.e-6;
      const int MAX_ITER = 50;
      const double eps = (uN - u0) * 1.e-5;
    
      // newton algorithm
      // N is the total number of points (3 for quadratic, 4 for cubic ...)
      // u[0] = u0;
      // u[N-1] = uN;
      // initialize as equidistant in parameter space
      u[0] = u0;
      double du = (uN - u0) / (N - 1);
      for(int i = 1; i < N; i++) {
        u[i] = u[i - 1] + du;
      }
    
      // return true;
    
      // create the tangent matrix
      const int M = N - 2;
      fullMatrix<double> J(M, M);
      fullVector<double> DU(M);
      fullVector<double> R(M);
      fullVector<double> Rp(M);
    
      int iter = 1;
    
      while(iter < MAX_ITER) {
        iter++;
        myresid(N, ge, u, R);
    
        for(int i = 0; i < M; i++) {
          u[i + 1] += eps;
          myresid(N, ge, u, Rp);
          for(int j = 0; j < M; j++) {
            J(i, j) = (Rp(j) - R(j)) / eps;
          }
          u[i + 1] -= eps;
        }
    
        if(M == 1)
          DU(0) = R(0) / J(0, 0);
        else
          J.luSolve(R, DU);
    
        for(int i = 0; i < M; i++) {
          u[i + 1] -= underRelax * DU(i);
        }
    
        if(u[1] < u0) break;
        if(u[N - 2] > uN) break;
    
        double newt_norm = DU.norm();
        if(newt_norm < PRECISION) {
          return true;
        }
      }
      return false;
    }
    
    static bool computeGLLParametersP6(GEdge *ge, double u0, double uN, int N,
                                       double *u, double underRelax)
    {
      static const double GLLQL[7] = {
        -1.000000000000000, -0.830223896278567, -0.468848793470714, 0,
        0.468848793470714,  0.830223896278567,  1.000000000000000};
      double weight[6];
      for(int i = 0; i < 6; ++i) {
        weight[i] = GLLQL[i + 1] - GLLQL[i];
      }
    
      const double PRECISION = 1.e-6;
      const int MAX_ITER = 50;
      const double eps = (uN - u0) * 1.e-5;
    
      // newton algorithm
      // N is the total number of points (3 for quadratic, 4 for cubic ...)
      // u[0] = u0;
      // u[N-1] = uN;
      // initialize as equidistant in parameter space
      u[0] = u0;
      u[N - 1] = uN;
      double uMiddle = .5 * (u0 + uN);
      double du = .5 * (uN - u0);
      for(int i = 1; i < N - 1; i++) {
        u[i] = uMiddle + GLLQL[i] * du;
      }
    
      // create the tangent matrix
      const int M = N - 2;
      fullMatrix<double> J(M, M);
      fullVector<double> DU(M);
      fullVector<double> R(M);
      fullVector<double> Rp(M);
    
      int iter = 1;
    
      while(iter < MAX_ITER) {
        iter++;
        myresid(N, ge, u, R, weight);
    
        for(int i = 0; i < M; i++) {
          u[i + 1] += eps;
          myresid(N, ge, u, Rp, weight);
          for(int j = 0; j < M; j++) {
            J(i, j) = (Rp(j) - R(j)) / eps;
          }
          u[i + 1] -= eps;
        }
    
        if(M == 1)
          DU(0) = R(0) / J(0, 0);
        else
          J.luSolve(R, DU);
    
        for(int i = 0; i < M; i++) {
          u[i + 1] -= underRelax * DU(i);
        }
    
        if(u[1] < u0) break;
        if(u[N - 2] > uN) break;
    
        double newt_norm = DU.norm();
        if(newt_norm < PRECISION) {
          return true;
        }
      }
      return false;
    }
    
    static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
                                             double v0, double vN, SPoint3 &p0,
                                             SPoint3 &pN, int N, bool geodesic,
                                             double *u, double *v)
    {
      const int NI = N - 1;
      u[0] = u0;
      v[0] = v0;
      u[NI] = uN;
      v[NI] = vN;
      const double fact = 1. / (double)NI;
      for(int i = 1; i < NI; i++) {
        const double t = i * fact;
        u[i] = u0 + (uN - u0) * t;
        v[i] = v0 + (vN - v0) * t;
        if(geodesic) {
          SPoint3 pc(t * pN + (1. - t) * p0);
          double guess[2] = {u[i], v[i]};
          GPoint gp = gf->closestPoint(pc, guess);
          if(gp.succeeded()) {
            u[i] = gp.u();
            v[i] = gp.v();
          }
        }
      }
      return true;
    }
    
    static fullMatrix<double> *lob2lagP6 = NULL;
    
    void createMatLob2LagP6()
    {
      const double lobPt[7] = {
        -1.000000000000000, -0.830223896278567, -0.468848793470714, 0,
        0.468848793470714,  0.830223896278567,  1.000000000000000};
      const double lagPt[7] = {
        -1.000000000000000, -0.666666666666666, -0.333333333333333, 0,
        0.333333333333333,  0.666666666666666,  1.000000000000000};
      const int ndofs = 7;
    
      const double monomial[7] = {0, 1, 2, 3, 4, 5, 6};
    
      fullMatrix<double> Vandermonde(ndofs, ndofs);
      for(int i = 0; i < ndofs; i++) {
        for(int j = 0; j < ndofs; j++) {
          Vandermonde(i, j) = pow_int(lobPt[i], monomial[j]);
        }
      }
    
      fullMatrix<double> coefficient(ndofs, ndofs);
      Vandermonde.invert(coefficient);
    
      lob2lagP6 = new fullMatrix<double>(ndofs, ndofs);
    
      for(int i = 0; i < ndofs; i++) {
        for(int j = 0; j < ndofs; j++) {
          Vandermonde(i, j) = pow_int(lagPt[i], monomial[j]);
        }
      }
    
      Vandermonde.mult(coefficient, (*lob2lagP6));
    }
    
    // Creation of high-order edge vertices
    
    static bool getEdgeVerticesOnGeo(GEdge *ge, MVertex *v0, MVertex *v1,
                                     std::vector<MVertex *> &ve, int nPts = 1)
    {
      static bool GLLquad = false;
      static const double relaxFail = 1e-2;
      double u0 = 0., u1 = 0., US[100];
      bool reparamOK = reparamMeshVertexOnEdge(v0, ge, u0);
      if(ge->periodic(0) && ge->getEndVertex() &&
         ge->getEndVertex()->getNumMeshVertices() > 0 &&
         v1 == ge->getEndVertex()->mesh_vertices[0])
        u1 = ge->parBounds(0).high();
      else
        reparamOK &= reparamMeshVertexOnEdge(v1, ge, u1);
    
      if(reparamOK) {
        double uMin = std::min(u0, u1), uMax = std::max(u0, u1);
        bool failed = true;
        double relax = 1.;
        while(failed && (relax > relaxFail)) {
          if(GLLquad)
            failed = !computeGLLParametersP6(ge, uMin, uMax, nPts + 2, US, relax);
          else
            failed =
              !computeEquidistantParameters(ge, uMin, uMax, nPts + 2, US, relax);
          relax *= 0.5;
        }
        if(failed) {
          Msg::Warning(
            "Failed to compute equidistant parameters (relax = %g, value = %g) "
            "for edge %d-%d parametrized with %g %g on curve %d",
            relax, US[1], v0->getNum(), v1->getNum(), u0, u1, ge->tag());
          US[0] = uMin;
          const double du = (uMax - uMin) / (nPts + 1);
          for(int i = 1; i <= nPts; i++) US[i] = US[i - 1] + du;
        }
      }
      else
        Msg::Error("Cannot reparametrize a mesh node in high order meshing");
    
      if(!reparamOK) return false;
    
      if(GLLquad) {
        fullMatrix<double> M(7, 3);
        M(0, 0) = u0 < u1 ? v0->x() : v1->x();
        M(0, 1) = u0 < u1 ? v0->y() : v1->y();
        M(0, 2) = u0 < u1 ? v0->z() : v1->z();
        M(6, 0) = u0 < u1 ? v1->x() : v0->x();
        M(6, 1) = u0 < u1 ? v1->y() : v0->y();
        M(6, 2) = u0 < u1 ? v1->z() : v0->z();
        for(int j = 0; j < nPts; j++) {
          int count = u0 < u1 ? j + 1 : nPts + 1 - (j + 1);
          GPoint pc = ge->point(US[count]);
          M(j + 1, 0) = pc.x();
          M(j + 1, 1) = pc.y();
          M(j + 1, 2) = pc.z();
        }
        fullMatrix<double> Mlag(7, 3);
        if(!lob2lagP6) createMatLob2LagP6();
        lob2lagP6->mult(M, Mlag);
    
        for(int j = 0; j < nPts; j++) {
          MVertex *v;
          int count = u0 < u1 ? j + 1 : nPts + 1 - (j + 1);
          // FIXME US[count] false!!!
          v = new MEdgeVertex(Mlag(count, 0), Mlag(count, 1), Mlag(count, 2), ge,
                              US[count]);
          // this destroys the ordering of the mesh vertices on the edge
          ve.push_back(v);
        }
      }
      else {
        for(int j = 0; j < nPts; j++) {
          MVertex *v;
          int count = u0 < u1 ? j + 1 : nPts + 1 - (j + 1);
          GPoint pc = ge->point(US[count]);
          v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[count]);
          // this destroys the ordering of the mesh vertices on the edge
          ve.push_back(v);
        }
      }
    
      //  GLLquad = false;
      return true;
    }
    
    static bool getEdgeVerticesOnGeo(GFace *gf, MVertex *v0, MVertex *v1,
                                     std::vector<MVertex *> &ve, int nPts = 1)
    {
      SPoint2 p0, p1;
      double US[100], VS[100];
      bool reparamOK = reparamMeshEdgeOnFace(v0, v1, gf, p0, p1);
      if(reparamOK) {
        SPoint3 pnt0, pnt1;
        if(nPts >= 30)
          computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], pnt0, pnt1,
                                       nPts + 2, false, US, VS);
        else {
          pnt0 = v0->point();
          pnt1 = v1->point();
          // FIXME: using the geodesic is sometimes a bad idea when the edge is "far
          // away" from the surface (e.g. on the diameter of a circle!)
          computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], pnt0, pnt1,
                                       nPts + 2, true, US, VS);
        }
      }
      else {
        Msg::Error("Cannot reparametrize mesh edge %lu-%lu on surface %d",
                   v0->getNum(), v1->getNum(), gf->tag());
        return false;
      }
    
      for(int j = 0; j < nPts; j++) {
        GPoint pc = gf->point(US[j + 1], VS[j + 1]);
        MVertex *v =
          new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j + 1], VS[j + 1]);
        ve.push_back(v);
      }
    
      return true;
    }
    
    static void interpVerticesInExistingEdge(GEntity *ge, const MElement *edgeEl,
                                             std::vector<MVertex *> &veEdge,
                                             int nPts)
    {
      fullMatrix<double> points;
      points = edgeEl->getFunctionSpace(nPts + 1)->points;
      for(int k = 2; k < nPts + 2; k++) {
        SPoint3 pos;
        edgeEl->pnt(points(k, 0), 0., 0., pos);
        MVertex *v = new MVertex(pos.x(), pos.y(), pos.z(), ge);
        veEdge.push_back(v);
      }
    }
    
    inline static bool getMinMaxVert(MVertex *v0, MVertex *v1, MVertex *&vMin,
                                     MVertex *&vMax)
    {
      const bool increasing = (v0->getNum() < v1->getNum());
      if(increasing) {
        vMin = v0;
        vMax = v1;
      }
      else {
        vMin = v1;
        vMax = v0;
      }
      return increasing;
    }
    
    // Get new interior vertices for a 1D element
    static void getEdgeVertices(GEdge *ge, MElement *ele,
                                std::vector<MVertex *> &ve,
                                std::vector<MVertex *> &newHOVert,
                                edgeContainer &edgeVertices, bool linear,
                                int nPts = 1)
    {
      if(!ge->haveParametrization()) linear = true;
    
      std::vector<MVertex *> veOld;
      ele->getVertices(veOld);
      MVertex *vMin, *vMax;
      const bool increasing = getMinMaxVert(veOld[0], veOld[1], vMin, vMax);
      std::pair<MVertex *, MVertex *> p(vMin, vMax);
      std::vector<MVertex *> veEdge;
      // Get vertices on geometry if asked
      bool gotVertOnGeo =
        linear ? false : getEdgeVerticesOnGeo(ge, veOld[0], veOld[1], veEdge, nPts);
      // If not on geometry, create from mesh interpolation
      if(!gotVertOnGeo) interpVerticesInExistingEdge(ge, ele, veEdge, nPts);
      newHOVert.insert(newHOVert.end(), veEdge.begin(), veEdge.end());
      if(edgeVertices.count(p) == 0) {
        if(increasing) // Add newly created vertices to list
          edgeVertices[p].insert(edgeVertices[p].end(), veEdge.begin(),
                                 veEdge.end());
        else
          edgeVertices[p].insert(edgeVertices[p].end(), veEdge.rbegin(),
                                 veEdge.rend());
      }
      else if(p.first != p.second) {
        // Vertices already exist and edge is not a degenerated edge
        Msg::Error("Mesh edges from different entities share nodes: create a finer mesh "
                   "(curve involved: %d)", ge->tag());
      }
      ve.insert(ve.end(), veEdge.begin(), veEdge.end());
    }
    
    // Get new interior vertices for an edge in a 2D element
    static void getEdgeVertices(GFace *gf, MElement *ele,
                                std::vector<MVertex *> &ve,
                                std::vector<MVertex *> &newHOVert,
                                edgeContainer &edgeVertices, bool linear,
                                int nPts = 1)
    {
      if(!gf->haveParametrization()) linear = true;
    
      for(int i = 0; i < ele->getNumEdges(); i++) {
        std::vector<MVertex *> veOld;
        ele->getEdgeVertices(i, veOld);
        MVertex *vMin, *vMax;
        const bool increasing = getMinMaxVert(veOld[0], veOld[1], vMin, vMax);
        std::pair<MVertex *, MVertex *> p(vMin, vMax);
        std::vector<MVertex *> veEdge;
    
        edgeContainer::iterator eIter = edgeVertices.find(p);
    
        if(eIter != edgeVertices.end()) { // Vertices already exist
          std::vector<MVertex *> &eVtcs = eIter->second;
          if(increasing)
            veEdge.assign(eVtcs.begin(), eVtcs.end());
          else
            veEdge.assign(eVtcs.rbegin(), eVtcs.rend());
        }
        else { // Vertices do not exist, create them
          // Get vertices on geometry if asked
          bool gotVertOnGeo =
            linear ? false :
                     getEdgeVerticesOnGeo(gf, veOld[0], veOld[1], veEdge, nPts);
          if(!gotVertOnGeo) {
            // If not on geometry, create from mesh interpolation
            const MLineN edgeEl(veOld, ele->getPolynomialOrder());
            interpVerticesInExistingEdge(gf, &edgeEl, veEdge, nPts);
          }
          newHOVert.insert(newHOVert.end(), veEdge.begin(), veEdge.end());
    
          std::vector<MVertex *> &eVtcs = edgeVertices[p];
    
          if(increasing) // Add newly created vertices to list
            eVtcs.insert(eVtcs.end(), veEdge.begin(), veEdge.end());
          else
            eVtcs.insert(eVtcs.end(), veEdge.rbegin(), veEdge.rend());
        }
        ve.insert(ve.end(), veEdge.begin(), veEdge.end());
      }
    }
    
    // Get new interior vertices for an edge in a 3D element
    static void getEdgeVertices(GRegion *gr, MElement *ele,
                                std::vector<MVertex *> &ve,
                                std::vector<MVertex *> &newHOVert,
                                edgeContainer &edgeVertices, int nPts = 1)
    {
      for(int i = 0; i < ele->getNumEdges(); i++) {
        std::vector<MVertex *> veOld;
        ele->getEdgeVertices(i, veOld);
        MVertex *vMin, *vMax;
        const bool increasing = getMinMaxVert(veOld[0], veOld[1], vMin, vMax);
        std::pair<MVertex *, MVertex *> p(vMin, vMax);
        std::vector<MVertex *> veEdge;
        if(edgeVertices.count(p)) { // Vertices already exist
          if(increasing)
            veEdge.assign(edgeVertices[p].begin(), edgeVertices[p].end());
          else
            veEdge.assign(edgeVertices[p].rbegin(), edgeVertices[p].rend());
        }
        else { // Vertices do not exist, create them
          const MLineN edgeEl(veOld, ele->getPolynomialOrder());
          interpVerticesInExistingEdge(gr, &edgeEl, veEdge, nPts);
          newHOVert.insert(newHOVert.end(), veEdge.begin(), veEdge.end());
          if(increasing) // Add newly created vertices to list
            edgeVertices[p].insert(edgeVertices[p].end(), veEdge.begin(),
                                   veEdge.end());
          else
            edgeVertices[p].insert(edgeVertices[p].end(), veEdge.rbegin(),
                                   veEdge.rend());
        }
        ve.insert(ve.end(), veEdge.begin(), veEdge.end());
      }
    }
    
    // Creation of high-order face vertices
    
    static void reorientTrianglePoints(std::vector<MVertex *> &vtcs,
                                       int orientation, bool swap)
    {
      int nbPts = vtcs.size();
      if(nbPts <= 1) return;
      std::vector<MVertex *> tmp(nbPts);
      int interiorOrder = (int)((sqrt(1. + 8. * nbPts) - 3) / 2);
      int pos = 0;
      for(int o = interiorOrder; o > 0; o -= 3) {
        if(swap) {
          tmp[pos] = vtcs[pos];
          tmp[pos + 1] = vtcs[pos + 2];
          tmp[pos + 2] = vtcs[pos + 1];
          for(int i = 0; i < 3 * (o - 1); i++)
            tmp[pos + 3 + i] = vtcs[pos + 3 * o - i - 1];
        }
        else {
          for(int i = 0; i < 3 * o; i++) tmp[pos + i] = vtcs[pos + i];
        }
        for(int i = 0; i < 3; i++) {
          int ri = (i + orientation) % 3;
          vtcs[pos + ri] = tmp[pos + i];
          for(int j = 0; j < o - 1; j++)
            vtcs[pos + 3 + (o - 1) * ri + j] = tmp[pos + 3 + (o - 1) * i + j];
        }
        pos += 3 * o;
      }
    }
    
    static void reorientQuadPoints(std::vector<MVertex *> &vtcs, int orientation,
                                   bool swap, int order)
    {
      int nbPts = vtcs.size();
      if(nbPts <= 1) return;
      std::vector<MVertex *> tmp(nbPts);
    
      int start = 0;
      while(1) {
        // CORNERS
        int index = 0;
        if(order == 0) {
          start++;
        }
        else {
          int i1(0), i2(0), i3(0), i4(0);
          if(!swap) {
            if(orientation == 0) {
              i1 = 0;
              i2 = 1;
              i3 = 2;
              i4 = 3;
            }
            else if(orientation == 1) {
              i1 = 3;
              i2 = 0;
              i3 = 1;
              i4 = 2;
            }
            else if(orientation == 2) {
              i1 = 2;
              i2 = 3;
              i3 = 0;
              i4 = 1;
            }
            else if(orientation == 3) {
              i1 = 1;
              i2 = 2;
              i3 = 3;
              i4 = 0;
            }
          }
          else {
            if(orientation == 0) {
              i1 = 0;
              i2 = 3;
              i3 = 2;
              i4 = 1;
            }
            else if(orientation == 3) {
              i1 = 3;
              i2 = 2;
              i3 = 1;
              i4 = 0;
            }
            else if(orientation == 2) {
              i1 = 2;
              i2 = 1;
              i3 = 0;
              i4 = 3;
            }
            else if(orientation == 1) {
              i1 = 1;
              i2 = 0;
              i3 = 3;
              i4 = 2;
            }
          }
    
          int indices[4] = {i1, i2, i3, i4};
          for(int i = 0; i < 4; i++) tmp[i] = vtcs[start + indices[i]];
          for(int i = 0; i < 4; i++) vtcs[start + i] = tmp[i];
    
          // EDGES
          start += 4;
          for(int iEdge = 0; iEdge < 4; iEdge++) {
            int p1 = indices[iEdge];
            int p2 = indices[(iEdge + 1) % 4];
            int nbP = order - 1;
            if(p1 == 0 && p2 == 1) {
              for(int i = start + 0 * nbP; i < start + 1 * nbP; i++)
                tmp[index++] = vtcs[i];
            }
            else if(p1 == 1 && p2 == 2) {
              for(int i = start + 1 * nbP; i < start + 2 * nbP; i++)
                tmp[index++] = vtcs[i];
            }
            else if(p1 == 2 && p2 == 3) {
              for(int i = start + 2 * nbP; i < start + 3 * nbP; i++)
                tmp[index++] = vtcs[i];
            }
            else if(p1 == 3 && p2 == 0) {
              for(int i = start + 3 * nbP; i < start + 4 * nbP; i++)
                tmp[index++] = vtcs[i];
            }
            else if(p1 == 1 && p2 == 0) {
              for(int i = start + 1 * nbP - 1; i >= start + 0 * nbP; i--)
                tmp[index++] = vtcs[i];
            }
            else if(p1 == 2 && p2 == 1) {
              for(int i = start + 2 * nbP - 1; i >= start + 1 * nbP; i--)
                tmp[index++] = vtcs[i];
            }
            else if(p1 == 3 && p2 == 2) {
              for(int i = start + 3 * nbP - 1; i >= start + 2 * nbP; i--)
                tmp[index++] = vtcs[i];
            }
            else if(p1 == 0 && p2 == 3) {
              for(int i = start + 4 * nbP - 1; i >= start + 3 * nbP; i--)
                tmp[index++] = vtcs[i];
            }
            else
              Msg::Error("Something wrong in reorientQuadPoints");
          }
          for(int i = 0; i < index; i++) vtcs[start + i] = tmp[i];
          start += index;
        }
    
        order -= 2;
        if(start >= (int)vtcs.size()) break;
      }
    }
    
    static void getFaceVerticesOnGeo(GFace *gf,
                                     const fullMatrix<double> &coefficients,
                                     const std::vector<MVertex *> &vertices,
                                     std::vector<MVertex *> &vf)
    {
      SPoint2 pts[1000];
      bool reparamOK = true;
      for(std::size_t k = 0; k < vertices.size(); ++k)
        reparamOK &= reparamMeshVertexOnFace(vertices[k], gf, pts[k]);
      for(int k = 0; k < coefficients.size1(); k++) {
        double X(0), Y(0), Z(0), GUESS[2] = {0, 0};
        for(int j = 0; j < coefficients.size2(); j++) {
          MVertex *vt = vertices[j];
          X += coefficients(k, j) * vt->x();
          Y += coefficients(k, j) * vt->y();
          Z += coefficients(k, j) * vt->z();
          if(reparamOK) {
            GUESS[0] += coefficients(k, j) * pts[j][0];
            GUESS[1] += coefficients(k, j) * pts[j][1];
          }
        }
        MVertex *v;
        if(reparamOK) {
          // GPoint gp = gf->point(SPoint2(GUESS[0], GUESS[1]));
          // closest point is not necessary (slow and for high quality HO
          // meshes it should be optimized anyway afterwards + closest point
          // is still buggy (e.g. BFGS for a plane "Ruled Surface")
          GPoint gp = gf->closestPoint(SPoint3(X, Y, Z), GUESS);
          // AJ: ClosestPoint is absolutely necessary when the parameterization
          // is degenerate...
          if(gp.g()) {
            v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
          }
          else {
            v = new MVertex(X, Y, Z, gf);
          }
        }
        else {
          GPoint gp = gf->closestPoint(SPoint3(X, Y, Z), GUESS);
          if(gp.succeeded())
            v = new MVertex(gp.x(), gp.y(), gp.z(), gf);
          else
            v = new MVertex(X, Y, Z, gf);
        }
        vf.push_back(v);
      }
    }
    
    static void interpVerticesInExistingFace(GEntity *ge,
                                             const fullMatrix<double> &coefficients,
                                             const std::vector<MVertex *> &vertices,
                                             std::vector<MVertex *> &vFace)
    {
      for(int k = 0; k < coefficients.size1(); k++) {
        double x(0), y(0), z(0);
        for(int j = 0; j < coefficients.size2(); j++) {
          MVertex *v = vertices[j];
          x += coefficients(k, j) * v->x();
          y += coefficients(k, j) * v->y();
          z += coefficients(k, j) * v->z();
        }
        vFace.push_back(new MVertex(x, y, z, ge));
      }
    }
    
    // Get new interior vertices for a 2D element
    static void getFaceVertices(GFace *gf, MElement *ele,
                                std::vector<MVertex *> &newVertices,
                                std::vector<MVertex *> &newHOVert,
                                faceContainer &faceVertices, bool linear,
                                int nPts = 1)
    {
      if(!gf->haveParametrization()) linear = true;
    
      std::vector<MVertex *> boundaryVertices;
      {
        int nCorner = ele->getNumPrimaryVertices();
        boundaryVertices.reserve(nCorner + newVertices.size());
        ele->getVertices(boundaryVertices);
        boundaryVertices.resize(nCorner);
        boundaryVertices.insert(boundaryVertices.end(), newVertices.begin(),
                                newVertices.end());
      }
      int type = ele->getType();
      fullMatrix<double> *coefficients = getInnerVertexPlacement(type, nPts + 1);
      std::vector<MVertex *> vFace;
      if(!linear) { // Get vertices on geometry if asked...
        getFaceVerticesOnGeo(gf, *coefficients, boundaryVertices, vFace);
      }
      else { // ... otherwise, create from mesh interpolation
        interpVerticesInExistingFace(gf, *coefficients, boundaryVertices, vFace);
      }
    
      MFace face = ele->getFace(0);
      faceVertices[face].insert(faceVertices[face].end(), vFace.begin(),
                                vFace.end());
      newVertices.insert(newVertices.end(), vFace.begin(), vFace.end());
      newHOVert.insert(newHOVert.end(), vFace.begin(), vFace.end());
    }
    
    static int retrieveFaceBoundaryVertices(int k, int type, int nPts,
                                            const std::vector<MVertex *> &vCorner,
                                            const std::vector<MVertex *> &vEdges,
                                            std::vector<MVertex *> &v)
    {
      v.clear();
      int nCorner;
      switch(type) {
      case TYPE_TET:
        v.push_back(vCorner[MTetrahedron::faces_tetra(k, 0)]);
        v.push_back(vCorner[MTetrahedron::faces_tetra(k, 1)]);
        v.push_back(vCorner[MTetrahedron::faces_tetra(k, 2)]);
        for(int i = 0; i < 3; ++i) {
          int edge = MTetrahedron::faces2edge_tetra(k, i);
          int n = std::abs(edge) - 1;
          v.insert(v.end(), vEdges.begin() + n * nPts,
                   vEdges.begin() + (n + 1) * nPts);
          if(edge < 0) std::reverse(v.end() - nPts, v.end());
        }
        return TYPE_TRI;
      case TYPE_HEX:
        v.push_back(vCorner[MHexahedron::faces_hexa(k, 0)]);
        v.push_back(vCorner[MHexahedron::faces_hexa(k, 1)]);
        v.push_back(vCorner[MHexahedron::faces_hexa(k, 2)]);
        v.push_back(vCorner[MHexahedron::faces_hexa(k, 3)]);
        for(int i = 0; i < 4; ++i) {
          int edge = MHexahedron::faces2edge_hexa(k, i);
          int n = std::abs(edge) - 1;
          v.insert(v.end(), vEdges.begin() + n * nPts,
                   vEdges.begin() + (n + 1) * nPts);
          if(edge < 0) std::reverse(v.end() - nPts, v.end());
        }
        return TYPE_QUA;
      case TYPE_PRI:
        nCorner = k < 2 ? 3 : 4;
        v.push_back(vCorner[MPrism::faces_prism(k, 0)]);
        v.push_back(vCorner[MPrism::faces_prism(k, 1)]);
        v.push_back(vCorner[MPrism::faces_prism(k, 2)]);
        if(nCorner == 4) v.push_back(vCorner[MPrism::faces_prism(k, 3)]);
        for(int i = 0; i < nCorner; ++i) {
          int edge = MPrism::faces2edge_prism(k, i);
          int n = std::abs(edge) - 1;
          v.insert(v.end(), vEdges.begin() + n * nPts,
                   vEdges.begin() + (n + 1) * nPts);
          if(edge < 0) std::reverse(v.end() - nPts, v.end());
        }
        return nCorner == 3 ? TYPE_TRI : TYPE_QUA;
      case TYPE_PYR:
        nCorner = k < 4 ? 3 : 4;
        v.push_back(vCorner[MPyramid::faces_pyramid(k, 0)]);
        v.push_back(vCorner[MPyramid::faces_pyramid(k, 1)]);
        v.push_back(vCorner[MPyramid::faces_pyramid(k, 2)]);
        if(nCorner == 4) v.push_back(vCorner[MPyramid::faces_pyramid(k, 3)]);
        for(int i = 0; i < nCorner; ++i) {
          int edge = MPyramid::faces2edge_pyramid(k, i);
          int n = std::abs(edge) - 1;
          v.insert(v.end(), vEdges.begin() + n * nPts,
                   vEdges.begin() + (n + 1) * nPts);
          if(edge < 0) std::reverse(v.end() - nPts, v.end());
        }
        return nCorner == 3 ? TYPE_TRI : TYPE_QUA;
      default: return -1;
      }
    }
    
    // Get new face (excluding edge) vertices for a face of a 3D element
    static void getFaceVertices(GRegion *gr, MElement *ele,
                                std::vector<MVertex *> &newVertices,
                                std::vector<MVertex *> &newHOVert,
                                faceContainer &faceVertices, int nPts = 1)
    {
      std::vector<MVertex *> vCorner;
      ele->getVertices(vCorner);
      // NB: We can get more than corner vertices but we use only corners
    
      for(int i = 0; i < ele->getNumFaces(); i++) {
        MFace face = ele->getFace(i);
        std::vector<MVertex *> vFace;
        faceContainer::iterator fIter = faceVertices.find(face);
        if(fIter != faceVertices.end()) { // Vertices already exist
          std::vector<MVertex *> vtcs = fIter->second;
          int orientation;
          bool swap;
          if(fIter->first.computeCorrespondence(face, orientation, swap)) {
            // Check correspondence and apply permutation if needed
            if(face.getNumVertices() == 3 && nPts > 1)
              reorientTrianglePoints(vtcs, orientation, swap);
            else if(face.getNumVertices() == 4)
              reorientQuadPoints(vtcs, orientation, swap, nPts - 1);
          }
          else
            Msg::Error(
              "Error in face lookup for retrieval of high order face nodes");
          vFace.assign(vtcs.begin(), vtcs.end());
        }
        else { // Vertices do not exist, create them by interpolation
          std::vector<MVertex *> faceBoundaryVertices;
          int type = retrieveFaceBoundaryVertices(
            i, ele->getType(), nPts, vCorner, newVertices, faceBoundaryVertices);
          fullMatrix<double> *coefficients =
            getInnerVertexPlacement(type, nPts + 1);
          interpVerticesInExistingFace(gr, *coefficients, faceBoundaryVertices,
                                       vFace);
          newHOVert.insert(newHOVert.end(), vFace.begin(), vFace.end());
          faceVertices[face].insert(faceVertices[face].end(), vFace.begin(),
                                    vFace.end());
        }
        newVertices.insert(newVertices.end(), vFace.begin(), vFace.end());
      }
    }
    
    // Get new interior vertices for a 3D element
    static void getVolumeVertices(GRegion *gr, MElement *ele,
                                  std::vector<MVertex *> &newVertices,
                                  std::vector<MVertex *> &newHOVert, int nPts = 1)
    {
      std::vector<MVertex *> boundaryVertices;
      {
        int nCorner = ele->getNumPrimaryVertices();
        boundaryVertices.reserve(nCorner + newVertices.size());
        ele->getVertices(boundaryVertices);
        boundaryVertices.resize(nCorner);
        boundaryVertices.insert(boundaryVertices.end(), newVertices.begin(),
                                newVertices.end());
      }
      int type = ele->getType();
      fullMatrix<double> &coefficients = *getInnerVertexPlacement(type, nPts + 1);
    
      for(int k = 0; k < coefficients.size1(); k++) {
        double x(0), y(0), z(0);
        for(int j = 0; j < coefficients.size2(); j++) {
          MVertex *v = boundaryVertices[j];
          x += coefficients(k, j) * v->x();
          y += coefficients(k, j) * v->y();
          z += coefficients(k, j) * v->z();
        }
        MVertex *v = new MVertex(x, y, z, gr);
        newHOVert.push_back(v);
        newVertices.push_back(v);
      }
    }
    
    // Creation of high-order elements
    
    static void setHighOrder(GEdge *ge, std::vector<MVertex *> &newHOVert,
                             edgeContainer &edgeVertices, bool linear,
                             int nbPts = 1)
    {
      std::vector<MLine *> lines2;
      for(std::size_t i = 0; i < ge->lines.size(); i++) {
        MLine *l = ge->lines[i];
        std::vector<MVertex *> ve;
        getEdgeVertices(ge, l, ve, newHOVert, edgeVertices, linear, nbPts);
        if(nbPts == 1)
          lines2.push_back(
            new MLine3(l->getVertex(0), l->getVertex(1), ve[0], l->getPartition()));
        else
          lines2.push_back(
            new MLineN(l->getVertex(0), l->getVertex(1), ve, l->getPartition()));
        delete l;
      }
      ge->lines = lines2;
      ge->deleteVertexArrays();
    }
    
    static MTriangle *setHighOrder(MTriangle *t, GFace *gf,
                                   std::vector<MVertex *> &newHOVert,
                                   edgeContainer &edgeVertices,
                                   faceContainer &faceVertices, bool linear,
                                   bool incomplete, int nPts)
    {
      std::vector<MVertex *> v;
      getEdgeVertices(gf, t, v, newHOVert, edgeVertices, linear, nPts);
      if(nPts == 1) {
        return new MTriangle6(t->getVertex(0), t->getVertex(1), t->getVertex(2),
                              v[0], v[1], v[2], 0, t->getPartition());
      }
      else {
        if(!incomplete)
          getFaceVertices(gf, t, v, newHOVert, faceVertices, linear, nPts);
        return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2), v,
                              nPts + 1, 0, t->getPartition());
      }
    }
    
    static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
                                     std::vector<MVertex *> &newHOVert,
                                     edgeContainer &edgeVertices,
                                     faceContainer &faceVertices, bool linear,
                                     bool incomplete, int nPts)
    {
      std::vector<MVertex *> v;
      getEdgeVertices(gf, q, v, newHOVert, edgeVertices, linear, nPts);
      if(incomplete) {
        if(nPts == 1) {
          return new MQuadrangle8(q->getVertex(0), q->getVertex(1), q->getVertex(2),
                                  q->getVertex(3), v[0], v[1], v[2], v[3], 0,
                                  q->getPartition());
        }
        else {
          return new MQuadrangleN(q->getVertex(0), q->getVertex(1), q->getVertex(2),
                                  q->getVertex(3), v, nPts + 1, 0,
                                  q->getPartition());
        }
      }
      else {
        getFaceVertices(gf, q, v, newHOVert, faceVertices, linear, nPts);
        if(nPts == 1) {
          return new MQuadrangle9(q->getVertex(0), q->getVertex(1), q->getVertex(2),
                                  q->getVertex(3), v[0], v[1], v[2], v[3], v[4], 0,
                                  q->getPartition());
        }
        else {
          return new MQuadrangleN(q->getVertex(0), q->getVertex(1), q->getVertex(2),
                                  q->getVertex(3), v, nPts + 1, 0,
                                  q->getPartition());
        }
      }
    }
    
    static void setHighOrder(GFace *gf, std::vector<MVertex *> &newHOVert,
                             edgeContainer &edgeVertices,
                             faceContainer &faceVertices, bool linear,
                             bool incomplete, int nPts = 1)
    {
      std::vector<MTriangle *> triangles2;
      for(std::size_t i = 0; i < gf->triangles.size(); i++) {
        MTriangle *t = gf->triangles[i];
        MTriangle *tNew = setHighOrder(t, gf, newHOVert, edgeVertices, faceVertices,
                                       linear, incomplete, nPts);
        triangles2.push_back(tNew);
        delete t;
      }
      gf->triangles = triangles2;
    
      std::vector<MQuadrangle *> quadrangles2;
      for(std::size_t i = 0; i < gf->quadrangles.size(); i++) {
        MQuadrangle *q = gf->quadrangles[i];
        MQuadrangle *qNew = setHighOrder(q, gf, newHOVert, edgeVertices,
                                         faceVertices, linear, incomplete, nPts);
        quadrangles2.push_back(qNew);
        delete q;
      }
      gf->quadrangles = quadrangles2;
      gf->deleteVertexArrays();
    }
    
    static MTetrahedron *setHighOrder(MTetrahedron *t, GRegion *gr,
                                      std::vector<MVertex *> &newHOVert,
                                      edgeContainer &edgeVertices,
                                      faceContainer &faceVertices, bool incomplete,
                                      int nPts)
    {
      std::vector<MVertex *> v;
      getEdgeVertices(gr, t, v, newHOVert, edgeVertices, nPts);
      if(nPts == 1) {
        return new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
                                  t->getVertex(3), v[0], v[1], v[2], v[3], v[4],
                                  v[5], 0, t->getPartition());
      }
      else {
        if(!incomplete) {
          getFaceVertices(gr, t, v, newHOVert, faceVertices, nPts);
          getVolumeVertices(gr, t, v, newHOVert, nPts);
        }
        return new MTetrahedronN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
                                 t->getVertex(3), v, nPts + 1, 0,
                                 t->getPartition());
      }
    }
    
    static MHexahedron *setHighOrder(MHexahedron *h, GRegion *gr,
                                     std::vector<MVertex *> &newHOVert,
                                     edgeContainer &edgeVertices,
                                     faceContainer &faceVertices, bool incomplete,
                                     int nPts)
    {
      std::vector<MVertex *> v;
      getEdgeVertices(gr, h, v, newHOVert, edgeVertices, nPts);
      if(incomplete) {
        if(nPts == 1) {
          return new MHexahedron20(
            h->getVertex(0), h->getVertex(1), h->getVertex(2), h->getVertex(3),
            h->getVertex(4), h->getVertex(5), h->getVertex(6), h->getVertex(7),
            v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10],
            v[11], 0, h->getPartition());
        }
        else {
          return new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
                                  h->getVertex(3), h->getVertex(4), h->getVertex(5),
                                  h->getVertex(6), h->getVertex(7), v, nPts + 1, 0,
                                  h->getPartition());
        }
      }
      else {
        getFaceVertices(gr, h, v, newHOVert, faceVertices, nPts);
        getVolumeVertices(gr, h, v, newHOVert, nPts);
        if(nPts == 1) {
          return new MHexahedron27(
            h->getVertex(0), h->getVertex(1), h->getVertex(2), h->getVertex(3),
            h->getVertex(4), h->getVertex(5), h->getVertex(6), h->getVertex(7),
            v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10],
            v[11], v[12], v[13], v[14], v[15], v[16], v[17], v[18], 0,
            h->getPartition());
        }
        else {
          return new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
                                  h->getVertex(3), h->getVertex(4), h->getVertex(5),
                                  h->getVertex(6), h->getVertex(7), v, nPts + 1, 0,
                                  h->getPartition());
        }
      }
    }
    
    static MPrism *setHighOrder(MPrism *p, GRegion *gr,
                                std::vector<MVertex *> &newHOVert,
                                edgeContainer &edgeVertices,
                                faceContainer &faceVertices, bool incomplete,
                                int nPts)
    {
      std::vector<MVertex *> v;
      getEdgeVertices(gr, p, v, newHOVert, edgeVertices, nPts);
      if(incomplete) {
        if(nPts == 1) {
          return new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2),
                              p->getVertex(3), p->getVertex(4), p->getVertex(5),
                              v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8],
                              0, p->getPartition());
        }
        else {
          return new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
                             p->getVertex(3), p->getVertex(4), p->getVertex(5), v,
                             nPts + 1, 0, p->getPartition());
        }
      }
      else {
        getFaceVertices(gr, p, v, newHOVert, faceVertices, nPts);
        if(nPts == 1) {
          return new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2),
                              p->getVertex(3), p->getVertex(4), p->getVertex(5),
                              v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8],
                              v[9], v[10], v[11], 0, p->getPartition());
        }
        else {
          getVolumeVertices(gr, p, v, newHOVert, nPts);
          return new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
                             p->getVertex(3), p->getVertex(4), p->getVertex(5), v,
                             nPts + 1, 0, p->getPartition());
        }
      }
    }
    
    static MPyramid *setHighOrder(MPyramid *p, GRegion *gr,
                                  std::vector<MVertex *> &newHOVert,
                                  edgeContainer &edgeVertices,
                                  faceContainer &faceVertices, bool incomplete,
                                  int nPts)
    {
      std::vector<MVertex *> v;
      getEdgeVertices(gr, p, v, newHOVert, edgeVertices, nPts);
      if(!incomplete) {
        getFaceVertices(gr, p, v, newHOVert, faceVertices, nPts);
        if(nPts > 1) {
          getVolumeVertices(gr, p, v, newHOVert, nPts);
        }
      }
      return new MPyramidN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
                           p->getVertex(3), p->getVertex(4), v, nPts + 1, 0,
                           p->getPartition());
    }
    
    static void setHighOrder(GRegion *gr, std::vector<MVertex *> &newHOVert,
                             edgeContainer &edgeVertices,
                             faceContainer &faceVertices, bool incomplete,
                             int nPts = 1)
    {
      std::vector<MTetrahedron *> tetrahedra2;
      for(std::size_t i = 0; i < gr->tetrahedra.size(); i++) {
        MTetrahedron *t = gr->tetrahedra[i];
        MTetrahedron *tNew = setHighOrder(t, gr, newHOVert, edgeVertices,
                                          faceVertices, incomplete, nPts);
        tetrahedra2.push_back(tNew);
        delete t;
      }
      gr->tetrahedra = tetrahedra2;
    
      std::vector<MHexahedron *> hexahedra2;
      for(std::size_t i = 0; i < gr->hexahedra.size(); i++) {
        MHexahedron *h = gr->hexahedra[i];
        MHexahedron *hNew = setHighOrder(h, gr, newHOVert, edgeVertices,
                                         faceVertices, incomplete, nPts);
        hexahedra2.push_back(hNew);
        delete h;
      }
      gr->hexahedra = hexahedra2;
    
      std::vector<MPrism *> prisms2;
      for(std::size_t i = 0; i < gr->prisms.size(); i++) {
        MPrism *p = gr->prisms[i];
        MPrism *pNew = setHighOrder(p, gr, newHOVert, edgeVertices, faceVertices,
                                    incomplete, nPts);
        prisms2.push_back(pNew);
        delete p;
      }
      gr->prisms = prisms2;
    
      std::vector<MPyramid *> pyramids2;
      for(std::size_t i = 0; i < gr->pyramids.size(); i++) {
        MPyramid *p = gr->pyramids[i];
        MPyramid *pNew = setHighOrder(p, gr, newHOVert, edgeVertices, faceVertices,
                                      incomplete, nPts);
        pyramids2.push_back(pNew);
        delete p;
      }
      gr->pyramids = pyramids2;
    
      gr->deleteVertexArrays();
    }
    
    // High-level functions
    
    template <class T>
    static void setFirstOrder(GEntity *e, std::vector<T *> &elements,
                              bool onlyVisible)
    {
      if(onlyVisible && !e->getVisibility()) return;
      std::vector<T *> elements1;
      for(std::size_t i = 0; i < elements.size(); i++) {
        T *ele = elements[i];
        int n = ele->getNumPrimaryVertices();
        std::vector<MVertex *> v1;
        v1.reserve(n);
        for(int j = 0; j < n; j++) v1.push_back(ele->getVertex(j));
        elements1.push_back(new T(v1, 0, ele->getPartition()));
        delete ele;
      }
      elements = elements1;
      e->deleteVertexArrays();
    }
    
    static void updateHighOrderVertices(GEntity *e,
                                        const std::vector<MVertex *> &newHOVert,
                                        bool onlyVisible)
    {
      if(onlyVisible && !e->getVisibility()) return;
      std::vector<MVertex *> v1;
      for(std::size_t i = 0; i < e->mesh_vertices.size(); i++) {
        if(e->mesh_vertices[i]->getPolynomialOrder() > 1)
          delete e->mesh_vertices[i];
        else
          v1.push_back(e->mesh_vertices[i]);
      }
      v1.insert(v1.end(), newHOVert.begin(), newHOVert.end());
      e->mesh_vertices = v1;
      e->deleteVertexArrays();
    }
    
    static void updatePeriodicEdgesAndFaces(GModel *m)
    {
      // Edges
    
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
        GEdge *tgt = *it;
    
        // non complete periodic info (e.g. through extrusion)
        if(tgt->vertexCounterparts.empty()) continue;
    
        GEdge *src = dynamic_cast<GEdge *>(tgt->getMeshMaster());
    
        if(src != NULL && src != tgt) {
          std::map<MVertex *, MVertex *> &v2v = tgt->correspondingVertices;
          std::map<MVertex *, MVertex *> &p2p = tgt->correspondingHOPoints;
          p2p.clear();
    
          Msg::Info("Reconstructing periodicity for curve connection %d - %d",
                    tgt->tag(), src->tag());
    
          std::map<MEdge, MLine *, Less_Edge> srcEdges;
          for(std::size_t i = 0; i < src->getNumMeshElements(); i++) {
            MLine *srcLine = dynamic_cast<MLine *>(src->getMeshElement(i));
            if(!srcLine) {
              Msg::Error("Master element %d is not a line",
                         src->getMeshElement(i)->getNum());
              return;
            }
            srcEdges[MEdge(srcLine->getVertex(0), srcLine->getVertex(1))] = srcLine;
          }
    
          for(std::size_t i = 0; i < tgt->getNumMeshElements(); ++i) {
            MLine *tgtLine = dynamic_cast<MLine *>(tgt->getMeshElement(i));
            MVertex *vtcs[2];
            if(!tgtLine) {
              Msg::Error("Slave element %d is not a line",
                         tgt->getMeshElement(i)->getNum());
              return;
            }
            for(int iVtx = 0; iVtx < 2; iVtx++) {
              MVertex *vtx = tgtLine->getVertex(iVtx);
              std::map<MVertex *, MVertex *>::iterator tIter = v2v.find(vtx);
              if(tIter == v2v.end()) {
                Msg::Error("Cannot find periodic counterpart of node %d"
                           " of curve %d on curve %d", vtx->getNum(), tgt->tag(),
                           src->tag());
                return;
              }
              else
                vtcs[iVtx] = tIter->second;
            }
    
            std::map<MEdge, MLine *, Less_Edge>::iterator srcIter =
              srcEdges.find(MEdge(vtcs[0], vtcs[1]));
            if(srcIter == srcEdges.end()) {
              Msg::Error("Can't find periodic counterpart of mesh edge %d-%d "
                         "on curve %d, connected to mesh edge %d-%d on curve %d",
                         tgtLine->getVertex(0)->getNum(),
                         tgtLine->getVertex(1)->getNum(), tgt->tag(),
                         vtcs[0]->getNum(), vtcs[1]->getNum(), src->tag());
              return;
            }
            else {
              MLine *srcLine = srcIter->second;
              if(tgtLine->getNumVertices() != srcLine->getNumVertices()) throw;
              for(std::size_t i = 2; i < tgtLine->getNumVertices(); i++)
                p2p[tgtLine->getVertex(i)] = srcLine->getVertex(i);
            }
          }
        }
      }
    
      if(CTX::instance()->mesh.hoPeriodic) {
    #if defined(HAVE_OPTHOM)
        std::vector<GEntity *> modelEdges(m->firstEdge(), m->lastEdge());
        HighOrderMeshPeriodicity edgePeriodicity(modelEdges);
        edgePeriodicity.fixPeriodicity();
        edgePeriodicity.fixPeriodicity(); // apply twice for operation order effects
    #else
        Msg::Error("High-order mesh optimization requires the OPTHOM module");
    #endif
      }
    
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        GFace *tgt = *it;
    
        // non complete periodic info (e.g. through extrusion)
        if(tgt->vertexCounterparts.empty()) continue;
    
        GFace *src = dynamic_cast<GFace *>(tgt->getMeshMaster());
        if(src != NULL && src != tgt) {
          Msg::Info("Reconstructing periodicity for surface connection %d - %d",
                    tgt->tag(), src->tag());
    
          std::map<MVertex *, MVertex *> &v2v = tgt->correspondingVertices;
          std::map<MVertex *, MVertex *> &p2p = tgt->correspondingHOPoints;
          p2p.clear();
    
          std::map<MFace, MElement *, Less_Face> srcFaces;
    
          for(std::size_t i = 0; i < src->getNumMeshElements(); ++i) {
            MElement *srcElmt = src->getMeshElement(i);
            int nbVtcs = 0;
            if(dynamic_cast<MTriangle *>(srcElmt)) nbVtcs = 3;
            if(dynamic_cast<MQuadrangle *>(srcElmt)) nbVtcs = 4;
            std::vector<MVertex *> vtcs;
            vtcs.reserve(nbVtcs);
            for(int iVtx = 0; iVtx < nbVtcs; iVtx++) {
              vtcs.push_back(srcElmt->getVertex(iVtx));
            }
            srcFaces[MFace(vtcs)] = srcElmt;
          }
    
          for(std::size_t i = 0; i < tgt->getNumMeshElements(); ++i) {
            MElement *tgtElmt = tgt->getMeshElement(i);
            int nbVtcs = 0;
            if(dynamic_cast<MTriangle *>(tgtElmt)) nbVtcs = 3;
            if(dynamic_cast<MQuadrangle *>(tgtElmt)) nbVtcs = 4;
            std::vector<MVertex *> vtcs;
            for(int iVtx = 0; iVtx < nbVtcs; iVtx++) {
              MVertex *vtx = tgtElmt->getVertex(iVtx);
    
              std::map<MVertex *, MVertex *>::iterator tIter = v2v.find(vtx);
              if(tIter == v2v.end()) {
                Msg::Error("Cannot find periodic counterpart of node %d "
                           "of surface %d on surface %d",
                           vtx->getNum(), tgt->tag(), src->tag());
                return;
              }
              else
                vtcs.push_back(tIter->second);
            }
    
            std::map<MFace, MElement *, Less_Face>::iterator srcIter =
              srcFaces.find(MFace(vtcs));
            if(srcIter == srcFaces.end()) {
              std::ostringstream faceDef;
              for(int iVtx = 0; iVtx < nbVtcs; iVtx++)
                faceDef << vtcs[iVtx]->getNum() << " ";
              Msg::Error("Cannot find periodic counterpart of mesh face %s in "
                         "surface %d on surface %d",
                         faceDef.str().c_str(), tgt->tag(), src->tag());
              return;
            }
            else {
              MElement *srcElmt = srcIter->second;
              for(std::size_t i = nbVtcs; i < srcElmt->getNumVertices(); i++) {
                p2p[tgtElmt->getVertex(i)] = srcElmt->getVertex(i);
              }
            }
          }
        }
      }
    
      if(CTX::instance()->mesh.hoPeriodic) {
    #if defined(HAVE_OPTHOM)
        std::vector<GEntity *> modelFaces;
        modelFaces.insert(modelFaces.end(), m->firstFace(), m->lastFace());
        HighOrderMeshPeriodicity facePeriodicity(modelFaces);
        facePeriodicity.fixPeriodicity();
    #else
        Msg::Error("High-order mesh optimization requires the OPTHOM module");
    #endif
      }
    
      Msg::Debug("Finalized topology of periodic connections");
    }
    
    void SetOrder1(GModel *m, bool onlyVisible)
    {
      m->destroyMeshCaches();
    
      // replace all elements with first order elements
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
        setFirstOrder(*it, (*it)->lines, onlyVisible);
      }
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        setFirstOrder(*it, (*it)->triangles, onlyVisible);
        setFirstOrder(*it, (*it)->quadrangles, onlyVisible);
      }
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
        setFirstOrder(*it, (*it)->tetrahedra, onlyVisible);
        setFirstOrder(*it, (*it)->hexahedra, onlyVisible);
        setFirstOrder(*it, (*it)->prisms, onlyVisible);
        setFirstOrder(*it, (*it)->pyramids, onlyVisible);
      }
    
      // remove all high order vertices
      std::vector<MVertex *> newHOVert;
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
        updateHighOrderVertices(*it, newHOVert, onlyVisible);
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
        updateHighOrderVertices(*it, newHOVert, onlyVisible);
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
        updateHighOrderVertices(*it, newHOVert, onlyVisible);
    
      updatePeriodicEdgesAndFaces(m);
    }
    
    void checkHighOrderTriangles(const char *cc, GModel *m,
                                 std::vector<MElement *> &bad, double &minJGlob)
    {
      bad.clear();
      minJGlob = 1.0;
      double minGGlob = 1.0;
      double avg = 0.0;
      int count = 0, nbfair = 0;
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        for(std::size_t i = 0; i < (*it)->triangles.size(); i++) {
          MTriangle *t = (*it)->triangles[i];
          double disto_ = t->distoShapeMeasure();
          double gamma_ = t->gammaShapeMeasure();
          double disto = disto_;
          minJGlob = std::min(minJGlob, disto);
          minGGlob = std::min(minGGlob, gamma_);
          avg += disto;
          count++;
          if(disto < 0)
            bad.push_back(t);
          else if(disto < 0.2)
            nbfair++;
        }
        /*
        for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++){
          MQuadrangle *t = (*it)->quadrangles[i];
          double disto_ = t->distoShapeMeasure();
          double gamma_ = t->gammaShapeMeasure();
          double disto = disto_;
          minJGlob = std::min(minJGlob, disto);
          minGGlob = std::min(minGGlob, gamma_);
          avg += disto; count++;
          if (disto < 0) bad.push_back(t);
          else if (disto < 0.2) nbfair++;
        }
        */
      }
      if(!count) return;
      if(minJGlob > 0)
        Msg::Info(
          "%s: worst distortion = %g (%d elements in ]0, 0.2]); worst gamma = %g",
          cc, minJGlob, nbfair, minGGlob);
      else
        Msg::Warning(
          "%s: worst distortion = %g (avg = %g, %d elements with jac. < 0); "
          "worst gamma = %g",
          cc, minJGlob, avg / (count ? count : 1), bad.size(), minGGlob);
    }
    
    void checkHighOrderTetrahedron(const char *cc, GModel *m,
                                   std::vector<MElement *> &bad, double &minJGlob)
    {
      bad.clear();
      minJGlob = 1.0;
      double avg = 0.0;
      int count = 0, nbfair = 0;
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
        for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++) {
          MTetrahedron *t = (*it)->tetrahedra[i];
          double disto_ = t->distoShapeMeasure();
          minJGlob = std::min(minJGlob, disto_);
          avg += disto_;
          count++;
          if(disto_ < 0)
            bad.push_back(t);
          else if(disto_ < 0.2)
            nbfair++;
        }
      }
      if(!count) return;
    
      if(minJGlob > 0)
        Msg::Info("%s: worst distortion = %g (%d elements in ]0, 0.2])", cc,
                  minJGlob, nbfair);
      else
        Msg::Warning(
          "%s: worst distortion = %g (avg = %g, %d elements with jac. < 0)", cc,
          minJGlob, avg / (count ? count : 1), bad.size());
    }
    
    void getMeshInfoForHighOrder(GModel *gm, int &meshOrder, bool &complete,
                                 bool &CAD)
    {
      meshOrder = -1;
      CAD = true;
      complete = true;
      // The following code does not allow to determine accurately if elements
      // are complete or not. To be more precise, we should try to find a hexahedron
      // if order = 2+, a prism or a pyramid if order = 3+ or a tet if order = 4+
      // and so on...
      // But it is more likely that we want complete elements so always
      // setting true to 'complete' variable is acceptable.
      for(GModel::riter itr = gm->firstRegion(); itr != gm->lastRegion(); ++itr) {
        if((*itr)->getNumMeshElements()) {
          meshOrder = (*itr)->getMeshElement(0)->getPolynomialOrder();
          // complete = (meshOrder <= 2) ? 1 :
          //   (*itr)->getMeshElement(0)->getNumVolumeVertices();
          break;
        }
      }
      for(GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
        if((*itf)->getNumMeshElements()) {
          if(meshOrder == -1) {
            meshOrder = (*itf)->getMeshElement(0)->getPolynomialOrder();
            // complete = (meshOrder <= 2) ? 1 :
            //   (*itf)->getMeshElement(0)->getNumFaceVertices();
            // if((*itf)->geomType() == GEntity::DiscreteSurface) CAD = false;
            break;
          }
        }
      }
    }
    
    void SetOrderN(GModel *m, int order, bool linear, bool incomplete,
                   bool onlyVisible)
    {
      // replace all the elements in the mesh with second order elements
      // by creating unique vertices on the edges/faces of the mesh:
      //
      // - if linear is set to true, new vertices are created by linear
      //   interpolation between existing ones. If not, new vertices are
      //   created on the exact geometry, provided that the geometrical
      //   edges/faces are discretized with 1D/2D elements. (I.e., if
      //   there are only 3D elements in the mesh then any new vertices on
      //   the boundary will always be created by linear interpolation,
      //   whether linear is set or not.)
      //
      // - if incomplete is set to true, we only create new vertices on
      //   edges (creating 8-node quads, 20-node hexas, etc., instead of
      //   9-node quads, 27-node hexas, etc.)
    
      // - if onlyVisible is true, then only the visible entities will be curved.
    
      int nPts = order - 1;
    
      char msg[256];
      sprintf(msg, "Meshing order %d (curvilinear %s)...", order,
              linear ? "off" : "on");
    
      Msg::StatusBar(true, msg);
    
      double t1 = Cpu();
    
      m->destroyMeshCaches();
    
      // Keep track of vertex/entities created
      edgeContainer edgeVertices;
      faceContainer faceVertices;
      std::map<GEntity *, std::vector<MVertex *> > newHOVert;
    
      Msg::ResetProgressMeter();
    
      int counter = 0,
          nTot = m->getNumEdges() + m->getNumFaces() + m->getNumRegions();
    
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
        Msg::Info("Meshing curve %d order %d", (*it)->tag(), order);
        Msg::ProgressMeter(++counter, nTot, false, msg);
        if(onlyVisible && !(*it)->getVisibility()) continue;
        setHighOrder(*it, newHOVert[*it], edgeVertices, linear, nPts);
      }
    
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        Msg::Info("Meshing surface %d order %d", (*it)->tag(), order);
        Msg::ProgressMeter(++counter, nTot, false, msg);
        if(onlyVisible && !(*it)->getVisibility()) continue;
        setHighOrder(*it, newHOVert[*it], edgeVertices, faceVertices, linear,
                     incomplete, nPts);
        if((*it)->getColumns() != 0) (*it)->getColumns()->clearElementData();
      }
    
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
        Msg::Info("Meshing volume %d order %d", (*it)->tag(), order);
        Msg::ProgressMeter(++counter, nTot, false, msg);
        if(onlyVisible && !(*it)->getVisibility()) continue;
        setHighOrder(*it, newHOVert[*it], edgeVertices, faceVertices, incomplete,
                     nPts);
        if((*it)->getColumns() != 0) (*it)->getColumns()->clearElementData();
      }
    
      // Update all high order vertices
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
        updateHighOrderVertices(*it, newHOVert[*it], onlyVisible);
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
        updateHighOrderVertices(*it, newHOVert[*it], onlyVisible);
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
        updateHighOrderVertices(*it, newHOVert[*it], onlyVisible);
    
      if(CTX::instance()->mesh.hoOptimize) {
    #if defined(HAVE_OPTHOM)
        // Determine mesh dimension and curve BL elements
        FastCurvingParameters p;
        p.dim = 0;
        p.curveOuterBL = FastCurvingParameters::OUTER_CURVE;
        p.thickness = false;
        // p.optimizeGeometry = true;
        for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
          if((*it)->getNumMeshElements() > 0) {
            p.dim = 3;
            break;
          }
        if(p.dim == 0)
          for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
            if((*it)->getNumMeshElements() > 0) {
              p.dim = 2;
              break;
            }
        if(p.dim == 2) HighOrderMeshFastCurving(GModel::current(), p, true);
    #else
      // Msg::Error("High-order mesh optimization requires the OPTHOM module");
    #endif
      }
    
      updatePeriodicEdgesAndFaces(m);
    
      double t2 = Cpu();
    
      std::vector<MElement *> bad;
      double worst;
      checkHighOrderTriangles("Surface mesh", m, bad, worst);
      checkHighOrderTetrahedron("Volume mesh", m, bad, worst);
      // FIXME : add other element check
    
      Msg::StatusBar(true, "Done meshing order %d (%g s)", order, t2 - t1);
    }
    
    void SetHighOrderComplete(GModel *m, bool onlyVisible)
    {
      faceContainer faceVertices;
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        if(onlyVisible && !(*it)->getVisibility()) continue;
        std::vector<MVertex *> dumNewHOVert;
        std::vector<MTriangle *> newT;
        for(std::size_t i = 0; i < (*it)->triangles.size(); i++) {
          MTriangle *t = (*it)->triangles[i];
          std::vector<MVertex *> vv;
          for(std::size_t j = 3; j < t->getNumVertices() - t->getNumFaceVertices();
              j++)
            vv.push_back(t->getVertex(j));
          int nPts = t->getPolynomialOrder() - 1;
          getFaceVertices(*it, t, vv, dumNewHOVert, faceVertices, false, nPts);
          newT.push_back(new MTriangleN(t->getVertex(0), t->getVertex(1),
                                        t->getVertex(2), vv, nPts + 1, 0,
                                        t->getPartition()));
          delete t;
        }
        (*it)->triangles = newT;
    
        std::vector<MQuadrangle *> newQ;
        for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++) {
          MQuadrangle *t = (*it)->quadrangles[i];
    
          std::vector<MVertex *> vv;
          vv.reserve(t->getNumVertices() - t->getNumFaceVertices() - 4);
    
          for(std::size_t j = 4; j < t->getNumVertices() - t->getNumFaceVertices();
              j++)
            vv.push_back(t->getVertex(j));
    
          int nPts = t->getPolynomialOrder() - 1;
          getFaceVertices(*it, t, vv, dumNewHOVert, faceVertices, false, nPts);
          newQ.push_back(new MQuadrangleN(t->getVertex(0), t->getVertex(1),
                                          t->getVertex(2), t->getVertex(3), vv,
                                          nPts + 1, 0, t->getPartition()));
          delete t;
        }
        (*it)->quadrangles = newQ;
    
        std::set<MVertex *> newV;
        for(std::size_t i = 0; i < (*it)->getNumMeshElements(); ++i) {
          MElement *e = (*it)->getMeshElement(i);
          for(std::size_t j = 0; j < e->getNumVertices(); j++)
            newV.insert(e->getVertex(j));
        }
        (*it)->mesh_vertices.clear();
        (*it)->mesh_vertices.insert((*it)->mesh_vertices.begin(), newV.begin(),
                                    newV.end());
      }
    
      updatePeriodicEdgesAndFaces(m);
    }
    
    void SetHighOrderIncomplete(GModel *m, bool onlyVisible)
    {
      std::set<MVertex *> toDelete;
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        if(onlyVisible && !(*it)->getVisibility()) continue;
        std::vector<MTriangle *> newT;
    
        for(std::size_t i = 0; i < (*it)->triangles.size(); i++) {
          MTriangle *t = (*it)->triangles[i];
          std::vector<MVertex *> vt;
          int order = t->getPolynomialOrder();
          for(std::size_t j = 3; j < t->getNumVertices() - t->getNumFaceVertices();
              j++)
            vt.push_back(t->getVertex(j));
          for(std::size_t j = t->getNumVertices() - t->getNumFaceVertices();
              j < t->getNumVertices(); j++)
            toDelete.insert(t->getVertex(j));
          newT.push_back(new MTriangleN(t->getVertex(0), t->getVertex(1),
                                        t->getVertex(2), vt, order, 0,
                                        t->getPartition()));
          delete t;
        }
        (*it)->triangles = newT;
    
        std::vector<MQuadrangle *> newQ;
        for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++) {
          MQuadrangle *q = (*it)->quadrangles[i];
          std::vector<MVertex *> vt;
          int nPts = q->getPolynomialOrder() - 1;
          for(std::size_t j = 4; j < q->getNumVertices() - q->getNumFaceVertices();
              j++)
            vt.push_back(q->getVertex(j));
          for(std::size_t j = q->getNumVertices() - q->getNumFaceVertices();
              j < q->getNumVertices(); j++)
            toDelete.insert(q->getVertex(j));
          newQ.push_back(new MQuadrangleN(q->getVertex(0), q->getVertex(1),
                                          q->getVertex(2), q->getVertex(3), vt,
                                          nPts + 1, 0, q->getPartition()));
          delete q;
        }
        (*it)->quadrangles = newQ;
    
        std::vector<MVertex *> newV;
        int numd = 0;
        for(std::size_t i = 0; i < (*it)->mesh_vertices.size(); ++i) {
          if(toDelete.find((*it)->mesh_vertices[i]) == toDelete.end())
            newV.push_back((*it)->mesh_vertices[i]);
          else {
            delete(*it)->mesh_vertices[i];
            numd++;
          }
        }
        (*it)->mesh_vertices = newV;
      }
    
      updatePeriodicEdgesAndFaces(m);
    }