meshGFaceOptimize.cpp

Go to the documentation of this file.

// 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);
  }
}