Skip to content
Snippets Groups Projects
Select Git revision
  • f577eca09fb7f8171bebd192407e0f0ece4cb473
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

Mesh.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    meshGFaceOptimize.cpp 46.20 KiB
    // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file in the Gmsh root directory for license information.
    // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
    
    #include <stack>
    #include "GmshConfig.h"
    #include "meshGFaceOptimize.h"
    #include "qualityMeasures.h"
    #include "GFace.h"
    #include "GEdge.h"
    #include "GVertex.h"
    #include "GModel.h"
    #include "MVertex.h"
    #include "MTriangle.h"
    #include "MQuadrangle.h"
    #include "MLine.h"
    #include "BackgroundMeshTools.h"
    #include "Numeric.h"
    #include "GmshMessage.h"
    #include "Context.h"
    #include "OS.h"
    #include "SVector3.h"
    #include "SPoint3.h"
    #include "meshRelocateVertex.h"
    #include "Field.h"
    
    #if defined(HAVE_BLOSSOM)
    extern "C" struct CCdatagroup;
    extern "C" int perfect_match(int ncount, CCdatagroup *dat, int ecount,
                                 int **elist, int **elen, char *blo_filename,
                                 char *mat_filename, int just_fractional,
                                 int no_fractional, int use_all_trees,
                                 int partialprice, double *totalzeit);
    #endif
    
    RecombineTriangle::RecombineTriangle(const MEdge &me, MElement *_t1,
                                         MElement *_t2, Field *cross_field)
      : t1(_t1), t2(_t2)
    {
      n1 = me.getVertex(0);
      n2 = me.getVertex(1);
      n3 = 0;
      n4 = 0;
    
      if(t1->getVertex(0) != n1 && t1->getVertex(0) != n2)
        n3 = t1->getVertex(0);
      else if(t1->getVertex(1) != n1 && t1->getVertex(1) != n2)
        n3 = t1->getVertex(1);
      else if(t1->getVertex(2) != n1 && t1->getVertex(2) != n2)
        n3 = t1->getVertex(2);
    
      if(t2->getVertex(0) != n1 && t2->getVertex(0) != n2)
        n4 = t2->getVertex(0);
      else if(t2->getVertex(1) != n1 && t2->getVertex(1) != n2)
        n4 = t2->getVertex(1);
      else if(t2->getVertex(2) != n1 && t2->getVertex(2) != n2)
        n4 = t2->getVertex(2);
    
      if(!n3) {
        Msg::Warning("Invalid quadrangle in recombination");
        n3 = n1;
      }
      if(!n4) {
        Msg::Warning("Invalid quadrangle in recombination");
        n4 = n2;
      }
    
      MQuadrangle q(n1, n3, n2, n4);
      angle = q.etaShapeMeasure();
    
      double const a1 = 180 * angle3Vertices(n1, n4, n2) / M_PI;
      double const a2 = 180 * angle3Vertices(n4, n2, n3) / M_PI;
      double const a3 = 180 * angle3Vertices(n2, n3, n1) / M_PI;
      double const a4 = 180 * angle3Vertices(n3, n1, n4) / M_PI;
      quality = std::abs(90. - a1);
      quality = std::max(std::abs(90. - a2), quality);
      quality = std::max(std::abs(90. - a3), quality);
      quality = std::max(std::abs(90. - a4), quality);
    
      if(cross_field) {
        SVector3 t1;
        SVector3 t2(me.getVertex(1)->x() - me.getVertex(0)->x(),
                    me.getVertex(1)->y() - me.getVertex(0)->y(),
                    me.getVertex(1)->z() - me.getVertex(0)->z());
        (*cross_field)(0.5 * (me.getVertex(0)->x() + me.getVertex(1)->x()),
                       0.5 * (me.getVertex(0)->y() + me.getVertex(1)->y()),
                       0.5 * (me.getVertex(0)->z() + me.getVertex(1)->z()), t1,
                       NULL);
        //    double L1 = t1.norm();
        t1.normalize();
        //    double L2 = t2.norm();
        t2.normalize();
        SVector3 n = crossprod(t1, t2);
        double ErrDir = std::min(fabs(dot(t1, t2)), n.norm());
        //    double ErrL = fabs(L1-L2)/(L1);
        //    printf("quality %g err %g L %g %g\n",quality,ErrL,L1,L2);
        quality /= ErrDir;
        if(n1->onWhat()->dim() < 2 || n2->onWhat()->dim() < 2) quality /= 10;
      }
    }
    
    edge_angle::edge_angle(MVertex *_v1, MVertex *_v2, MElement *t1, MElement *t2)
      : v1(_v1), v2(_v2)
    {
      if(!t2)
        angle = 0;
      else {
        double c1[3];
        double c2[3];
        double c3[3];
        {
          MVertex *p1 = t1->getVertex(0);
          MVertex *p2 = t1->getVertex(1);
          MVertex *p3 = t1->getVertex(2);
          double a[3] = {p1->x() - p2->x(), p1->y() - p2->y(), p1->z() - p2->z()};
          double b[3] = {p1->x() - p3->x(), p1->y() - p3->y(), p1->z() - p3->z()};
          c1[2] = a[0] * b[1] - a[1] * b[0];
          c1[1] = -a[0] * b[2] + a[2] * b[0];
          c1[0] = a[1] * b[2] - a[2] * b[1];
        }
        {
          MVertex *p1 = t2->getVertex(0);
          MVertex *p2 = t2->getVertex(1);
          MVertex *p3 = t2->getVertex(2);
          double a[3] = {p1->x() - p2->x(), p1->y() - p2->y(), p1->z() - p2->z()};
          double b[3] = {p1->x() - p3->x(), p1->y() - p3->y(), p1->z() - p3->z()};
          c2[2] = a[0] * b[1] - a[1] * b[0];
          c2[1] = -a[0] * b[2] + a[2] * b[0];
          c2[0] = a[1] * b[2] - a[2] * b[1];
        }
        norme(c1);
        norme(c2);
        prodve(c1, c2, c3);
        double const cosa = prosca(c1, c2);
        double const sina = norme(c3);
        angle = std::atan2(sina, cosa);
      }
    }
    
    static void setLcsInit(MTriangle *t, std::map<MVertex *, double> &vSizes)
    {
      for(int i = 0; i < 3; i++) {
        for(int j = i + 1; j < 3; j++) {
          MVertex *vi = t->getVertex(i);
          MVertex *vj = t->getVertex(j);
          vSizes[vi] = -1;
          vSizes[vj] = -1;
        }
      }
    }
    
    static void setLcs(MTriangle *t, std::map<MVertex *, double> &vSizes,
                       bidimMeshData &data)
    {
      for(int i = 0; i < 3; i++) {
        for(int j = i + 1; j < 3; j++) {
          MVertex *vi = t->getVertex(i);
          MVertex *vj = t->getVertex(j);
          if(vi != data.equivalent(vj) && vj != data.equivalent(vi)) {
            double dx = vi->x() - vj->x();
            double dy = vi->y() - vj->y();
            double dz = vi->z() - vj->z();
            double l = sqrt(dx * dx + dy * dy + dz * dz);
            auto iti = vSizes.find(vi);
            auto itj = vSizes.find(vj);
            if(iti->second < 0 || iti->second > l) iti->second = l;
            if(itj->second < 0 || itj->second > l) itj->second = l;
          }
        }
      }
    }
    
    bool buildMeshGenerationDataStructures(
      GFace *gf, std::set<MTri3 *, compareTri3Ptr> &AllTris, bidimMeshData &data)
    {
      std::map<MVertex *, double> vSizesMap;
    
      for(std::size_t i = 0; i < gf->triangles.size(); i++)
        setLcsInit(gf->triangles[i], vSizesMap);
    
      auto itfind = vSizesMap.find(nullptr);
      if(itfind != vSizesMap.end()) {
        Msg::Error("Some NULL points exist in 2D mesh");
        return false;
      }
    
      for(std::size_t i = 0; i < gf->triangles.size(); i++)
        setLcs(gf->triangles[i], vSizesMap, data);
    
      // take care of embedded vertices
      std::set<MVertex *> embeddedVertices;
      {
        std::vector<GVertex *> emb_vertx = gf->getEmbeddedVertices();
        auto itvx = emb_vertx.begin();
        while(itvx != emb_vertx.end()) {
          if((*itvx)->mesh_vertices.size()) {
            MVertex *v = *((*itvx)->mesh_vertices.begin());
            vSizesMap[v] =
              std::min(vSizesMap[v], (*itvx)->prescribedMeshSizeAtVertex());
            embeddedVertices.insert(v);
          }
          ++itvx;
        }
      }
    
      // take good care of embedded edges
      {
        std::vector<GEdge *> embedded_edges = gf->getEmbeddedEdges();
        auto ite = embedded_edges.begin();
        while(ite != embedded_edges.end()) {
          if(!(*ite)->isMeshDegenerated()) {
            for(std::size_t i = 0; i < (*ite)->lines.size(); i++)
              data.internalEdges.insert(MEdge((*ite)->lines[i]->getVertex(0),
                                              (*ite)->lines[i]->getVertex(1)));
          }
          ++ite;
        }
      }
    
      // take care of small edges in  order not to "pollute" the size field
      {
        std::vector<GEdge *> _edges = gf->edges();
        auto ite = _edges.begin();
        while(ite != _edges.end()) {
          if(!(*ite)->isMeshDegenerated()) {
            for(std::size_t i = 0; i < (*ite)->lines.size(); i++) {
              double d = distance((*ite)->lines[i]->getVertex(0),
                                  (*ite)->lines[i]->getVertex(1));
              double d0 = vSizesMap[(*ite)->lines[i]->getVertex(0)];
              double d1 = vSizesMap[(*ite)->lines[i]->getVertex(1)];
              if(d0 < .5 * d) vSizesMap[(*ite)->lines[i]->getVertex(0)] = .5 * d;
              if(d1 < .5 * d) vSizesMap[(*ite)->lines[i]->getVertex(1)] = .5 * d;
            }
          }
          ++ite;
        }
      }
    
      for(auto it = vSizesMap.begin(); it != vSizesMap.end(); ++it) {
        SPoint2 param;
        reparamMeshVertexOnFace(it->first, gf, param);
        // Add size of background mesh to embedded vertices. For the other nodes,
        // use the size in vSizesMap
        const double lcBGM = (embeddedVertices.count(it->first) > 0) ?
                               BGM_MeshSize(gf, param[0], param[1], it->first->x(),
                                            it->first->y(), it->first->z()) :
                               it->second;
        data.addVertex(it->first, param[0], param[1], it->second, lcBGM);
      }
      for(std::size_t i = 0; i < gf->triangles.size(); i++) {
        double lc = 0.3333333333 *
                    (data.vSizes[data.getIndex(gf->triangles[i]->getVertex(0))] +
                     data.vSizes[data.getIndex(gf->triangles[i]->getVertex(1))] +
                     data.vSizes[data.getIndex(gf->triangles[i]->getVertex(2))]);
        double lcBGM =
          0.3333333333 *
          (data.vSizesBGM[data.getIndex(gf->triangles[i]->getVertex(0))] +
           data.vSizesBGM[data.getIndex(gf->triangles[i]->getVertex(1))] +
           data.vSizesBGM[data.getIndex(gf->triangles[i]->getVertex(2))]);
    
        double LL = Extend1dMeshIn2dSurfaces(gf) ? std::min(lc, lcBGM) : lcBGM;
        AllTris.insert(new MTri3(gf->triangles[i], LL, nullptr, &data, gf));
      }
      gf->triangles.clear();
      connectTriangles(AllTris);
    
      return true;
    }
    
    void computeEquivalences(GFace *gf, bidimMeshData &data)
    {
      if(data.equivalence) {
        std::vector<MTriangle *> newT;
        for(std::size_t i = 0; i < gf->triangles.size(); i++) {
          MTriangle *t = gf->triangles[i];
          MVertex *v[3];
          for(int j = 0; j < 3; j++) {
            v[j] = t->getVertex(j);
            auto it = data.equivalence->find(v[j]);
            if(it != data.equivalence->end()) { v[j] = it->second; }
          }
          if(v[0] != v[1] && v[0] != v[2] && v[2] != v[1])
            newT.push_back(new MTriangle(v[0], v[1], v[2]));
          delete t;
        }
        gf->triangles = newT;
      }
    }
    
    struct equivalentTriangle {
      MTriangle *_t;
      MVertex *_v[3];
      equivalentTriangle(MTriangle *t, std::map<MVertex *, MVertex *> *equivalence)
        : _t(t)
      {
        for(int i = 0; i < 3; i++) {
          MVertex *v = t->getVertex(i);
          auto it = equivalence->find(v);
          if(it == equivalence->end())
            _v[i] = v;
          else
            _v[i] = it->second;
        }
        std::sort(_v, _v + 3);
      }
      bool operator<(const equivalentTriangle &other) const
      {
        for(int i = 0; i < 3; i++) {
          if(other._v[i] > _v[i]) return true;
          if(other._v[i] < _v[i]) return false;
        }
        return false;
      }
    };
    
    bool computeEquivalentTriangles(GFace *gf,
                                    std::map<MVertex *, MVertex *> *equivalence)
    {
      if(!equivalence) return false;
      std::vector<MTriangle *> WTF;
      if(!equivalence) return false;
      std::set<equivalentTriangle> eqTs;
      for(std::size_t i = 0; i < gf->triangles.size(); i++) {
        equivalentTriangle et(gf->triangles[i], equivalence);
        auto iteq = eqTs.find(et);
        if(iteq == eqTs.end())
          eqTs.insert(et);
        else {
          WTF.push_back(iteq->_t);
          WTF.push_back(gf->triangles[i]);
        }
      }
    
      if(WTF.size()) {
        Msg::Info("%d triangles are equivalent", WTF.size());
        for(std::size_t i = 0; i < WTF.size(); i++) {}
        return true;
      }
      return false;
    }
    
    void splitEquivalentTriangles(GFace *gf, bidimMeshData &data)
    {
      computeEquivalentTriangles(gf, data.equivalence);
    }
    
    void transferDataStructure(GFace *gf,
                               std::set<MTri3 *, compareTri3Ptr> &AllTris,
                               bidimMeshData &data)
    {
      while(1) {
        if(AllTris.begin() == AllTris.end()) break;
        MTri3 *worst = *AllTris.begin();
        if(worst->isDeleted())
          delete worst->tri();
        else
          gf->triangles.push_back(worst->tri());
        delete worst;
        AllTris.erase(AllTris.begin());
      }
    
      // make sure all the triangles are oriented in the same way in
      // parameter space (it would be nicer to change the actual algorithm
      // to ensure that we create correctly-oriented triangles in the
      // first place)
    
      // if BL triangles are considered, then all that is WRONG !
    
      if(gf->triangles.size() > 1) {
        bool BL = !gf->getColumns()->_toFirst.empty();
    
        double n1[3], n2[3];
        MTriangle *t = gf->triangles[0];
        MVertex *v0 = t->getVertex(0), *v1 = t->getVertex(1), *v2 = t->getVertex(2);
    
        if(!BL) {
          int index0 = data.getIndex(v0);
          int index1 = data.getIndex(v1);
          int index2 = data.getIndex(v2);
          normal3points(data.Us[index0], data.Vs[index0], 0., data.Us[index1],
                        data.Vs[index1], 0., data.Us[index2], data.Vs[index2], 0.,
                        n1);
        }
        else {
          // BL --> PLANAR FACES !!!
          normal3points(v0->x(), v0->y(), v0->z(), v1->x(), v1->y(), v1->z(),
                        v2->x(), v2->y(), v2->z(), n1);
        }
        for(std::size_t j = 1; j < gf->triangles.size(); j++) {
          t = gf->triangles[j];
          v0 = t->getVertex(0);
          v1 = t->getVertex(1);
          v2 = t->getVertex(2);
          if(!BL) {
            int index0 = data.getIndex(v0);
            int index1 = data.getIndex(v1);
            int index2 = data.getIndex(v2);
            normal3points(data.Us[index0], data.Vs[index0], 0., data.Us[index1],
                          data.Vs[index1], 0., data.Us[index2], data.Vs[index2], 0.,
                          n2);
          }
          else {
            // BL --> PLANAR FACES !!!
            normal3points(v0->x(), v0->y(), v0->z(), v1->x(), v1->y(), v1->z(),
                          v2->x(), v2->y(), v2->z(), n2);
          }
          // orient the bignou
          if(prosca(n1, n2) < 0.0) t->reverse();
        }
      }
      splitEquivalentTriangles(gf, data);
      computeEquivalences(gf, data);
    }
    
    template <class T>
    void buildEdgeToElement(std::vector<T *> &elements, e2t_cont &adj)
    {
      for(std::size_t i = 0; i < elements.size(); i++) {
        T *t = elements[i];
        for(int j = 0; j < t->getNumEdges(); j++) {
          MEdge e = t->getEdge(j);
          auto it = adj.find(e);
          if(it == adj.end()) {
            std::pair<MElement *, MElement *> one =
              std::make_pair(t, (MElement *)nullptr);
            adj[e] = one;
          }
          else {
            it->second.second = t;
          }
        }
      }
    }
    
    void buildEdgeToElement(GFace *gf, e2t_cont &adj)
    {
      adj.clear();
      buildEdgeToElement(gf->triangles, adj);
      buildEdgeToElement(gf->quadrangles, adj);
    }
    
    void buildEdgeToTriangle(std::vector<MTriangle *> &tris, e2t_cont &adj)
    {
      adj.clear();
      buildEdgeToElement(tris, adj);
    }
    
    void buildEdgeToElements(std::vector<MElement *> &tris, e2t_cont &adj)
    {
      adj.clear();
      buildEdgeToElement(tris, adj);
    }
    
    void buildListOfEdgeAngle(e2t_cont adj, std::vector<edge_angle> &edges_detected,
                              std::vector<edge_angle> &edges_lonly)
    {
      auto it = adj.begin();
      for(; it != adj.end(); ++it) {
        if(it->second.second)
          edges_detected.push_back(edge_angle(it->first.getVertex(0),
                                              it->first.getVertex(1),
                                              it->second.first, it->second.second));
        else
          edges_lonly.push_back(edge_angle(it->first.getVertex(0),
                                           it->first.getVertex(1), it->second.first,
                                           it->second.second));
      }
      std::sort(edges_detected.begin(), edges_detected.end());
    }
    
    static void parametricCoordinates(MElement *t, GFace *gf, double u[4],
                                      double v[4], MVertex *close = nullptr)
    {
      for(std::size_t j = 0; j < t->getNumVertices(); j++) {
        MVertex *ver = t->getVertex(j);
        SPoint2 param, dummy;
        if(!close)
          reparamMeshVertexOnFace(ver, gf, param);
        else
          reparamMeshEdgeOnFace(ver, close, gf, param, dummy);
        u[j] = param[0];
        v[j] = param[1];
      }
    }
    
    double surfaceFaceUV(MElement *t, GFace *gf, bool maximal = true)
    {
      const int N = t->getNumVertices();
      if(N > 4) {
        Msg::Warning("surfaceFaceUV only for first order elements");
        return 0;
      }
    
      double u[4], v[4];
      parametricCoordinates(t, gf, u, v);
      if(N == 3)
        return 0.5 *
               fabs((u[1] - u[0]) * (v[2] - v[0]) - (u[2] - u[0]) * (v[1] - v[0]));
      else {
        const double a1 =
          0.5 *
            fabs((u[1] - u[0]) * (v[2] - v[0]) - (u[2] - u[0]) * (v[1] - v[0])) +
          0.5 * fabs((u[3] - u[2]) * (v[0] - v[2]) - (u[0] - u[2]) * (v[3] - v[2]));
        const double a2 =
          0.5 *
            fabs((u[2] - u[1]) * (v[3] - v[1]) - (u[3] - u[1]) * (v[2] - v[1])) +
          0.5 * fabs((u[0] - u[3]) * (v[1] - v[3]) - (u[1] - u[3]) * (v[0] - v[3]));
        return maximal ? std::max(a2, a1) : std::min(a2, a1);
      }
    }
    
    static int _removeTwoQuadsNodes(GFace *gf)
    {
      v2t_cont adj;
      buildVertexToElement(gf->triangles, adj);
      buildVertexToElement(gf->quadrangles, adj);
      auto it = adj.begin();
      std::set<MElement *> touched;
      std::set<MVertex *> vtouched;
      while(it != adj.end()) {
        MVertex *v = it->first;
        if(it->second.size() == 2 && v->onWhat() == gf) {
          MElement *q1 = it->second[0];
          MElement *q2 = it->second[1];
          if(q1->getNumVertices() == 4 && q2->getNumVertices() == 4 &&
             touched.find(q1) == touched.end() &&
             touched.find(q2) == touched.end()) {
            int comm = 0;
            for(int i = 0; i < 4; i++) {
              if(q1->getVertex(i) == v) {
                comm = i;
                break;
              }
            }
            MVertex *v1 = q1->getVertex((comm + 1) % 4);
            MVertex *v2 = q1->getVertex((comm + 2) % 4);
            MVertex *v3 = q1->getVertex((comm + 3) % 4);
            MVertex *v4 = nullptr;
            for(int i = 0; i < 4; i++) {
              if(q2->getVertex(i) != v1 && q2->getVertex(i) != v3 &&
                 q2->getVertex(i) != v) {
                v4 = q2->getVertex(i);
                break;
              }
            }
            if(!v4) {
              Msg::Error("Bug in removeTwoQuadsNodes %p,%p,%p", v1, v2, v3);
              q1->writePOS(stdout, true, false, false, false, false, false);
              q2->writePOS(stdout, true, false, false, false, false, false);
              return 0;
            }
            MQuadrangle *q = new MQuadrangle(v1, v2, v3, v4);
            touched.insert(q1);
            touched.insert(q2);
            gf->quadrangles.push_back(q);
            vtouched.insert(v);
          }
        }
        it++;
      }
      std::vector<MQuadrangle *> quadrangles2;
      quadrangles2.reserve(gf->quadrangles.size() - touched.size());
      for(std::size_t i = 0; i < gf->quadrangles.size(); i++) {
        if(touched.find(gf->quadrangles[i]) == touched.end()) {
          quadrangles2.push_back(gf->quadrangles[i]);
        }
        else {
          delete gf->quadrangles[i];
        }
      }
      gf->quadrangles = quadrangles2;
    
      std::vector<MVertex *> mesh_vertices2;
      mesh_vertices2.reserve(gf->mesh_vertices.size() - vtouched.size());
      for(std::size_t i = 0; i < gf->mesh_vertices.size(); i++) {
        if(vtouched.find(gf->mesh_vertices[i]) == vtouched.end()) {
          mesh_vertices2.push_back(gf->mesh_vertices[i]);
        }
        else {
          delete gf->mesh_vertices[i];
        }
      }
      gf->mesh_vertices = mesh_vertices2;
    
      return vtouched.size();
    }
    
    int removeTwoQuadsNodes(GFace *gf)
    {
      int nbRemove = 0;
      while(1) {
        int x = _removeTwoQuadsNodes(gf);
        if(!x) break;
        nbRemove += x;
      }
      Msg::Debug("%i two-quadrangles nodes removed", nbRemove);
      return nbRemove;
    }
    
    static bool _tryToCollapseThatVertex2(GFace *gf, std::vector<MElement *> &e1,
                                          std::vector<MElement *> &e2, MElement *q,
                                          MVertex *v1, MVertex *v2)
    {
      std::vector<MElement *> e = e1;
      e.insert(e.end(), e2.begin(), e2.end());
    
      double x1 = v1->x();
      double y1 = v1->y();
      double z1 = v1->z();
    
      double x2 = v2->x();
      double y2 = v2->y();
      double z2 = v2->z();
    
      // new position of v1 && v2
      double initialGuess[2] = {0, 0};
      GPoint pp = gf->closestPoint(
        SPoint3(0.5 * (x1 + x2), 0.5 * (y1 + y2), 0.5 * (z1 + z2)), initialGuess);
    
      double worst_quality_old = 1.0;
      double worst_quality_new = 1.0;
    
      int count = 0;
      for(std::size_t j = 0; j < e.size(); ++j) {
        if(e[j] != q) {
          count++;
          worst_quality_old = std::min(worst_quality_old, e[j]->etaShapeMeasure());
          v1->x() = pp.x();
          v1->y() = pp.y();
          v1->z() = pp.z();
          v2->x() = pp.x();
          v2->y() = pp.y();
          v2->z() = pp.z();
          worst_quality_new = std::min(worst_quality_new, e[j]->etaShapeMeasure());
          v1->x() = x1;
          v1->y() = y1;
          v1->z() = z1;
          v2->x() = x2;
          v2->y() = y2;
          v2->z() = z2;
        }
      }
    
      if(worst_quality_new > worst_quality_old) {
        v1->x() = pp.x();
        v1->y() = pp.y();
        v1->z() = pp.z();
        for(std::size_t j = 0; j < e.size(); ++j) {
          if(e[j] != q) {
            for(int k = 0; k < 4; k++) {
              if(e[j]->getVertex(k) == v2) { e[j]->setVertex(k, v1); }
            }
          }
        }
        return true;
      }
      return false;
    }
    
    // collapse v1 & v2 to their middle and replace into e1 & e2
    static bool _tryToCollapseThatVertex(GFace *gf, std::vector<MElement *> &e1,
                                         std::vector<MElement *> &e2, MElement *q,
                                         MVertex *v1, MVertex *v2)
    {
      std::vector<MElement *> e = e1;
      e.insert(e.end(), e2.begin(), e2.end());
    
      double uu1, vv1;
      if(!v1->getParameter(0, uu1)) {
        return _tryToCollapseThatVertex2(gf, e1, e2, q, v1, v2);
      }
      v1->getParameter(1, vv1);
      double x1 = v1->x();
      double y1 = v1->y();
      double z1 = v1->z();
    
      double uu2, vv2;
      v2->getParameter(0, uu2);
      v2->getParameter(1, vv2);
      double x2 = v2->x();
      double y2 = v2->y();
      double z2 = v2->z();
    
      // new position of v1 && v2
      GPoint pp = gf->point(0.5 * (uu1 + uu2), 0.5 * (vv1 + vv2));
      double worst_quality_old = 1.0;
      double worst_quality_new = 1.0;
      int count = 0;
      for(std::size_t j = 0; j < e.size(); ++j) {
        if(e[j] != q) {
          count++;
          worst_quality_old = std::min(worst_quality_old, e[j]->etaShapeMeasure());
          v1->x() = pp.x();
          v1->y() = pp.y();
          v1->z() = pp.z();
          v1->setParameter(0, pp.u());
          v1->setParameter(1, pp.v());
          v2->x() = pp.x();
          v2->y() = pp.y();
          v2->z() = pp.z();
          v2->setParameter(0, pp.u());
          v2->setParameter(1, pp.v());
          worst_quality_new = std::min(worst_quality_new, e[j]->etaShapeMeasure());
          v1->x() = x1;
          v1->y() = y1;
          v1->z() = z1;
          v1->setParameter(0, uu1);
          v1->setParameter(1, vv1);
          v2->x() = x2;
          v2->y() = y2;
          v2->z() = z2;
          v2->setParameter(0, uu2);
          v1->setParameter(1, vv2);
        }
      }
    
      if(worst_quality_new > worst_quality_old) {
        v1->x() = pp.x();
        v1->y() = pp.y();
        v1->z() = pp.z();
        v1->setParameter(0, pp.u());
        v1->setParameter(1, pp.v());
        for(std::size_t j = 0; j < e.size(); ++j) {
          if(e[j] != q) {
            for(int k = 0; k < 4; k++) {
              if(e[j]->getVertex(k) == v2) { e[j]->setVertex(k, v1); }
            }
          }
        }
        return true;
      }
      return false;
    }
    
    static bool _isItAGoodIdeaToMoveThatVertex(GFace *gf,
                                               const std::vector<MElement *> &e1,
                                               MVertex *v1, const SPoint2 &before,
                                               const SPoint2 &after)
    {
      double surface_old = 0;
      double surface_new = 0;
    
      GPoint gp = gf->point(after);
      if(!gp.succeeded()) return false;
      SPoint3 pafter(gp.x(), gp.y(), gp.z());
      SPoint3 pbefore(v1->x(), v1->y(), v1->z());
    
      double minq = 1.0;
      for(std::size_t j = 0; j < e1.size(); ++j) {
        surface_old += surfaceFaceUV(e1[j], gf, false);
        minq = std::min(e1[j]->etaShapeMeasure(), minq);
      }
    
      v1->setParameter(0, after.x());
      v1->setParameter(1, after.y());
      v1->setXYZ(pafter.x(), pafter.y(), pafter.z());
    
      double minq_new = 1.0;
      for(std::size_t j = 0; j < e1.size(); ++j) {
        surface_new += surfaceFaceUV(e1[j], gf, false);
        minq_new = std::min(e1[j]->etaShapeMeasure(), minq_new);
      }
    
      v1->setParameter(0, before.x());
      v1->setParameter(1, before.y());
      v1->setXYZ(pbefore.x(), pbefore.y(), pbefore.z());
      if((1. + 1.e-10) * surface_old < surface_new || minq_new < minq) {
        return false;
      }
      return true;
    }
    
    static bool has_none_of(std::set<MVertex *> const &touched, MVertex *const v1,
                            MVertex *const v2, MVertex *const v3, MVertex *const v4)
    {
      return touched.find(v1) == touched.end() &&
             touched.find(v2) == touched.end() &&
             touched.find(v3) == touched.end() && touched.find(v4) == touched.end();
    }
    
    static bool are_all_on_surface(MVertex *const v1, MVertex *const v2,
                                   MVertex *const v3, MVertex *const v4, GFace *gf)
    {
      return v1->onWhat() == gf && v2->onWhat() == gf && v3->onWhat() == gf &&
             v4->onWhat() == gf;
    }
    
    template <class InputIterator>
    bool are_size_three(InputIterator iterator1, InputIterator iterator2)
    {
      return iterator1->second.size() == 3 && iterator2->second.size() == 3;
    }
    
    static int _removeDiamonds(GFace *const gf)
    {
      v2t_cont adj;
      buildVertexToElement(gf->quadrangles, adj);
    
      std::vector<MElement *> diamonds;
      std::set<MVertex *> touched;
      std::vector<MVertex *> deleted;
    
      std::vector<MQuadrangle *> quadrangles2;
      quadrangles2.reserve(gf->quadrangles.size());
    
      for(std::size_t i = 0; i < gf->triangles.size(); i++) {
        touched.insert(gf->triangles[i]->getVertex(0));
        touched.insert(gf->triangles[i]->getVertex(1));
        touched.insert(gf->triangles[i]->getVertex(2));
      }
    
      for(std::size_t i = 0; i < gf->quadrangles.size(); i++) {
        MQuadrangle *const q = gf->quadrangles[i];
    
        MVertex *const v1 = q->getVertex(0);
        MVertex *const v2 = q->getVertex(1);
        MVertex *const v3 = q->getVertex(2);
        MVertex *const v4 = q->getVertex(3);
    
        if(has_none_of(touched, v1, v2, v3, v4)) {
          auto it1 = adj.find(v1);
          auto it2 = adj.find(v2);
          auto it3 = adj.find(v3);
          auto it4 = adj.find(v4);
    
          if(are_all_on_surface(v1, v2, v3, v4, gf) && are_size_three(it1, it3) &&
             _tryToCollapseThatVertex(gf, it1->second, it3->second, q, v1, v3)) {
            touched.insert(v1);
            touched.insert(v2);
            touched.insert(v3);
            touched.insert(v4);
    
            deleted.push_back(v3);
    
            diamonds.push_back(q);
          }
          else if(are_all_on_surface(v1, v2, v3, v4, gf) &&
                  are_size_three(it2, it4) &&
                  _tryToCollapseThatVertex(gf, it2->second, it4->second, q, v2,
                                           v4)) {
            touched.insert(v1);
            touched.insert(v2);
            touched.insert(v3);
            touched.insert(v4);
    
            deleted.push_back(v4);
    
            diamonds.push_back(q);
          }
          else {
            quadrangles2.push_back(q);
          }
        }
        else {
          quadrangles2.push_back(q);
        }
      }
      std::sort(deleted.begin(), deleted.end());
      deleted.erase(std::unique(deleted.begin(), deleted.end()), deleted.end());
    
      gf->quadrangles = quadrangles2;
    
      std::vector<MVertex *> mesh_vertices2;
      mesh_vertices2.reserve(gf->mesh_vertices.size());
    
      for(std::size_t i = 0; i < gf->mesh_vertices.size(); i++) {
        if(!std::binary_search(deleted.begin(), deleted.end(),
                               gf->mesh_vertices[i])) {
          mesh_vertices2.push_back(gf->mesh_vertices[i]);
        }
        else {
          delete gf->mesh_vertices[i];
        }
      }
      gf->mesh_vertices = mesh_vertices2;
    
      std::sort(diamonds.begin(), diamonds.end());
    
      return std::distance(diamonds.begin(),
                           std::unique(diamonds.begin(), diamonds.end()));
    }
    
    int removeDiamonds(GFace *gf)
    {
      std::size_t nbRemove = 0;
      while(1) {
        int x = _removeDiamonds(gf);
        if(!x) break;
        nbRemove += x;
      }
      Msg::Debug("%i diamond quads removed", nbRemove);
      return nbRemove;
    }
    
    struct p1p2p3 {
      MVertex *p1, *p2;
    };
    
    static void _relocate(GFace *gf, MVertex *ver,
                          const std::vector<MElement *> &lt)
    {
      if(ver->onWhat() != gf) return;
      MFaceVertex *fv = dynamic_cast<MFaceVertex *>(ver);
      if(fv && fv->bl_data) return;
    
      double initu, initv;
      ver->getParameter(0, initu);
      ver->getParameter(1, initv);
    
      // compute the vertices connected to that one
      std::map<MVertex *, SPoint2, MVertexPtrLessThan> pts;
      for(std::size_t i = 0; i < lt.size(); i++) {
        for(int j = 0; j < lt[i]->getNumEdges(); j++) {
          MEdge e = lt[i]->getEdge(j);
          SPoint2 param0, param1;
          if(e.getVertex(0) == ver) {
            reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0,
                                  param1);
            pts[e.getVertex(1)] = param1;
          }
          else if(e.getVertex(1) == ver) {
            reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0,
                                  param1);
            pts[e.getVertex(0)] = param0;
          }
        }
      }
    
      SPoint2 before(initu, initv);
      double metric[3];
      SPoint2 after(0, 0);
      double COUNT = 0.0;
      for(auto it = pts.begin(); it != pts.end(); ++it) {
        SPoint2 adj = it->second;
        SVector3 d(adj.x() - before.x(), adj.y() - before.y(), 0.0);
        d.normalize();
        buildMetric(gf, adj, metric);
        const double F =
          sqrt(metric[0] * d.x() * d.x() + 2 * metric[1] * d.x() * d.y() +
               metric[2] * d.y() * d.y());
        after += adj * F;
        COUNT += F;
      }
      after *= (1.0 / COUNT);
      double FACTOR = 1.0;
      const int MAXITER = 5;
      SPoint2 actual = before;
      for(int ITER = 0; ITER < MAXITER; ITER++) {
        SPoint2 trial = after * FACTOR + before * (1. - FACTOR);
        bool success = _isItAGoodIdeaToMoveThatVertex(gf, lt, ver, actual, trial);
        if(success) {
          GPoint pt = gf->point(trial);
          if(pt.succeeded()) {
            actual = trial;
            ver->setParameter(0, trial.x());
            ver->setParameter(1, trial.y());
            ver->x() = pt.x();
            ver->y() = pt.y();
            ver->z() = pt.z();
          }
        }
        FACTOR /= 1.4;
      }
    }
    
    void getAllBoundaryLayerVertices(GFace *gf, std::set<MVertex *> &vs)
    {
      vs.clear();
      BoundaryLayerColumns *_columns = gf->getColumns();
      if(!_columns) return;
      for(auto it = _columns->_data.begin(); it != _columns->_data.end(); it++) {
        BoundaryLayerData &data = it->second;
        for(size_t i = 0; i < data._column.size(); i++) vs.insert(data._column[i]);
      }
    }
    
    void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
    {
      if((gf->triangles.size() > 0 && gf->triangles[0]->getPolynomialOrder() > 1) ||
         (gf->quadrangles.size() > 0 &&
          gf->quadrangles[0]->getPolynomialOrder() > 1)) {
        Msg::Error(
          "Surface mesh smoothing only valid for first order mesh (use the high-"
          "order optimization tools for high-order meshes)");
        return;
      }
    
      if(!niter) return;
    
      Msg::Debug("laplacian smoothing ...");
    
      std::set<MVertex *> vs;
      getAllBoundaryLayerVertices(gf, vs);
      v2t_cont adj;
      buildVertexToElement(gf->triangles, adj);
      buildVertexToElement(gf->quadrangles, adj);
      for(int i = 0; i < niter; i++) {
        auto it = adj.begin();
        while(it != adj.end()) {
          if(vs.find(it->first) == vs.end()) {
            _relocate(gf, it->first, it->second);
          }
          ++it;
        }
      }
    }
    
    static void _recombineIntoQuads(GFace *gf, bool blossom, bool cubicGraph = 1)
    {
      if(gf->triangles.empty()) return;
      if(gf->compound.size()) return;
    
      std::vector<MVertex *> emb_edgeverts;
      {
        std::vector<GEdge *> emb_edges = gf->getEmbeddedEdges();
        auto ite = emb_edges.begin();
        while(ite != emb_edges.end()) {
          if(!(*ite)->isMeshDegenerated()) {
            emb_edgeverts.insert(emb_edgeverts.end(), (*ite)->mesh_vertices.begin(),
                                 (*ite)->mesh_vertices.end());
            if((*ite)->getBeginVertex())
              emb_edgeverts.insert(emb_edgeverts.end(),
                                   (*ite)->getBeginVertex()->mesh_vertices.begin(),
                                   (*ite)->getBeginVertex()->mesh_vertices.end());
            if((*ite)->getEndVertex())
              emb_edgeverts.insert(emb_edgeverts.end(),
                                   (*ite)->getEndVertex()->mesh_vertices.begin(),
                                   (*ite)->getEndVertex()->mesh_vertices.end());
          }
          ++ite;
        }
      }
    
      {
        std::vector<GEdge *> const &_edges = gf->edges();
        auto ite = _edges.begin();
        while(ite != _edges.end()) {
          if(!(*ite)->isMeshDegenerated()) {
            if((*ite)->isSeam(gf)) {
              emb_edgeverts.insert(emb_edgeverts.end(),
                                   (*ite)->mesh_vertices.begin(),
                                   (*ite)->mesh_vertices.end());
              if((*ite)->getBeginVertex())
                emb_edgeverts.insert(
                  emb_edgeverts.end(),
                  (*ite)->getBeginVertex()->mesh_vertices.begin(),
                  (*ite)->getBeginVertex()->mesh_vertices.end());
              if((*ite)->getEndVertex())
                emb_edgeverts.insert(emb_edgeverts.end(),
                                     (*ite)->getEndVertex()->mesh_vertices.begin(),
                                     (*ite)->getEndVertex()->mesh_vertices.end());
            }
          }
          ++ite;
        }
      }
    
      // sort and erase the duplicates
      std::sort(emb_edgeverts.begin(), emb_edgeverts.end());
      emb_edgeverts.erase(std::unique(emb_edgeverts.begin(), emb_edgeverts.end()),
                          emb_edgeverts.end());
    
      e2t_cont adj;
      buildEdgeToElement(gf->triangles, adj);
    
      FieldManager *fields = gf->model()->getFields();
      Field *cross_field = NULL;
      SVector3 t1;
    
      if(fields->getBackgroundField() > 0) {
        cross_field = fields->get(fields->getBackgroundField());
        if(cross_field->numComponents() != 3) { // we hae a true scaled cross fiel
          cross_field = NULL;
        }
      }
    
      std::vector<RecombineTriangle> pairs;
    
      std::map<MVertex *, std::pair<MElement *, MElement *> > makeGraphPeriodic;
    
      for(auto it = adj.begin(); it != adj.end(); ++it) {
        if(it->second.second && it->second.first->getNumVertices() == 3 &&
           it->second.second->getNumVertices() == 3 &&
           (!std::binary_search(emb_edgeverts.begin(), emb_edgeverts.end(),
                                it->first.getVertex(0)) ||
            !std::binary_search(emb_edgeverts.begin(), emb_edgeverts.end(),
                                it->first.getVertex(1)))) {
          pairs.push_back(RecombineTriangle(it->first, it->second.first,
                                            it->second.second, cross_field));
        }
        else if(!it->second.second && it->second.first->getNumVertices() == 3) {
          for(int i = 0; i < 2; i++) {
            MVertex *const v = it->first.getVertex(i);
            auto itv = makeGraphPeriodic.find(v);
            if(itv == makeGraphPeriodic.end()) {
              makeGraphPeriodic[v] =
                std::make_pair(it->second.first, static_cast<MElement *>(nullptr));
            }
            else {
              if(itv->second.first != it->second.first)
                itv->second.second = it->second.first;
              else
                makeGraphPeriodic.erase(itv);
            }
          }
        }
      }
    
      std::sort(pairs.begin(), pairs.end());
      std::set<MElement *> touched;
    
      if(blossom) {
    #if defined(HAVE_BLOSSOM)
        int ncount = gf->triangles.size();
        if(ncount % 2 != 0) {
          Msg::Warning("Cannot apply Blossom: odd number of triangles (%d) in "
                       "surface %d",
                       ncount, gf->tag());
        }
        else {
          int ecount =
            cubicGraph ? (pairs.size() + makeGraphPeriodic.size()) : pairs.size();
          Msg::Info("Blossom: %d internal %d closed", (int)pairs.size(),
                    (int)makeGraphPeriodic.size());
          Msg::Debug("Perfect Match Starts %d edges %d nodes", ecount, ncount);
          std::map<MElement *, int> t2n;
          std::map<int, MElement *> n2t;
          for(std::size_t i = 0; i < gf->triangles.size(); ++i) {
            t2n[gf->triangles[i]] = i;
            n2t[i] = gf->triangles[i];
          }
          // do not use new[] here, blossom will free it with free() and not with
          // delete
          int *elist = (int *)malloc(sizeof(int) * 2 * ecount);
          int *elen = (int *)malloc(sizeof(int) * ecount);
    
          for(std::size_t i = 0; i < pairs.size(); ++i) {
            elist[2 * i] = t2n[pairs[i].t1];
            elist[2 * i + 1] = t2n[pairs[i].t2];
            elen[i] = (int)1000 * std::exp(-pairs[i].angle);
            int NB = 0;
            if(pairs[i].n1->onWhat()->dim() < 2) NB++;
            if(pairs[i].n2->onWhat()->dim() < 2) NB++;
            if(pairs[i].n3->onWhat()->dim() < 2) NB++;
            if(pairs[i].n4->onWhat()->dim() < 2) NB++;
            if(elen[i] > static_cast<int>(1000 * std::exp(0.1)) && NB > 2) {
              elen[i] = 5000;
            }
            else if(elen[i] >= 1000 && NB > 2) {
              elen[i] = 10000;
            }
          }
    
          if(cubicGraph) {
            auto itv = makeGraphPeriodic.begin();
            std::size_t CC = pairs.size();
            for(; itv != makeGraphPeriodic.end(); ++itv) {
              elist[2 * CC] = t2n[itv->second.first];
              elist[2 * CC + 1] = t2n[itv->second.second];
              elen[CC++] = 100000;
            }
          }
    
          double matzeit = 0.0;
          char MATCHFILE[256];
          sprintf(MATCHFILE, ".face.match");
          if(perfect_match(ncount, nullptr, ecount, &elist, &elen, nullptr,
                           MATCHFILE, 0, 0, 0, 0, &matzeit)) {
            Msg::Error(
              "Perfect Match failed in quadrangulation, try something else");
            free(elist);
            pairs.clear();
          }
          else {
            // TEST
            for(int k = 0; k < elist[0]; k++) {
              int i1 = elist[1 + 3 * k], i2 = elist[1 + 3 * k + 1],
                  an = elist[1 + 3 * k + 2];
              // FIXME !
              if(an == 100000 /*|| an == 1000*/) {
                // toProcess.push_back(std::make_pair(n2t[i1],n2t[i2]));
                // Msg::Warning("Extra edge found in blossom algorithm,
                // optimization"
                //              "will be required");
              }
              else {
                MElement *t1 = n2t[i1];
                MElement *t2 = n2t[i2];
                touched.insert(t1);
                touched.insert(t2);
                MVertex *other = nullptr;
                for(int i = 0; i < 3; i++) {
                  if(t1->getVertex(0) != t2->getVertex(i) &&
                     t1->getVertex(1) != t2->getVertex(i) &&
                     t1->getVertex(2) != t2->getVertex(i)) {
                    other = t2->getVertex(i);
                    break;
                  }
                }
                int start = 0;
                for(int i = 0; i < 3; i++) {
                  if(t2->getVertex(0) != t1->getVertex(i) &&
                     t2->getVertex(1) != t1->getVertex(i) &&
                     t2->getVertex(2) != t1->getVertex(i)) {
                    start = i;
                    break;
                  }
                }
                MQuadrangle *q = new MQuadrangle(
                  t1->getVertex(start), t1->getVertex((start + 1) % 3), other,
                  t1->getVertex((start + 2) % 3));
                gf->quadrangles.push_back(q);
              }
            }
            free(elist);
            pairs.clear();
            Msg::Debug("Perfect Match Succeeded in Quadrangulation (%g sec)",
                       matzeit);
          }
        }
    #else
        Msg::Warning("Gmsh should be compiled with the Blossom IV code and "
                     " CONCORDE in order to allow the Blossom optimization");
    #endif
      }
    
      // simple greedy recombination
      auto itp = pairs.begin();
      while(itp != pairs.end()) {
        if(itp->angle < gf->meshAttributes.recombineAngle) {
          MElement *t1 = itp->t1;
          MElement *t2 = itp->t2;
          if(touched.find(t1) == touched.end() &&
             touched.find(t2) == touched.end()) {
            touched.insert(t1);
            touched.insert(t2);
            int orientation = 0;
            for(int i = 0; i < 3; i++) {
              if(t1->getVertex(i) == itp->n1) {
                if(t1->getVertex((i + 1) % 3) == itp->n2)
                  orientation = 1;
                else
                  orientation = -1;
                break;
              }
            }
            gf->quadrangles.push_back(
              new MQuadrangle(itp->n1, orientation < 0 ? itp->n3 : itp->n4, itp->n2,
                              orientation < 0 ? itp->n4 : itp->n3));
          }
        }
        ++itp;
      }
    
      std::vector<MTriangle *> triangles2;
      triangles2.reserve(gf->triangles.size());
      for(std::size_t i = 0; i < gf->triangles.size(); i++) {
        if(touched.find(gf->triangles[i]) == touched.end()) {
          triangles2.push_back(gf->triangles[i]);
        }
        else {
          delete gf->triangles[i];
        }
      }
      gf->triangles = triangles2;
    }
    
    static double printStats(GFace *gf, const char *message)
    {
      int nbBad = 0;
      int nbInv = 0;
      double Qav = 0;
      double Qmin = 1;
      for(std::size_t i = 0; i < gf->quadrangles.size(); i++) {
        double Q = gf->quadrangles[i]->etaShapeMeasure();
        if(Q <= 0.0) nbInv++;
        if(Q <= 0.1) nbBad++;
        Qav += Q;
        Qmin = std::min(Q, Qmin);
      }
      Msg::Info(
        "%s: %d quads, %d triangles, %d invalid quads, %d quads with Q < 0.1, "
        "avg Q = %g, min Q = %g",
        message, gf->quadrangles.size(), gf->triangles.size(), nbInv, nbBad,
        Qav / gf->quadrangles.size(), Qmin);
      return Qmin;
    }
    
    // The topological optimization routines assume that a full topology of the
    // model exists. When reading multi-surface STL files for example, if
    // CreateTopology or ReclassifySurfaces is not called, quads can have nodes
    // owned by an adjacent surface. Since the topological optimization routines
    // remove nodes, this will produce an invalide model mesh (and crash).
    static bool _isModelOkForTopologicalOpti(GModel *m)
    {
      for(auto it = m->firstFace(); it != m->lastFace(); it++) {
        GFace *gf = *it;
        for(std::size_t j = 0; j < gf->getNumMeshElements(); j++) {
          MElement *e = gf->getMeshElement(j);
          for(std::size_t k = 0; k < e->getNumVertices(); k++) {
            GEntity *ge = e->getVertex(k)->onWhat();
            if(!ge) return false;
            if(ge->dim() == 2 && ge != gf) return false;
          }
        }
      }
      return true;
    }
    
    void recombineIntoQuads(GFace *gf, bool blossom, int topologicalOptiPasses,
                            bool nodeRepositioning, double minqual)
    {
      double t1 = Cpu(), w1 = TimeOfDay();
    
      bool haveParam = (gf->geomType() != GEntity::DiscreteSurface);
      bool debug = (Msg::GetVerbosity() == 99);
    
      if(debug) gf->model()->writeMSH("recombine_0before.msh");
    
      _recombineIntoQuads(gf, blossom);
    
      if(debug) gf->model()->writeMSH("recombine_1raw.msh");
    
      if(haveParam && nodeRepositioning) {
        RelocateVertices(gf, CTX::instance()->mesh.nbSmoothing);
        if(debug) gf->model()->writeMSH("recombine_2smoothed.msh");
      }
    
    // FIXME: not thread-safe
    #pragma omp critical
      {
        if(topologicalOptiPasses > 0) {
          if(!_isModelOkForTopologicalOpti(gf->model())) {
            Msg::Info
              ("Skipping topological optimization - mesh topology is not complete");
          }
          else {
            int iter = 0, nbTwoQuadNodes = 1, nbDiamonds = 1;
            while(nbTwoQuadNodes || nbDiamonds) {
              Msg::Debug("Topological optimization of quad mesh: pass %d", iter);
              nbTwoQuadNodes = removeTwoQuadsNodes(gf);
              // removeDiamonds uses the parametrization or searches for closest point
              nbDiamonds = haveParam ? removeDiamonds(gf) : 0;
              if(haveParam && nodeRepositioning)
                RelocateVertices(gf, CTX::instance()->mesh.nbSmoothing);
              iter++;
              if(iter > topologicalOptiPasses) break;
            }
            if(debug) gf->model()->writeMSH("recombine_3topo.msh");
          }
        }
      }
    
      // re-split bad quads into triangles
      if(minqual > 0)
        quadsToTriangles(gf, minqual);
    
      if(debug) gf->model()->writeMSH("recombine_4quality.msh");
    
      if(haveParam && nodeRepositioning)
        RelocateVertices(gf, CTX::instance()->mesh.nbSmoothing);
    
      double t2 = Cpu(), w2 = TimeOfDay();
    
      char name[256];
      sprintf(name, "%s recombination completed (Wall %gs, CPU %gs)",
              blossom ? "Blossom" : "Simple", w2 - w1, t2 - t1);
      printStats(gf, name);
    
      if(debug) gf->model()->writeMSH("recombine_5final.msh");
    }
    
    void quadsToTriangles(GFace *gf, double minqual)
    {
      std::vector<MQuadrangle *> qds;
      std::map<MElement *, std::pair<MElement *, MElement *> > change;
      for(std::size_t i = 0; i < gf->quadrangles.size(); i++) {
        MQuadrangle *q = gf->quadrangles[i];
        if(q->etaShapeMeasure() < minqual + 1e-12) {
          MTriangle *t11 =
            new MTriangle(q->getVertex(0), q->getVertex(1), q->getVertex(2));
          MTriangle *t12 =
            new MTriangle(q->getVertex(2), q->getVertex(3), q->getVertex(0));
          MTriangle *t21 =
            new MTriangle(q->getVertex(1), q->getVertex(2), q->getVertex(3));
          MTriangle *t22 =
            new MTriangle(q->getVertex(3), q->getVertex(0), q->getVertex(1));
          double qual1 =
            std::min(t11->gammaShapeMeasure(), t12->gammaShapeMeasure());
          double qual2 =
            std::min(t21->gammaShapeMeasure(), t22->gammaShapeMeasure());
          double ori2 = dot(t21->getFace(0).normal(), t22->getFace(0).normal());
          // choose (t11, t12) if it leads to the best quality OR if choosing (t21,
          // t22) would revert the orientation (which can happen if q is not convex)
          if(qual1 > qual2 || ori2 < 0) {
            gf->triangles.push_back(t11);
            gf->triangles.push_back(t12);
            change[q] = std::make_pair(t11, t12);
            delete t21;
            delete t22;
          }
          else {
            gf->triangles.push_back(t21);
            gf->triangles.push_back(t22);
            change[q] = std::make_pair(t21, t22);
            delete t11;
            delete t12;
          }
          delete q;
        }
        else {
          qds.push_back(q);
        }
      }
      gf->quadrangles = qds;
    
      BoundaryLayerColumns *_columns = gf->getColumns();
      if(!_columns) return;
    
      // Update the data struture for boundary layers
      // WARNING: First quad element is replaced by one of the two triangles,
      // without taking care of if it is the truly the first one or not.
    
      //  std::map<MElement*,MElement*> _toFirst;
      std::map<MElement *, std::vector<MElement *> > newElemColumns;
    
      for(auto it = _columns->_elemColumns.begin();
          it != _columns->_elemColumns.end(); it++) {
        MElement *firstEl = it->first;
        auto it2 = change.find(firstEl);
        if(it2 != change.end()) firstEl = it2->second.first;
        // it2->second.first may be the one that touch boundary or not...
    
        std::vector<MElement *> &newColumn = newElemColumns[firstEl];
        std::vector<MElement *> &oldColumn = it->second;
    
        for(std::size_t i = 0; i < oldColumn.size(); i++) {
          MElement *oldEl = oldColumn[i];
          it2 = change.find(oldEl);
          if(it2 == change.end()) {
            newColumn.push_back(oldEl);
            _columns->_toFirst[oldEl] = firstEl;
          }
          else {
            newColumn.push_back(it2->second.first);
            newColumn.push_back(it2->second.second);
            _columns->_toFirst.erase(oldEl);
            _columns->_toFirst[it2->second.first] = firstEl;
            _columns->_toFirst[it2->second.second] = firstEl;
          }
        }
      }
      _columns->_elemColumns = newElemColumns;
    }
    
    void splitElementsInBoundaryLayerIfNeeded(GFace *gf)
    {
      if(!CTX::instance()->mesh.recombineAll && !gf->meshAttributes.recombine) {
        int numSplit = 0;
        int numNoSplit = 0;
        FieldManager *fields = gf->model()->getFields();
        int n = fields->getNumBoundaryLayerFields();
        for(int i = 0; i < n; ++i) {
          Field *bl_field = fields->get(fields->getBoundaryLayerField(i));
          if(bl_field == nullptr) continue;
          BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField *>(bl_field);
          if(blf->iRecombine)
            ++numNoSplit;
          else
            ++numSplit;
        }
    
        if(numSplit > 0 && numNoSplit > 0)
          Msg::Warning("Cannot generate simplicial and non-simplicial boundary "
                       "layers together. Keeping them non-simplicial...");
        if(numNoSplit == 0 && numSplit > 0) quadsToTriangles(gf, 10000);
      }
    }