text: // Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
  //
  // See the LICENSE.txt file for license information. Please report all
  // issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
  
  #include <sstream>
  #include "GmshConfig.h"
  #include "GmshMessage.h"
  #include "GModel.h"
  #include "GFace.h"
  #include "GEdge.h"
  #include "MTriangle.h"
  #include "MQuadrangle.h"
  #include "MElementCut.h"
  #include "VertexArray.h"
  #include "fullMatrix.h"
  #include "Numeric.h"
  #include "GaussLegendre1D.h"
  #include "Context.h"
  #include "OS.h"
  #include "discreteFace.h"
  
  #if defined(HAVE_MESH)
  #include "meshGFace.h"
  #include "meshGFaceOptimize.h"
  #include "BackgroundMeshTools.h"
  #endif
  
  #if defined(HAVE_ALGLIB)
  #include <stdafx.h>
  #include <optimization.h>
  #endif
  
  // TODO C++11 remove macro
  #define SQU(a) ((a) * (a))
  
  GFace::GFace(GModel *model, int tag)
      : GEntity(model, tag), r1(0), r2(0), va_geom_triangles(0)
  {
      meshStatistics.status = GFace::PENDING;
      meshStatistics.refineAllEdges = false;
      GFace::resetMeshAttributes();
  }
  
  GFace::~GFace()
  {
      std::vector<GEdge *>::iterator it = l_edges.begin();
      while(it != l_edges.end()) {
          (*it)->delFace(this);
          ++it;
      }
  
      if(va_geom_triangles) delete va_geom_triangles;
  
      GFace::deleteMesh();
  }
  
  int GFace::getCurvatureControlParameter() const
  {
      std::map<int, int>::iterator it =
              CTX::instance()->mesh.curvatureControlPerFace.find(tag());
      return it == CTX::instance()->mesh.curvatureControlPerFace.end() ?
                  CTX::instance()->mesh.minCircPoints :
                  it->second;
  }
  
  void GFace::setCurvatureControlParameter(int n)
  {
      CTX::instance()->mesh.curvatureControlPerFace[tag()] = n;
  }
  
  int GFace::getMeshingAlgo() const
  {
      std::map<int, int>::iterator it =
              CTX::instance()->mesh.algo2dPerFace.find(tag());
      return it == CTX::instance()->mesh.algo2dPerFace.end() ?
                  CTX::instance()->mesh.algo2d :
                  it->second;
  }
  
  void GFace::setMeshingAlgo(int algo) const
  {
      CTX::instance()->mesh.algo2dPerFace[tag()] = algo;
  }
  
  void GFace::unsetMeshingAlgo() const
  {
      CTX::instance()->mesh.algo2dPerFace.erase(tag());
  }
  
  void GFace::delFreeEdge(GEdge *edge)
  {
      // delete the edge from the edge list and the orientation list
      std::vector<GEdge *>::iterator ite = l_edges.begin();
      std::vector<int>::iterator itd = l_dirs.begin();
      while(ite != l_edges.end()) {
          if(edge == *ite) {
              Msg::Debug("Erasing edge %d from edge list in face %d", edge->tag(),
                         tag());
              l_edges.erase(ite);
              if(itd != l_dirs.end()) l_dirs.erase(itd);
              break;
          }
          ite++;
          if(itd != l_dirs.end()) itd++;
      }
  
      // delete the edge from the edge loops
      for(std::list<GEdgeLoop>::iterator it = edgeLoops.begin();
          it != edgeLoops.end(); it++) {
          for(GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); it2++) {
              if(edge == it2->ge) {
                  Msg::Debug("Erasing edge %d from edge loop in face %d", edge->tag(),
                             tag());
                  it->erase(it2);
                  break;
              }
          }
      }
  }
  
  int GFace::delEdge(GEdge *edge)
  {
      // BUG If the iterator is equal to end() then the erase will be ill-formed
      // TODO C++11 fix this UB
      std::vector<GEdge *>::iterator it;
      int pos = 0;
      for(it = l_edges.begin(); it != l_edges.end(); ++it) {
          if(*it == edge) break;
          pos++;
      }
      l_edges.erase(it);
  
      std::vector<int>::iterator itOri;
      int posOri = 0, orientation = 0;
      for(itOri = l_dirs.begin(); itOri != l_dirs.end(); ++itOri) {
          if(posOri == pos) {
              orientation = *itOri;
              break;
          }
          posOri++;
      }
      l_dirs.erase(itOri);
  
      return orientation;
  }
  
  void GFace::deleteMesh(bool onlyDeleteElements)
  {
      if(!onlyDeleteElements) {
          for(std::size_t i = 0; i < mesh_vertices.size(); i++)
              delete mesh_vertices[i];
          mesh_vertices.clear();
          transfinite_vertices.clear();
      }
      for(std::size_t i = 0; i < triangles.size(); i++) delete triangles[i];
      triangles.clear();
      for(std::size_t i = 0; i < quadrangles.size(); i++) delete quadrangles[i];
      quadrangles.clear();
      for(std::size_t i = 0; i < polygons.size(); i++) delete polygons[i];
      polygons.clear();
      deleteVertexArrays();
      model()->destroyMeshCaches();
  }
  
  std::size_t GFace::getNumMeshElements() const
  {
      return triangles.size() + quadrangles.size() + polygons.size();
  }
  
  std::size_t GFace::getNumMeshElementsByType(const int familyType) const
  {
      if(familyType == TYPE_TRI)
          return triangles.size();
      else if(familyType == TYPE_QUA)
          return quadrangles.size();
      else if(familyType == TYPE_POLYG)
          return polygons.size();
  
      return 0;
  }
  
  std::size_t GFace::getNumMeshParentElements()
  {
      std::size_t n = 0;
      for(std::size_t i = 0; i < polygons.size(); i++)
          if(polygons[i]->ownsParent()) n++;
      return n;
  }
  
  void GFace::getNumMeshElements(unsigned *const c) const
  {
      c[0] += triangles.size();
      c[1] += quadrangles.size();
      c[2] += polygons.size();
  }
  
  MElement *const *GFace::getStartElementType(int type) const
  {
      switch(type) {
      case 0:
          if(triangles.empty()) return 0; // msvc would throw an exception
          return reinterpret_cast<MElement *const *>(&triangles[0]);
      case 1:
          if(quadrangles.empty()) return 0; // msvc would throw an exception
          return reinterpret_cast<MElement *const *>(&quadrangles[0]);
      case 2:
          if(polygons.empty()) return 0;
          return reinterpret_cast<MElement *const *>(&polygons[0]);
      }
      return 0;
  }
  
  MElement *GFace::getMeshElement(std::size_t index) const
  {
      if(index < triangles.size())
          return triangles[index];
      else if(index < triangles.size() + quadrangles.size())
          return quadrangles[index - triangles.size()];
      else if(index < triangles.size() + quadrangles.size() + polygons.size())
          return polygons[index - triangles.size() - quadrangles.size()];
      return 0;
  }
  
  MElement *GFace::getMeshElementByType(const int familyType,
                                        const std::size_t index) const
  {
      if(familyType == TYPE_TRI)
          return triangles[index];
      else if(familyType == TYPE_QUA)
          return quadrangles[index];
      else if(familyType == TYPE_POLYG)
          return polygons[index];
  
      return 0;
  }
  
  void GFace::resetMeshAttributes()
  {
      meshAttributes.recombine = 0;
      meshAttributes.recombineAngle = 45.;
      meshAttributes.method = MESH_UNSTRUCTURED;
      meshAttributes.transfiniteArrangement = 0;
      meshAttributes.transfiniteSmoothing = -1;
      meshAttributes.extrude = 0;
      meshAttributes.reverseMesh = false;
      meshAttributes.meshSize = MAX_LC;
  }
  
  SBoundingBox3d GFace::bounds(bool fast) const
  {
      SBoundingBox3d res;
      if(geomType() != DiscreteSurface && geomType() != PartitionSurface) {
          // TODO C++11 std::accumulate
          std::vector<GEdge *>::const_iterator it = l_edges.begin();
          for(; it != l_edges.end(); it++) res += (*it)->bounds(fast);
      }
      else {
          int ipp = getNumMeshElements() / 20;
          if(ipp < 1) ipp = 1;
          for(std::size_t i = 0; i < getNumMeshElements(); i += ipp)
              for(std::size_t j = 0; j < getMeshElement(i)->getNumVertices(); j++)
                  res += getMeshElement(i)->getVertex(j)->point();
      }
      return res;
  }
  
  SOrientedBoundingBox GFace::getOBB()
  {
      if(!_obb) {
          std::vector<SPoint3> vertices;
          if(getNumMeshVertices() > 0) {
              int N = getNumMeshVertices();
              for(int i = 0; i < N; i++) {
                  MVertex *mv = getMeshVertex(i);
                  vertices.push_back(mv->point());
              }
              std::vector<GEdge *> const &eds = edges();
              for(std::vector<GEdge *>::const_iterator ed = eds.begin();
                  ed != eds.end(); ed++) {
                  int N2 = (*ed)->getNumMeshVertices();
                  for(int i = 0; i < N2; i++) {
                      MVertex *mv = (*ed)->getMeshVertex(i);
                      vertices.push_back(mv->point());
                  }
                  // Don't forget to add the first and last vertices...
                  if((*ed)->getBeginVertex()){
                      SPoint3 pt1((*ed)->getBeginVertex()->x(),
                                  (*ed)->getBeginVertex()->y(),
                                  (*ed)->getBeginVertex()->z());
                      vertices.push_back(pt1);
                  }
                  if((*ed)->getEndVertex()){
                      SPoint3 pt2((*ed)->getEndVertex()->x(),
                                  (*ed)->getEndVertex()->y(),
                                  (*ed)->getEndVertex()->z());
                      vertices.push_back(pt2);
                  }
              }
          }
          else if(buildSTLTriangulation()) {
              vertices = stl_vertices_xyz;
          }
          else {
              // Fallback, if we can't make a STL triangulation of the surface, use its
              // edges..
              std::vector<GEdge *> const &b_edges = edges();
              int N = 10;
              for(std::vector<GEdge *>::const_iterator b_edge = b_edges.begin();
                  b_edge != b_edges.end(); b_edge++) {
                  Range<double> tr = (*b_edge)->parBounds(0);
                  for(int j = 0; j < N; j++) {
                      double t =
                              tr.low() + (double)j / (double)(N - 1) * (tr.high() - tr.low());
                      GPoint p = (*b_edge)->point(t);
                      SPoint3 pt(p.x(), p.y(), p.z());
                      vertices.push_back(pt);
                  }
              }
          }
          _obb = SOrientedBoundingBox::buildOBB(vertices);
      }
      return SOrientedBoundingBox(_obb);
  }
  
  std::vector<MVertex *> GFace::getEmbeddedMeshVertices() const
  {
      std::set<MVertex *> tmp;
      for(std::vector<GEdge *>::const_iterator it = embedded_edges.begin();
          it != embedded_edges.end(); it++) {
          tmp.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
          if((*it)->getBeginVertex())
              tmp.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
                         (*it)->getBeginVertex()->mesh_vertices.end());
          if((*it)->getEndVertex())
              tmp.insert((*it)->getEndVertex()->mesh_vertices.begin(),
                         (*it)->getEndVertex()->mesh_vertices.end());
      }
      for(std::set<GVertex *>::const_iterator it = embedded_vertices.begin();
          it != embedded_vertices.end(); it++) {
          tmp.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
      }
      return std::vector<MVertex *>(tmp.begin(), tmp.end());
  }
  
  std::vector<GVertex *> GFace::vertices() const
  {
      std::set<GVertex *> v;
      for(std::vector<GEdge *>::const_iterator it = l_edges.begin();
          it != l_edges.end(); ++it) {
          GVertex *const v1 = (*it)->getBeginVertex();
          if(v1) v.insert(v1);
  
          GVertex *const v2 = (*it)->getEndVertex();
          if(v2) v.insert(v2);
      }
      return std::vector<GVertex *>(v.begin(), v.end());
  }
  
  void GFace::setVisibility(char val, bool recursive)
  {
      GEntity::setVisibility(val);
      if(recursive) {
          for(std::vector<GEdge *>::iterator it = l_edges.begin();
              it != l_edges.end(); ++it)
              (*it)->setVisibility(val, recursive);
      }
  }
  
  void GFace::setColor(unsigned int val, bool recursive)
  {
      GEntity::setColor(val);
      if(recursive) {
          for(std::vector<GEdge *>::iterator it = l_edges.begin();
              it != l_edges.end(); ++it)
              (*it)->setColor(val, recursive);
      }
  }
  
  std::string GFace::getAdditionalInfoString(bool multline)
  {
      std::ostringstream sstream;
      if(l_edges.size() > 20) {
          sstream << "Boundary curves: " << l_edges.front()->tag() << ", ...,"
                  << l_edges.back()->tag();
          if(multline)
              sstream << "\n";
          else
              sstream << " ";
      }
      else if(l_edges.size()) {
          sstream << "Boundary curves: ";
          for(std::vector<GEdge *>::iterator it = l_edges.begin();
              it != l_edges.end(); ++it) {
              if(it != l_edges.begin()) sstream << ", ";
              sstream << (*it)->tag();
          }
          if(multline)
              sstream << "\n";
          else
              sstream << " ";
      }
  
      if(embedded_edges.size()) {
          sstream << "Embedded curves: ";
          for(std::vector<GEdge *>::iterator it = embedded_edges.begin();
              it != embedded_edges.end(); ++it) {
              if(it != embedded_edges.begin()) sstream << ", ";
              sstream << (*it)->tag();
          }
          if(multline)
              sstream << "\n";
          else
              sstream << " ";
      }
  
      if(embedded_vertices.size()) {
          sstream << "Embedded points: ";
          for(std::set<GVertex *>::iterator it = embedded_vertices.begin();
              it != embedded_vertices.end(); ++it) {
              if(it != embedded_vertices.begin()) sstream << ", ";
              sstream << (*it)->tag();
          }
          if(multline)
              sstream << "\n";
          else
              sstream << " ";
      }
  
      if(meshAttributes.recombine || meshAttributes.method == MESH_TRANSFINITE ||
              meshAttributes.extrude || meshAttributes.reverseMesh) {
          sstream << "Mesh attributes:";
          if(meshAttributes.recombine) sstream << " recombined";
          if(meshAttributes.method == MESH_TRANSFINITE) sstream << " transfinite";
          if(meshAttributes.extrude) sstream << " extruded";
          if(meshAttributes.reverseMesh) sstream << " reverse";
      }
      std::string str = sstream.str();
      if(str.size() && (str[str.size() - 1] == '\n' || str[str.size() - 1] == ' '))
          str.resize(str.size() - 1);
      return str;
  }
  
  void GFace::writeGEO(FILE *fp)
  {
      if(geomType() == DiscreteSurface) return;
  
      std::vector<GEdge *> const &edg = edges();
      std::vector<int> const &dir = orientations();
      if(edg.size() && dir.size() == edg.size()) {
          std::vector<int> num, ori;
          for(std::vector<GEdge *>::const_iterator it = edg.begin(); it != edg.end();
              it++)
              num.push_back((*it)->tag());
          for(std::vector<int>::const_iterator it = dir.begin(); it != dir.end();
              it++)
              ori.push_back((*it) > 0 ? 1 : -1);
          fprintf(fp, "Line Loop(%d) = ", tag());
          for(std::size_t i = 0; i < num.size(); i++) {
              if(i)
                  fprintf(fp, ", %d", num[i] * ori[i]);
              else
                  fprintf(fp, "{%d", num[i] * ori[i]);
          }
          fprintf(fp, "};\n");
          if(geomType() == GEntity::Plane) {
              fprintf(fp, "Plane Surface(%d) = {%d};\n", tag(), tag());
          }
          else if(edg.size() == 3 || edg.size() == 4) {
              fprintf(fp, "Surface(%d) = {%d};\n", tag(), tag());
          }
          else {
              Msg::Error("Skipping surface %d in export", tag());
          }
      }
  
      for(std::vector<GEdge *>::iterator it = embedded_edges.begin();
          it != embedded_edges.end(); it++)
          fprintf(fp, "Line {%d} In Surface {%d};\n", (*it)->tag(), tag());
  
      for(std::set<GVertex *>::iterator it = embedded_vertices.begin();
          it != embedded_vertices.end(); it++)
          fprintf(fp, "Point {%d} In Surface {%d};\n", (*it)->tag(), tag());
  
      if(meshAttributes.method == MESH_TRANSFINITE) {
          fprintf(fp, "Transfinite Surface {%d}", tag());
          if(meshAttributes.corners.size()) {
              fprintf(fp, " = {");
              for(std::size_t i = 0; i < meshAttributes.corners.size(); i++) {
                  if(i) fprintf(fp, ",");
                  fprintf(fp, "%d", meshAttributes.corners[i]->tag());
              }
              fprintf(fp, "}");
          }
          fprintf(fp, ";\n");
      }
  
      if(meshAttributes.recombine) fprintf(fp, "Recombine Surface {%d};\n", tag());
  
      if(meshAttributes.reverseMesh) fprintf(fp, "Reverse Surface {%d};\n", tag());
  }
  
  void GFace::computeMeanPlane()
  {
      std::vector<SPoint3> pts;
  
      if(geomType() == Plane) {
          // Special case for planar CAD surfaces: we first try to compute the
          // parametrization in a way that is more robust with respect to
          // perturbations of the boundary. This is crucial for sensitivity analyses
          // and GFace::relocateMeshVertices(): after perturbation of the boundary, we
          // want a parametrization of the surface that is "close" to the original
          // one. If this fails, we fallback to the classical (SVD-based) algorithm.
          std::vector<GEdge *> const &edg = edges();
          for(std::vector<GEdge *>::const_iterator ite = edg.begin();
              ite != edg.end(); ite++) {
              const GEdge *e = *ite;
              if(e->geomType() == GEntity::DiscreteCurve ||
                      e->geomType() == GEntity::BoundaryLayerCurve) {
                  pts.clear();
                  break;
              }
              Range<double> b = e->parBounds(0);
              GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low()));
              pts.push_back(SPoint3(p1.x(), p1.y(), p1.z()));
              GPoint p2 = e->point(b.low() + 0.666 * (b.high() - b.low()));
              pts.push_back(SPoint3(p2.x(), p2.y(), p2.z()));
          }
          bool ok = false;
          double res[4] = {0., 0., 0., 0.}, xm = 0., ym = 0., zm = 0.;
          if(pts.size() >= 3) {
              SVector3 d01(pts[0], pts[1]);
              for(std::size_t i = 2; i < pts.size(); i++) {
                  SVector3 d0i(pts[0], pts[i]);
                  SVector3 n = crossprod(d01, d0i);
                  // if too small, the points are almost colinear; tolerance is relatively
                  // high so that we don't accept points on plane surfaces defined by
                  // lines that are not exactly co-planar
                  if(norm(n) >
                          sqrt(CTX::instance()->geom.tolerance) * CTX::instance()->lc) {
                      res[0] = n.x();
                      res[1] = n.y();
                      res[2] = n.z();
                      xm = pts[0].x();
                      ym = pts[0].y();
                      zm = pts[0].z();
                      ok = true;
                      break;
                  }
              }
          }
          if(ok) {
              double ex[3], t1[3], t2[3];
              ex[0] = ex[1] = ex[2] = 0.0;
              if(res[0] == 0.)
                  ex[0] = 1.0;
              else if(res[1] == 0.)
                  ex[1] = 1.0;
              else
                  ex[2] = 1.0;
              prodve(res, ex, t1);
              norme(t1);
              prodve(t1, res, t2);
              norme(t2);
              res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
              fillMeanPlane(res, t1, t2, meanPlane);
              return;
          }
      }
  
      std::vector<GVertex *> const &verts = vertices();
      std::vector<GVertex *>::const_iterator itv = verts.begin();
      for(; itv != verts.end(); itv++) {
          const GVertex *v = *itv;
          pts.push_back(SPoint3(v->x(), v->y(), v->z()));
      }
  
      bool colinear = (pts.size() < 3);
      if(pts.size() > 2) {
          SVector3 d01(pts[0], pts[1]), d02(pts[0], pts[2]);
          if(norm(crossprod(d01, d02)) < 1e-12) colinear = true;
      }
  
      if(colinear) {
          Msg::Debug("Adding edge points (%d) to compute mean plane of face %d",
                     pts.size(), tag());
          std::vector<GEdge *> const &edg = edges();
          for(std::vector<GEdge *>::const_iterator ite = edg.begin();
              ite != edg.end(); ite++) {
              const GEdge *e = *ite;
              if(e->mesh_vertices.size() > 1) {
                  for(std::size_t i = 0; i < e->mesh_vertices.size(); i++)
                      pts.push_back(e->mesh_vertices[i]->point());
              }
              else {
                  Range<double> b = e->parBounds(0);
                  GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low()));
                  pts.push_back(SPoint3(p1.x(), p1.y(), p1.z()));
                  GPoint p2 = e->point(b.low() + 0.666 * (b.high() - b.low()));
                  pts.push_back(SPoint3(p2.x(), p2.y(), p2.z()));
              }
          }
      }
  
      computeMeanPlane(pts);
  }
  
  void GFace::computeMeanPlane(const std::vector<MVertex *> &points)
  {
      std::vector<SPoint3> pts;
      for(std::size_t i = 0; i < points.size(); i++)
          pts.push_back(SPoint3(points[i]->x(), points[i]->y(), points[i]->z()));
      computeMeanPlane(pts);
  }
  
  void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
  {
      if(points.empty()) return;
  
      // The concept of a mean plane computed in the sense of least
      // squares is fine for plane surfaces(!), but not really the best
      // one for non-plane surfaces. Indeed, imagine a quarter of a circle
      // extruded along a very small height: the mean plane will be in the
      // plane of the circle... Until we implement something better, there
      // is a test in this routine that checks the coherence of the
      // computation for non-plane surfaces and reverts to a basic mean
      // plane approximatation if the check fails.
  
      double xm = 0., ym = 0., zm = 0.;
      int ndata = points.size();
      int na = 3;
      for(int i = 0; i < ndata; i++) {
          xm += points[i].x();
          ym += points[i].y();
          zm += points[i].z();
      }
      xm /= (double)ndata;
      ym /= (double)ndata;
      zm /= (double)ndata;
  
      fullMatrix<double> U(ndata, na), V(na, na);
      fullVector<double> sigma(na);
      for(int i = 0; i < ndata; i++) {
          U(i, 0) = points[i].x() - xm;
          U(i, 1) = points[i].y() - ym;
          U(i, 2) = points[i].z() - zm;
      }
      U.svd(V, sigma);
      double res[4], svd[3];
      svd[0] = sigma(0);
      svd[1] = sigma(1);
      svd[2] = sigma(2);
      int min;
      if(std::abs(svd[0]) < std::abs(svd[1]) && std::abs(svd[0]) < std::abs(svd[2]))
          min = 0;
      else if(std::abs(svd[1]) < std::abs(svd[0]) &&
              std::abs(svd[1]) < std::abs(svd[2]))
          min = 1;
      else
          min = 2;
      res[0] = V(0, min);
      res[1] = V(1, min);
      res[2] = V(2, min);
      norme(res);
  
      double ex[3], t1[3], t2[3];
  
      // check coherence of results for non-plane surfaces
      if(geomType() != Plane && geomType() != DiscreteSurface) {
          double res2[3], c[3], sinc, angplan;
          double eps = 1.e-3;
  
          GPoint v1 = point(0.5, 0.5);
          GPoint v2 = point(0.5 + eps, 0.5);
          GPoint v3 = point(0.5, 0.5 + eps);
          t1[0] = v2.x() - v1.x();
          t1[1] = v2.y() - v1.y();
          t1[2] = v2.z() - v1.z();
          t2[0] = v3.x() - v1.x();
          t2[1] = v3.y() - v1.y();
          t2[2] = v3.z() - v1.z();
          norme(t1);
          norme(t2);
          // prodve(t1, t2, res2);
          // WARNING: the rest of the code assumes res = t2 x t1, not t1 x t2 (WTF?)
          prodve(t2, t1, res2);
          norme(res2);
  
          prodve(res, res2, c);
          double const cosc = prosca(res, res2);
          sinc = std::sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
          angplan = myatan2(sinc, cosc);
          angplan = angle_02pi(angplan) * 180. / M_PI;
          if((angplan > 70 && angplan < 110) || (angplan > 260 && angplan < 280)) {
              Msg::Info("SVD failed (angle=%g): using rough algo...", angplan);
              res[0] = res2[0];
              res[1] = res2[1];
              res[2] = res2[2];
              goto end;
          }
      }
  
      ex[0] = ex[1] = ex[2] = 0.0;
      if(res[0] == 0.)
          ex[0] = 1.0;
      else if(res[1] == 0.)
          ex[1] = 1.0;
      else
          ex[2] = 1.0;
  
      prodve(res, ex, t1);
      norme(t1);
      prodve(t1, res, t2);
      norme(t2);
  
  end:
      res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
  
      fillMeanPlane(res, t1, t2, meanPlane);
  
      Msg::Debug("Surface: %d", tag());
      Msg::Debug("SVD    : %g,%g,%g (min=%d)", svd[0], svd[1], svd[2], min);
      Msg::Debug("Plane  : (%g x + %g y + %g z = %g)", meanPlane.a, meanPlane.b,
                 meanPlane.c, meanPlane.d);
      Msg::Debug("Normal : (%g , %g , %g )", meanPlane.a, meanPlane.b, meanPlane.c);
      Msg::Debug("t1     : (%g , %g , %g )", t1[0], t1[1], t1[2]);
      Msg::Debug("t2     : (%g , %g , %g )", t2[0], t2[1], t2[2]);
      Msg::Debug("pt     : (%g , %g , %g )", meanPlane.x, meanPlane.y, meanPlane.z);
  
      // check coherence for plane surfaces
      if(geomType() == Plane) {
          SBoundingBox3d bb = bounds();
          double lc = norm(SVector3(bb.max(), bb.min()));
          std::vector<GVertex *> const &verts = vertices();
          std::vector<GVertex *>::const_iterator itv = verts.begin();
          for(; itv != verts.end(); itv++) {
              const GVertex *v = *itv;
              double const d = meanPlane.a * v->x() + meanPlane.b * v->y() +
                      meanPlane.c * v->z() - meanPlane.d;
              if(std::abs(d) > lc * 1.e-3) {
                  Msg::Debug("Plane surface %d (%gx+%gy+%gz=%g) is not plane!", tag(),
                             meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d);
                  Msg::Debug("Control point %d = (%g,%g,%g), val=%g", v->tag(), v->x(),
                             v->y(), v->z(), d);
                  break;
              }
          }
      }
  }
  
  void GFace::getMeanPlaneData(double VX[3], double VY[3], double &x, double &y,
  double &z) const
  {
      VX[0] = meanPlane.plan[0][0];
      VX[1] = meanPlane.plan[0][1];
      VX[2] = meanPlane.plan[0][2];
      VY[0] = meanPlane.plan[1][0];
      VY[1] = meanPlane.plan[1][1];
      VY[2] = meanPlane.plan[1][2];
      x = meanPlane.x;
      y = meanPlane.y;
      z = meanPlane.z;
  }
  
  void GFace::getMeanPlaneData(double plan[3][3]) const
  {
      for(int i = 0; i < 3; i++)
          for(int j = 0; j < 3; j++) plan[i][j] = meanPlane.plan[i][j];
  }
  
  void GFace::computeMeshSizeFieldAccuracy(double &avg, double &max_e,
                                           double &min_e, int &nE, int &GS)
  {
  #if defined(HAVE_MESH)
      std::set<MEdge, Less_Edge> es;
      for(std::size_t i = 0; i < getNumMeshElements(); i++) {
          MElement *e = getMeshElement(i);
          for(int j = 0; j < e->getNumEdges(); j++) es.insert(e->getEdge(j));
      }
  
      avg = 0;
      min_e = 1.e22;
      max_e = 0;
      nE = es.size();
      GS = 0;
      double oneoversqr2 = 1. / sqrt(2.);
      double sqr2 = sqrt(2.);
      for(std::set<MEdge, Less_Edge>::const_iterator it = es.begin();
          it != es.end(); ++it) {
          double u1, v1, u2, v2;
          MVertex *vert1 = it->getVertex(0);
          vert1->getParameter(0, u1);
          vert1->getParameter(1, v1);
          MVertex *vert2 = it->getVertex(1);
          vert2->getParameter(0, u2);
          vert2->getParameter(1, v2);
          double l1 = BGM_MeshSize(this, u1, v1, vert1->x(), vert1->y(), vert1->z());
          double l2 = BGM_MeshSize(this, u2, v2, vert2->x(), vert2->y(), vert2->z());
          double correctLC = 0.5 * (l1 + l2);
          double lone = it->length() / correctLC;
          if(lone > oneoversqr2 && lone < sqr2) GS++;
          avg += lone > 1 ? (1. / lone) - 1. : lone - 1.;
          max_e = std::max(max_e, lone);
          min_e = std::min(min_e, lone);
      }
  #endif
  }
  
  double GFace::curvatureDiv(const SPoint2 &param) const
  {
      if(geomType() == Plane) return 0;
  
      // X=X(u,v) Y=Y(u,v) Z=Z(u,v)
      // curv = div n = dnx/dx + dny/dy + dnz/dz
      // dnx/dx = dnx/du du/dx + dnx/dv dv/dx
  
      const double eps = 1.e-5;
  
      Pair<SVector3, SVector3> der = firstDer(param);
  
      SVector3 du = der.first();
      SVector3 dv = der.second();
      SVector3 nml = crossprod(du, dv);
  
      double detJ = norm(nml);
  
      du.normalize();
      dv.normalize();
  
      SVector3 n1, n2, n3, n4;
      if(param.x() - eps < 0.0) {
          n1 = normal(SPoint2(param.x(), param.y()));
          n2 = normal(SPoint2(param.x() + eps, param.y()));
      }
      else {
          n1 = normal(SPoint2(param.x() - eps, param.y()));
          n2 = normal(SPoint2(param.x(), param.y()));
      }
      if(param.y() - eps < 0.0) {
          n3 = normal(SPoint2(param.x(), param.y()));
          n4 = normal(SPoint2(param.x(), param.y() + eps));
      }
      else {
          n3 = normal(SPoint2(param.x(), param.y() - eps));
          n4 = normal(SPoint2(param.x(), param.y()));
      }
  
      SVector3 dndu = 100000 * (n2 - n1);
      SVector3 dndv = 100000 * (n4 - n3);
  
      SVector3 dudu = SVector3();
      SVector3 dvdv = SVector3();
      SVector3 dudv = SVector3();
      secondDer(param, dudu, dvdv, dudv);
  
      double ddu = dot(dndu, du);
      double ddv = dot(dndv, dv);
  
      return (fabs(ddu) + fabs(ddv)) / detJ;
  }
  
  double GFace::curvatureMax(const SPoint2 &param) const
  {
      if(geomType() == Plane) return 0.;
  
      double eigVal[2], eigVec[8];
      getMetricEigenVectors(param, eigVal, eigVec);
  
      return fabs(eigVal[1]);
  }
  
  double GFace::curvatures(const SPoint2 &param, SVector3 &dirMax,
                           SVector3 &dirMin, double &curvMax,
                           double &curvMin) const
  {
      Pair<SVector3, SVector3> D1 = firstDer(param);
  
      if(geomType() == Plane) {
          dirMax = D1.first();
          dirMin = D1.second();
          curvMax = 0.;
          curvMin = 0.;
          return 0.;
      }
  
      if(geomType() == Sphere) {
          dirMax = D1.first();
          dirMin = D1.second();
          curvMax = curvatureDiv(param);
          curvMin = curvMax;
          return curvMax;
      }
  
      double eigVal[2], eigVec[8];
      getMetricEigenVectors(param, eigVal, eigVec);
  
      // curvatures and main directions
      curvMax = fabs(eigVal[1]);
      curvMin = fabs(eigVal[0]);
      dirMax = eigVec[1] * D1.first() + eigVec[3] * D1.second();
      dirMin = eigVec[0] * D1.first() + eigVec[2] * D1.second();
  
      return curvMax;
  }
  
  double GFace::getMetricEigenvalue(const SPoint2 &)
  {
      Msg::Error("Metric eigenvalue is not implemented for this type of surface");
      return 0.;
  }
  
  // eigen values are absolute values and sorted from min to max of absolute
  // values eigen vectors are the corresponding COLUMNS of eigVec
  void GFace::getMetricEigenVectors(const SPoint2 &param, double eigVal[2],
  double eigVec[4]) const
  {
      // first derivatives
      Pair<SVector3, SVector3> D1 = firstDer(param);
      SVector3 du = D1.first();
      SVector3 dv = D1.second();
      SVector3 nor = crossprod(du, dv);
      nor.normalize();
  
      // second derivatives
      SVector3 dudu = SVector3();
      SVector3 dvdv = SVector3();
      SVector3 dudv = SVector3();
      secondDer(param, dudu, dvdv, dudv);
  
      // first form
      double form1[2][2];
      form1[0][0] = normSq(du);
      form1[1][1] = normSq(dv);
      form1[0][1] = form1[1][0] = dot(du, dv);
  
      // second form
      double form2[2][2];
      form2[0][0] = dot(nor, dudu);
      form2[1][1] = dot(nor, dvdv);
      form2[0][1] = form2[1][0] = dot(nor, dudv);
  
      // inverse of first form
      double inv_form1[2][2];
      double inv_det_form1 =
              1. / (form1[0][0] * form1[1][1] - form1[1][0] * form1[0][1]);
      inv_form1[0][0] = inv_det_form1 * form1[1][1];
      inv_form1[1][1] = inv_det_form1 * form1[0][0];
      inv_form1[1][0] = inv_form1[0][1] = -1 * inv_det_form1 * form1[0][1];
  
      // N = (inverse of form1) X (form2)
      fullMatrix<double> N(2, 2);
      N(0, 0) = inv_form1[0][0] * form2[0][0] + inv_form1[0][1] * form2[1][0];
      N(0, 1) = inv_form1[0][0] * form2[0][1] + inv_form1[0][1] * form2[1][1];
      N(1, 0) = inv_form1[1][0] * form2[0][0] + inv_form1[1][1] * form2[1][0];
      N(1, 1) = inv_form1[1][0] * form2[0][1] + inv_form1[1][1] * form2[1][1];
  
      // eigen values and vectors of N
      fullMatrix<double> vl(2, 2), vr(2, 2);
      fullVector<double> dr(2), di(2);
      if(N.eig(dr, di, vl, vr, true)) {
          eigVal[0] = fabs(dr(0));
          eigVal[1] = fabs(dr(1));
          eigVec[0] = vr(0, 0);
          eigVec[2] = vr(1, 0);
          eigVec[1] = vr(0, 1);
          eigVec[3] = vr(1, 1);
      }
      else {
          Msg::Error("Problem in eigen vectors computation");
          Msg::Error(" N = [ %f %f ]", N(0, 0), N(0, 1));
          Msg::Error("     [ %f %f ]", N(1, 0), N(1, 1));
          for(int i = 0; i < 2; i++) eigVal[i] = 0.;
          for(int i = 0; i < 4; i++) eigVec[i] = 0.;
      }
      if(fabs(di(0)) > 1.e-12 || fabs(di(1)) > 1.e-12) {
          Msg::Error("Found imaginary eigenvalues");
      }
  }
  
  void GFace::XYZtoUV(double X, double Y, double Z, double &U, double &V,
                      double relax, bool onSurface) const
  {
      const double Precision = onSurface ? 1.e-8 : 1.e-3;
      const int MaxIter = onSurface ? 25 : 10;
      const int NumInitGuess = 9;
  
      double Unew = 0., Vnew = 0., err, err2;
      int iter;
      double mat[3][3], jac[3][3];
      double umin, umax, vmin, vmax;
      // don't use 0.9, 0.1 it fails with ruled surfaces
      double initu[NumInitGuess] = {0.5, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 1.0, 0.0};
      double initv[NumInitGuess] = {0.5, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 1.0, 0.0};
  
      Range<double> ru = parBounds(0);
      Range<double> rv = parBounds(1);
      umin = ru.low();
      umax = ru.high();
      vmin = rv.low();
      vmax = rv.high();
  
      const double tol = Precision * (SQU(umax - umin) + SQU(vmax - vmin));
      for(int i = 0; i < NumInitGuess; i++) {
          initu[i] = umin + initu[i] * (umax - umin);
          initv[i] = vmin + initv[i] * (vmax - vmin);
      }
  
      for(int i = 0; i < NumInitGuess; i++) {
          for(int j = 0; j < NumInitGuess; j++) {
              U = initu[i];
              V = initv[j];
              err = 1.0;
              iter = 1;
  
              GPoint P = point(U, V);
              err2 = sqrt(SQU(X - P.x()) + SQU(Y - P.y()) + SQU(Z - P.z()));
              if(err2 < 1.e-8 * CTX::instance()->lc) return;
  
              while(err > tol && iter < MaxIter) {
                  P = point(U, V);
                  Pair<SVector3, SVector3> der = firstDer(SPoint2(U, V));
                  mat[0][0] = der.left().x();
                  mat[0][1] = der.left().y();
                  mat[0][2] = der.left().z();
                  mat[1][0] = der.right().x();
                  mat[1][1] = der.right().y();
                  mat[1][2] = der.right().z();
                  mat[2][0] = 0.;
                  mat[2][1] = 0.;
                  mat[2][2] = 0.;
                  invert_singular_matrix3x3(mat, jac);
  
                  Unew = U + relax * (jac[0][0] * (X - P.x()) + jac[1][0] * (Y - P.y()) +
                          jac[2][0] * (Z - P.z()));
                  Vnew = V + relax * (jac[0][1] * (X - P.x()) + jac[1][1] * (Y - P.y()) +
                          jac[2][1] * (Z - P.z()));
  
                  // don't remove this test: it is important
                  if((Unew > umax + tol || Unew < umin - tol) &&
                          (Vnew > vmax + tol || Vnew < vmin - tol))
                      break;
  
                  err = SQU(Unew - U) + SQU(Vnew - V);
                  err2 = sqrt(SQU(X - P.x()) + SQU(Y - P.y()) + SQU(Z - P.z()));
  
                  iter++;
                  U = Unew;
                  V = Vnew;
              }
  
              if(iter < MaxIter && err <= tol && Unew <= umax && Vnew <= vmax &&
                      Unew >= umin && Vnew >= vmin) {
                  if(onSurface && err2 > 1.e-4 * CTX::instance()->lc &&
                          !CTX::instance()->mesh.NewtonConvergenceTestXYZ) {
                      Msg::Warning("Converged for i=%d j=%d (err=%g iter=%d) BUT "
                                   "xyz error = %g in point (%e,%e,%e) on surface %d",
                                   i, j, err, iter, err2, X, Y, Z, tag());
                  }
  
                  if(onSurface && err2 > 1.e-4 * CTX::instance()->lc &&
                          CTX::instance()->mesh.NewtonConvergenceTestXYZ) {
                      // not converged in XYZ coordinates
                  }
                  else {
                      return;
                  }
              }
          }
      }
  
      if(!onSurface) return;
  
      if(relax < 1.e-6)
          Msg::Error("Could not converge: surface mesh will be wrong");
      else {
          Msg::Info("point %g %g %g : Relaxation factor = %g", X, Y, Z, 0.75 * relax);
          XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
      }
  }
  
  SPoint2 GFace::parFromPoint(const SPoint3 &p, bool onSurface) const
  {
      double U = 0., V = 0.;
      XYZtoUV(p.x(), p.y(), p.z(), U, V, 1.0, onSurface);
      return SPoint2(U, V);
  }
  
  #if defined(HAVE_ALGLIB)
  
  class data_wrapper {
  private:
      const GFace *gf;
      SPoint3 point;
  
  public:
      data_wrapper()
      {
          gf = nullptr;
          point = SPoint3();
      }
      ~data_wrapper() {}
      const GFace *get_face() { return gf; }
      void set_face(const GFace *face) { gf = face; }
      SPoint3 get_point() { return point; }
      void set_point(const SPoint3 &_point) { point = SPoint3(_point); }
  };
  
  // Callback function for ALGLIB
  void bfgs_callback(const alglib::real_1d_array &x, double &func,
                     alglib::real_1d_array &grad, void *ptr)
  {
      data_wrapper *w = static_cast<data_wrapper *>(ptr);
      SPoint3 p = w->get_point();
      const GFace *gf = w->get_face();
  
      // Value of the objective
      GPoint pnt = gf->point(x[0], x[1]);
      func = 0.5 * (p.x() - pnt.x()) * (p.x() - pnt.x()) +
              (p.y() - pnt.y()) * (p.y() - pnt.y()) +
              (p.z() - pnt.z()) * (p.z() - pnt.z());
      // printf("func : %f\n", func);
  
      // Value of the gradient
      Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(x[0], x[1]));
      grad[0] = -(p.x() - pnt.x()) * der.left().x() -
              (p.y() - pnt.y()) * der.left().y() -
              (p.z() - pnt.z()) * der.left().z();
      grad[1] = -(p.x() - pnt.x()) * der.right().x() -
              (p.y() - pnt.y()) * der.right().y() -
              (p.z() - pnt.z()) * der.right().z();
      //  printf("func %22.15E Gradients %22.15E %22.15E der %g %g %g\n",
      //         func, grad[0],
      //         grad[1],der.left().x(),der.left().y(),der.left().z());
  }
  #endif
  
  GPoint GFace::closestPoint(const SPoint3 &queryPoint,
                             const double initialGuess[2]) const
  {
  #if defined(HAVE_ALGLIB)
      // Test initial guess
      double min_u = initialGuess[0];
      double min_v = initialGuess[1];
      GPoint pnt = point(min_u, min_v);
      SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
      double min_dist = queryPoint.distance(spnt);
  
      // Try to find a better initial guess by sampling full parameter range
      const double nGuesses = 10.;
      const Range<double> uu = parBounds(0);
      const Range<double> vv = parBounds(1);
      const double ru = uu.high() - uu.low(), rv = vv.high() - vv.low();
      const double epsU = 1e-5 * ru, epsV = 1e-5 * rv;
      const double du = ru / nGuesses, dv = rv / nGuesses;
      for(double u = uu.low(); u <= uu.high() + epsU; u += du) {
          for(double v = vv.low(); v <= vv.high() + epsV; v += dv) {
              GPoint pnt = point(u, v);
              SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
              double dist = queryPoint.distance(spnt);
              if(dist < min_dist) {
                  min_dist = dist;
                  min_u = u;
                  min_v = v;
              }
          }
      }
  
      // Set up optimisation problem
      alglib::ae_int_t dim = 2;
      alglib::ae_int_t corr = 2; // Num of corrections in the scheme in [3,7]
      alglib::minlbfgsstate state;
      alglib::real_1d_array x;
      const double initialCond[2] = {min_u, min_v};
      x.setcontent(dim, initialCond);
      minlbfgscreate(2, corr, x, state);
  
      // Set stopping criteria
      const double epsg = 1.e-12;
      const double epsf = 0.;
      const double epsx = 0.;
      const alglib::ae_int_t maxits = 500;
      minlbfgssetcond(state, epsg, epsf, epsx, maxits);
  
      // Solve problem
      data_wrapper w;
      w.set_point(queryPoint);
      w.set_face(this);
      minlbfgsoptimize(state, bfgs_callback, nullptr, &w);
  
      // Get results
      alglib::minlbfgsreport rep;
      minlbfgsresults(state, x, rep);
      GPoint pntF = point(x[0], x[1]);
      return pntF;
  
  #else
      Msg::Error("Closest point not implemented for this type of surface");
      SPoint2 p = parFromPoint(queryPoint, false);
      return point(p);
  #endif
  }
  
  bool GFace::containsParam(const SPoint2 &pt)
  {
      Range<double> uu = parBounds(0);
      Range<double> vv = parBounds(1);
      if((pt.x() >= uu.low() && pt.x() <= uu.high()) &&
              (pt.y() >= vv.low() && pt.y() <= vv.high()))
          return true;
      else
          return false;
  }
  
  SVector3 GFace::normal(const SPoint2 &param) const
  {
      Pair<SVector3, SVector3> der = firstDer(param);
      SVector3 n = crossprod(der.first(), der.second());
      n.normalize();
      return n;
  }
  
  bool GFace::buildRepresentationCross(bool force)
  {
      if(cross[0].size()) {
          if(force) {
              cross[0].clear();
              cross[0].clear();
          }
          else
              return true;
      }
  
      Range<double> ubounds = parBounds(0);
      Range<double> vbounds = parBounds(1);
      // try to compute something better for Gmsh surfaces
      if(getNativeType() == GmshModel && geomType() == Plane) {
          SBoundingBox3d bb;
          std::vector<GEdge *> ed = edges();
          for(std::vector<GEdge *>::const_iterator it = ed.begin(); it != ed.end();
              it++) {
              GEdge *ge = *it;
              if(ge->geomType() != DiscreteCurve &&
                      ge->geomType() != BoundaryLayerCurve) {
                  Range<double> t_bounds = ge->parBounds(0);
                  const int N = 5;
                  double t0 = t_bounds.low(), dt = t_bounds.high() - t_bounds.low();
                  for(int i = 0; i < N; i++) {
                      double t = t0 + dt / (double)(N - 1) * i;
                      GPoint p = ge->point(t);
                      SPoint2 uv = parFromPoint(SPoint3(p.x(), p.y(), p.z()));
                      bb += SPoint3(uv.x(), uv.y(), 0.);
                  }
              }
          }
          if(!bb.empty()) {
              ubounds = Range<double>(bb.min().x(), bb.max().x());
              vbounds = Range<double>(bb.min().y(), bb.max().y());
          }
      }
  
      bool tri = (geomType() == RuledSurface && edges().size() == 3);
      if(CTX::instance()->geom.oldRuledSurface) tri = false;
      double c = tri ? 0.75 : 0.5;
      double uav = c * (ubounds.high() + ubounds.low());
      double vav = (1 - c) * (vbounds.high() + vbounds.low());
      double u2 = 0.5 * (ubounds.high() + ubounds.low());
      double v2 = 0.5 * (vbounds.high() + vbounds.low());
      double ud = (ubounds.high() - ubounds.low());
      double vd = (vbounds.high() - vbounds.low());
      const int N = 100;
      for(int dir = 0; dir < 2; dir++) {
          cross[dir].push_back(std::vector<SPoint3>());
          for(int i = 0; i < N; i++) {
              double t = (double)i / (double)(N - 1);
              SPoint2 uv;
              if(!dir) {
                  if(tri)
                      uv.setPosition(u2 + u2 * t, vbounds.low() + v2 * t);
                  else
                      uv.setPosition(ubounds.low() + ud * t, vav);
              }
              else {
                  if(tri)
                      uv.setPosition(u2 + u2 * t, v2 - v2 * t);
                  else
                      uv.setPosition(uav, vbounds.low() + vd * t);
              }
              GPoint p = point(uv);
              SPoint3 pt(p.x(), p.y(), p.z());
              bool inside =
                      (geomType() == Plane) ? containsPoint(pt) : containsParam(uv);
              if(inside) { cross[dir].back().push_back(pt); }
              else {
                  if(cross[dir].back().size())
                      cross[dir].push_back(std::vector<SPoint3>());
              }
          }
      }
      // if we couldn't determine a cross, add a dummy one so that we won't try
      // again unless we force the recomputation
      if(cross[0].empty()) {
          cross[0].push_back(std::vector<SPoint3>());
          return false;
      }
      return true;
  }
  
  bool GFace::buildSTLTriangulation(bool force)
  {
      if(stl_triangles.size() && !force) return true;
      stl_vertices_uv.clear();
      stl_vertices_xyz.clear();
      stl_triangles.clear();
  
      // Build a simple triangulation for surfaces which we know are not trimmed
      if(geomType() == ParametricSurface) {
          const int nu = 64, nv = 64;
          Range<double> ubounds = parBounds(0);
          Range<double> vbounds = parBounds(1);
          double umin = ubounds.low(), umax = ubounds.high();
          double vmin = vbounds.low(), vmax = vbounds.high();
          for(int i = 0; i < nu; i++) {
              for(int j = 0; j < nv; j++) {
                  double u = umin + (double)i / (double)(nu - 1) * (umax - umin);
                  double v = vmin + (double)j / (double)(nv - 1) * (vmax - vmin);
                  stl_vertices_uv.push_back(SPoint2(u, v));
                  GPoint gp = point(u, v);
                  stl_vertices_xyz.push_back(SPoint3(gp.x(), gp.y(), gp.z()));
              }
          }
          for(int i = 0; i < nu - 1; i++) {
              for(int j = 0; j < nv - 1; j++) {
                  stl_triangles.push_back(i * nv + j);
                  stl_triangles.push_back((i + 1) * nv + j);
                  stl_triangles.push_back((i + 1) * nv + j + 1);
                  stl_triangles.push_back(i * nv + j);
                  stl_triangles.push_back((i + 1) * nv + j + 1);
                  stl_triangles.push_back(i * nv + j + 1);
              }
          }
          return true;
      }
  
      return false;
  }
  
  bool GFace::fillVertexArray(bool force)
  {
      if(va_geom_triangles) {
          if(force)
              delete va_geom_triangles;
          else
              return true;
      }
  
      if(!buildSTLTriangulation(force)) return false;
      if(stl_triangles.empty()) return false;
  
      va_geom_triangles = new VertexArray(3, stl_triangles.size() / 3);
      unsigned int c = useColor() ? getColor() : CTX::instance()->color.geom.surface;
      unsigned int col[4] = {c, c, c, c};
      if(stl_vertices_xyz.size() &&
              (stl_vertices_xyz.size() == stl_normals.size())) {
          for(std::size_t i = 0; i < stl_triangles.size(); i += 3) {
              SPoint3 &p1(stl_vertices_xyz[stl_triangles[i]]);
              SPoint3 &p2(stl_vertices_xyz[stl_triangles[i + 1]]);
              SPoint3 &p3(stl_vertices_xyz[stl_triangles[i + 2]]);
              double x[3] = {p1.x(), p2.x(), p3.x()};
              double y[3] = {p1.y(), p2.y(), p3.y()};
              double z[3] = {p1.z(), p2.z(), p3.z()};
              SVector3 n[3] = {stl_normals[stl_triangles[i]],
                               stl_normals[stl_triangles[i + 1]],
                               stl_normals[stl_triangles[i + 2]]};
              va_geom_triangles->add(x, y, z, n, col);
          }
      }
      else if(stl_vertices_uv.size()) {
          for(std::size_t i = 0; i < stl_triangles.size(); i += 3) {
              SPoint2 &p1(stl_vertices_uv[stl_triangles[i]]);
              SPoint2 &p2(stl_vertices_uv[stl_triangles[i + 1]]);
              SPoint2 &p3(stl_vertices_uv[stl_triangles[i + 2]]);
              GPoint gp1 = point(p1);
              GPoint gp2 = point(p2);
              GPoint gp3 = point(p3);
              double x[3] = {gp1.x(), gp2.x(), gp3.x()};
              double y[3] = {gp1.y(), gp2.y(), gp3.y()};
              double z[3] = {gp1.z(), gp2.z(), gp3.z()};
              SVector3 n[3] = {normal(p1), normal(p2), normal(p3)};
              va_geom_triangles->add(x, y, z, n, col);
          }
      }
      va_geom_triangles->finalize();
      return true;
  }
  
  int GFace::genusGeom() const
  {
      int nSeams = 0;
      std::set<GEdge *> single_seams;
      for(std::vector<GEdge *>::const_iterator it = l_edges.begin();
          it != l_edges.end(); ++it) {
          if((*it)->isSeam(this)) {
              nSeams++;
              std::set<GEdge *>::iterator it2 = single_seams.find(*it);
              if(it2 != single_seams.end())
                  single_seams.erase(it2);
              else
                  single_seams.insert(*it);
          }
      }
      return nSeams - single_seams.size();
  }
  
  bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
                             std::vector<SPoint2> *uvpoints,
                             std::vector<SVector3> *normals)
  {
      if(!points) return false;
  
      if(buildSTLTriangulation() && stl_vertices_uv.size()) {
          for(std::size_t i = 0; i < stl_triangles.size(); i += 3) {
              SPoint2 &p0(stl_vertices_uv[stl_triangles[i]]);
              SPoint2 &p1(stl_vertices_uv[stl_triangles[i + 1]]);
              SPoint2 &p2(stl_vertices_uv[stl_triangles[i + 2]]);
              GPoint gp0 = point(p0);
              GPoint gp1 = point(p1);
              GPoint gp2 = point(p2);
              double maxEdge = std::max(gp0.distance(gp1),
                                        std::max(gp1.distance(gp2), gp2.distance(gp0)));
              int N = (int)(maxEdge / maxDist);
              for(double u = 0.; u < 1.; u += 1. / N) {
                  for(double v = 0.; v < 1 - u; v += 1. / N) {
                      SPoint2 p = p0 * (1. - u - v) + p1 * u + p2 * v;
                      GPoint gp(point(p));
                      points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
                      if(uvpoints) uvpoints->push_back(p);
                      if(normals) normals->push_back(normal(p));
                  }
              }
          }
      }
      else {
          int N = 1000; //(int)(maxDX / maxDist);
          Range<double> b1 = parBounds(0);
          Range<double> b2 = parBounds(1);
          for(int i = 0; i < N; i++) {
              for(int j = 0; j < N; j++) {
                  double u = (double)i / (N - 1);
                  double v = (double)j / (N - 1);
                  double t1 = b1.low() + u * (b1.high() - b1.low());
                  double t2 = b2.low() + v * (b2.high() - b2.low());
                  GPoint gp = point(t1, t2);
                  points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
                  if(uvpoints) uvpoints->push_back(SPoint2(t1, t2));
                  if(normals) normals->push_back(normal(SPoint2(t1, t2)));
              }
          }
      }
      return true;
  }
  
  #if defined(HAVE_MESH)
  
  static void meshCompound(GFace *gf, bool verbose)
  {
      discreteFace *df = new discreteFace(gf->model(), gf->tag() + 100000);
  
      std::vector<GFace *> triangles_tag;
      std::vector<SPoint2> triangles_uv;
  
      for(std::size_t i = 0; i < gf->_compound.size(); i++) {
          GFace *c = (GFace *)gf->_compound[i];
          df->triangles.insert(df->triangles.end(), c->triangles.begin(),
                               c->triangles.end());
          df->mesh_vertices.insert(df->mesh_vertices.end(), c->mesh_vertices.begin(),
                                   c->mesh_vertices.end());
          for(std::size_t j = 0; j < c->triangles.size(); j++) {
              triangles_tag.push_back(c);
              for(int k = 0; k < 3; k++) {
                  SPoint2 param;
                  reparamMeshVertexOnFace(c->triangles[j]->getVertex(k), c, param);
                  triangles_uv.push_back(param);
              }
          }
          c->triangles.clear();
          c->mesh_vertices.clear();
      }
  
      df->createGeometry();
      df->mesh(verbose);
  
      for(std::size_t i = 0; i < df->mesh_vertices.size(); i++) {
          double u, v;
          df->mesh_vertices[i]->getParameter(0, u);
          df->mesh_vertices[i]->getParameter(1, v);
          double U, V;
          int position = df->trianglePosition(u, v, U, V);
          if(position != -1) {
              triangles_tag[position]->mesh_vertices.push_back(df->mesh_vertices[i]);
              df->mesh_vertices[i]->setEntity(triangles_tag[position]);
              if(0 && triangles_tag[position]->geomType() != GEntity::DiscreteSurface) {
                  SPoint2 p0 = triangles_uv[3 * position + 0];
                  SPoint2 p1 = triangles_uv[3 * position + 1];
                  SPoint2 p2 = triangles_uv[3 * position + 2];
                  SPoint2 p = p0 * (1 - U - V) + p1 * U + p2 * V;
                  GPoint gp = triangles_tag[position]->point(p);
                  df->mesh_vertices[i]->setParameter(0, p.x());
                  df->mesh_vertices[i]->setParameter(1, p.y());
                  df->mesh_vertices[i]->x() = gp.x();
                  df->mesh_vertices[i]->y() = gp.y();
                  df->mesh_vertices[i]->z() = gp.z();
              }
          }
          else {
              df->mesh_vertices.push_back(df->mesh_vertices[i]);
              df->mesh_vertices[i]->setEntity(gf);
          }
      }
  
      for(std::size_t i = 0; i < df->triangles.size(); i++) {
          MTriangle *t = df->triangles[i];
          if(t->getVertex(0)->onWhat()->dim() == 2)
              ((GFace *)t->getVertex(0)->onWhat())->triangles.push_back(t);
          else if(t->getVertex(1)->onWhat()->dim() == 2)
              ((GFace *)t->getVertex(1)->onWhat())->triangles.push_back(t);
          else if(t->getVertex(2)->onWhat()->dim() == 2)
              ((GFace *)t->getVertex(2)->onWhat())->triangles.push_back(t);
          else
              gf->triangles.push_back(t); // FIXME could be better!
      }
      // gf->triangles = df->triangles;
      df->triangles.clear();
      df->mesh_vertices.clear();
      delete df;
  }
  #endif
  
  void GFace::mesh(bool verbose)
  {
      if(CTX::instance()->debugSurface > 0 &&
              tag() != CTX::instance()->debugSurface) {
          meshStatistics.status = GFace::DONE;
          return;
      }
  
  #if defined(HAVE_MESH)
      meshGFace mesher;
      mesher(this, verbose);
      if(_compound.size()) { // Some faces are meshed together
          if(_compound[0] == this) { // I'm the one that makes the compound job
              bool ok = true;
              for(std::size_t i = 0; i < _compound.size(); i++) {
                  GFace *gf = (GFace *)_compound[i];
                  ok &= (gf->meshStatistics.status == GFace::DONE);
              }
              if(!ok) { meshStatistics.status = GFace::PENDING; }
              else {
                  meshCompound(this, verbose);
                  meshStatistics.status = GFace::DONE;
                  return;
              }
          }
      }
  #endif
  }
  
  void GFace::moveToValidRange(SPoint2 &pt) const
  {
      for(int i = 0; i < 2; i++) {
          if(periodic(i)) {
              Range<double> range = parBounds(i);
              double tol = 1e-6 * (range.high() - range.low());
              if(pt[i] < range.low() - tol) pt[i] += period(i);
              if(pt[i] > range.high() + tol) pt[i] -= period(i);
              if(pt[i] < range.low()) pt[i] = range.low();
              if(pt[i] > range.high()) pt[i] = range.high();
          }
      }
  }
  
  void GFace::relocateMeshVertices()
  {
      for(std::size_t i = 0; i < mesh_vertices.size(); i++) {
          MVertex *v = mesh_vertices[i];
          double u0 = 0., u1 = 0.;
          if(v->getParameter(0, u0) && v->getParameter(1, u1)) {
              GPoint p = point(u0, u1);
              v->x() = p.x();
              v->y() = p.y();
              v->z() = p.z();
          }
      }
  }
  
  void GFace::setMeshMaster(GFace *master, const std::vector<double> &tfo)
  {
      std::vector<GEdge *>::const_iterator eIter;
      std::list<GVertex *>::iterator vIter;
  
      Msg::Info("Setting mesh master using transformation");
  
      // list all vertices and construct vertex to edge correspondence for local
      // edge
  
      std::set<GVertex *> l_vertices;
      std::map<std::pair<GVertex *, GVertex *>, GEdge *> l_vtxToEdge;
  
      for(eIter = l_edges.begin(); eIter != l_edges.end(); ++eIter) {
          GVertex *v0 = (*eIter)->getBeginVertex();
          GVertex *v1 = (*eIter)->getEndVertex();
          if(v0 && v1){
              l_vertices.insert(v0);
              l_vertices.insert(v1);
              l_vtxToEdge[std::make_pair(v0, v1)] = (*eIter);
          }
      }
  
      for(eIter = embedded_edges.begin(); eIter != embedded_edges.end(); ++eIter) {
          GVertex *v0 = (*eIter)->getBeginVertex();
          GVertex *v1 = (*eIter)->getEndVertex();
          if(v0 && v1){
              l_vertices.insert(v0);
              l_vertices.insert(v1);
              l_vtxToEdge[std::make_pair(v0, v1)] = (*eIter);
          }
      }
  
      l_vertices.insert(embedded_vertices.begin(), embedded_vertices.end());
  
      // list all vertices and vertex to edge correspondence for remote edge
  
      std::vector<GEdge *> const &m_edges = master->edges();
      std::set<GVertex *> m_vertices;
      std::map<std::pair<GVertex *, GVertex *>, GEdge *> m_vtxToEdge;
      for(eIter = m_edges.begin(); eIter != m_edges.end(); ++eIter) {
          GVertex *v0 = (*eIter)->getBeginVertex();
          GVertex *v1 = (*eIter)->getEndVertex();
          if(v0 && v1){
              m_vertices.insert(v0);
              m_vertices.insert(v1);
              m_vtxToEdge[std::make_pair(v0, v1)] = (*eIter);
          }
      }
  
      std::vector<GEdge *> const &m_embedded_edges = master->embeddedEdges();
  
      for(eIter = m_embedded_edges.begin(); eIter != m_embedded_edges.end();
          eIter++) {
          GVertex *v0 = (*eIter)->getBeginVertex();
          GVertex *v1 = (*eIter)->getEndVertex();
          if(v0 && v1){
              m_vertices.insert(v0);
              m_vertices.insert(v1);
              m_vtxToEdge[std::make_pair(v0, v1)] = (*eIter);
          }
      }
  
      std::set<GVertex *, GEntityLessThan> m_embedded_vertices = master->embeddedVertices();
      m_vertices.insert(m_embedded_vertices.begin(), m_embedded_vertices.end());
  
      // check topological correspondence
  
      if(l_vertices.size() != m_vertices.size()) {
          Msg::Error(
                      "Periodic connection specified between topologically "
                      "incompatible surfaces %d and %d (that have %d vs %d model vertices)",
                      master->tag(), tag(), l_vertices.size(), m_vertices.size());
          return;
      }
  
      if(l_vtxToEdge.size() != m_vtxToEdge.size()) {
          Msg::Error(
                      "Periodic connection specified between topologically "
                      "incompatible surfaces %d and %d (that have %d vs %d model edges)",
                      master->tag(), tag(), l_vtxToEdge.size(), m_vtxToEdge.size());
          return;
      }
  
      // compute corresponding vertices
  
      std::map<GVertex *, GVertex *> gVertexCounterparts;
  
      std::set<GVertex *>::iterator mvIter;
      for(mvIter = m_vertices.begin(); mvIter != m_vertices.end(); ++mvIter) {
          GVertex *m_vertex = *mvIter;
  
          SPoint3 xyzOri((*mvIter)->x(), (*mvIter)->y(), (*mvIter)->z());
          SPoint3 xyzTfo(0, 0, 0);
  
          int idx = 0;
          for(int i = 0; i < 3; i++) {
              for(int j = 0; j < 3; j++) xyzTfo[i] += xyzOri[j] * tfo[idx++];
              xyzTfo[i] += tfo[idx++];
          }
  
          GVertex *l_vertex = nullptr;
  
          double dist_min = 1.e22;
          std::set<GVertex *>::iterator lvIter = l_vertices.begin();
          for(; lvIter != l_vertices.end(); ++lvIter) {
              SPoint3 xyz((*lvIter)->x(), (*lvIter)->y(), (*lvIter)->z());
              SVector3 dist = xyz - xyzTfo;
              dist_min = std::min(dist_min, dist.norm());
              if(dist.norm() < CTX::instance()->geom.tolerance * CTX::instance()->lc) {
                  l_vertex = *lvIter;
                  break;
              }
          }
  
          if(l_vertex == nullptr) {
              Msg::Error("Was not able to find corresponding node %d "
                         "for periodic connection of surface %d to %d "
                         "(min distance is %g with a tolerance of %g)",
                         m_vertex->tag(), master->tag(), tag(), dist_min,
                         CTX::instance()->geom.tolerance * CTX::instance()->lc);
              return;
          }
          gVertexCounterparts[l_vertex] = m_vertex;
      }
  
      if(gVertexCounterparts.size() != m_vertices.size()) {
          Msg::Error("Could not find all node correspondances "
                     "for the periodic connection from surface %d to %d",
                     master->tag(), tag());
          return;
      }
  
      // construct edge correspondence and update the edge masters
  
      std::map<GEdge *, std::pair<GEdge *, int> > gEdgeCounterparts;
  
      std::map<std::pair<GVertex *, GVertex *>, GEdge *>::iterator lv2eIter;
      for(lv2eIter = l_vtxToEdge.begin(); lv2eIter != l_vtxToEdge.end();
          lv2eIter++) {
          std::pair<GVertex *, GVertex *> lPair = lv2eIter->first;
          GEdge *localEdge = lv2eIter->second;
  
          std::pair<GVertex *, GVertex *> mPair(gVertexCounterparts[lPair.first],
                  gVertexCounterparts[lPair.second]);
          int sign = 1;
          std::map<std::pair<GVertex *, GVertex *>, GEdge *>::iterator mv2eIter =
                  m_vtxToEdge.find(mPair);
          if(mv2eIter == m_vtxToEdge.end()) {
              sign *= -1;
              std::pair<GVertex *, GVertex *> backward(mPair.second, mPair.first);
              mv2eIter = m_vtxToEdge.find(backward);
          }
  
          if(mv2eIter == m_vtxToEdge.end()) {
              Msg::Error("Could not find periodic copy of edge %d-%d "
                         "(corresponding to vertices %d %d) in face %d",
                         lPair.first->tag(), lPair.second->tag(), mPair.first->tag(),
                         mPair.second->tag(), master->tag());
              return;
          }
          GEdge *masterEdge = mv2eIter->second;
  
          if(masterEdge->getMeshMaster() != localEdge) {
              localEdge->setMeshMaster(masterEdge, tfo);
              Msg::Info("Setting edge master %d - %d", localEdge->tag(),
                        masterEdge->tag());
          }
          gEdgeCounterparts[localEdge] = std::make_pair(masterEdge, sign);
      }
  
      // complete the information at the edge level
      edgeCounterparts = gEdgeCounterparts;
      vertexCounterparts = gVertexCounterparts;
      GEntity::setMeshMaster(master, tfo);
  }
  
  inline double myAngle(const SVector3 &a, const SVector3 &b, const SVector3 &d)
  {
      double const cosTheta = dot(a, b);
      double const sinTheta = dot(crossprod(a, b), d);
      return std::atan2(sinTheta, cosTheta);
  }
  
  struct myPlane {
      SPoint3 p;
      SVector3 n;
      double a;
      // nx x + ny y + nz z + a = 0
      myPlane(const SPoint3 &_p, const SVector3 &_n) : p(_p), n(_n)
      {
          n.normalize();
          a = -(n.x() * p.x() + n.y() * p.y() + n.z() * p.z());
      }
      double eval(double x, double y, double z)
      {
          return n.x() * x + n.y() * y + n.z() * z + a;
      }
  };
  
  struct myLine {
      SPoint3 p;
      SVector3 t;
      myLine() : p(0, 0, 0), t(0, 0, 1) {}
      myLine(myPlane &p1, myPlane &p2)
      {
          t = crossprod(p1.n, p2.n);
          if(t.norm() == 0.0) { Msg::Error("parallel planes do not intersect"); }
          else
              t.normalize();
          // find a point, assume z = 0
          double a[2][2] = {{p1.n.x(), p1.n.y()}, {p2.n.x(), p2.n.y()}};
          double b[2] = {-p1.a, -p2.a}, x[2];
          if(!sys2x2(a, b, x)) {
              // assume x = 0
              double az[2][2] = {{p1.n.y(), p1.n.z()}, {p2.n.y(), p2.n.z()}};
              double bz[2] = {-p1.a, -p2.a};
              if(!sys2x2(az, bz, x)) {
                  // assume y = 0
                  double ay[2][2] = {{p1.n.x(), p1.n.z()}, {p2.n.x(), p2.n.z()}};
                  double by[2] = {-p1.a, -p2.a};
                  if(!sys2x2(ay, by, x)) {
                      Msg::Error("parallel planes do not intersect");
                  }
                  else {
                      p = SPoint3(x[0], 0., x[1]);
                  }
              }
              else {
                  p = SPoint3(0., x[0], x[1]);
              }
          }
          else {
              p = SPoint3(x[0], x[1], 0.);
          }
      }
      SPoint3 orthogonalProjection(SPoint3 &a)
      {
          // (x - a) * t = 0 -->
          // x = p + u t --> (p + ut - a) * t = 0 --> u = (a-p) * t
          const double u = dot(a - p, t);
          return SPoint3(p.x() + t.x() * u, p.y() + t.y() * u, p.z() + t.z() * u);
      }
  };
  
  void GFace::setMeshMaster(GFace *master, const std::map<int, int> &edgeCopies)
  {
      std::map<GVertex *, GVertex *> vs2vt;
  
      for(std::vector<GEdge *>::iterator it = l_edges.begin(); it != l_edges.end();
          ++it) {
          // slave edge
          GEdge *le = *it;
  
          int sign = 1;
          std::map<int, int>::const_iterator adnksd = edgeCopies.find(le->tag());
          int source_e;
          if(adnksd != edgeCopies.end())
              source_e = adnksd->second;
          else {
              sign = -1;
              adnksd = edgeCopies.find(-(*it)->tag());
              if(adnksd != edgeCopies.end())
                  source_e = adnksd->second;
              else {
                  Msg::Error("Could not find edge counterpart %d in slave surface %d",
                             (*it)->tag(), master->tag());
                  return;
              }
          }
  
          // master edge
          GEdge *me = master->model()->getEdgeByTag(abs(source_e));
  
          if(le->getBeginVertex() && le->getEndVertex() &&
                  me->getBeginVertex() && me->getEndVertex()){
              if(source_e * sign > 0) {
                  vs2vt[me->getBeginVertex()] = le->getBeginVertex();
                  vs2vt[me->getEndVertex()] = le->getEndVertex();
              }
              else {
                  vs2vt[me->getBeginVertex()] = le->getEndVertex();
                  vs2vt[me->getEndVertex()] = le->getBeginVertex();
              }
          }
      }
  
      // --- find out the transformation
  
      bool translation = true;
      SVector3 DX;
  
      int count = 0;
      for(std::map<GVertex *, GVertex *>::iterator it = vs2vt.begin();
          it != vs2vt.end(); ++it) {
          GVertex *vs = it->first;
          GVertex *vt = it->second;
          if(count == 0)
              DX = SVector3(vt->x() - vs->x(), vt->y() - vs->y(), vt->z() - vs->z());
          else {
              SVector3 DX2(vt->x() - vs->x(), vt->y() - vs->y(), vt->z() - vs->z());
              SVector3 DDX(DX2 - DX);
              if(DDX.norm() > DX.norm() * 1.e-5) translation = false;
          }
          count++;
      }
  
      std::vector<double> tfo(16);
  
      if(translation) {
          Msg::Info(
                      "Periodic mesh translation found between %d and %d: dx = (%g,%g,%g)",
                      tag(), master->tag(), DX.x(), DX.y(), DX.z());
  
          for(size_t i = 0; i < 16; i++) tfo[i] = 0;
          for(size_t i = 0; i < 4; i++) tfo[i * 4 + i] = 1;
          tfo[3] = DX.x();
          tfo[7] = DX.y();
          tfo[11] = DX.z();
      }
  
      else {
          bool rotation = false;
          myLine LINE;
          double ANGLE = 0;
  
          count = 0;
          rotation = true;
          std::vector<SPoint3> mps, mpt;
          for(std::map<GVertex *, GVertex *>::iterator it = vs2vt.begin();
              it != vs2vt.end(); ++it) {
              GVertex *vs = it->first;
              GVertex *vt = it->second;
              mps.push_back(SPoint3(vs->x(), vs->y(), vs->z()));
              mpt.push_back(SPoint3(vt->x(), vt->y(), vt->z()));
          }
          mean_plane mean_source, mean_target;
          computeMeanPlaneSimple(mps, mean_source);
          computeMeanPlaneSimple(mpt, mean_target);
  
          myPlane PLANE_SOURCE(SPoint3(mean_source.x, mean_source.y, mean_source.z),
                               SVector3(mean_source.a, mean_source.b, mean_source.c));
          myPlane PLANE_TARGET(SPoint3(mean_target.x, mean_target.y, mean_target.z),
                               SVector3(mean_target.a, mean_target.b, mean_target.c));
  
          LINE = myLine(PLANE_SOURCE, PLANE_TARGET);
  
          // LINE is the axis of rotation
          // let us compute the angle of rotation
          count = 0;
          for(std::map<GVertex *, GVertex *>::iterator it = vs2vt.begin();
              it != vs2vt.end(); ++it) {
              GVertex *vs = it->first;
              GVertex *vt = it->second;
              // project both points on the axis: that should be the same point !
              SPoint3 ps = SPoint3(vs->x(), vs->y(), vs->z());
              SPoint3 pt = SPoint3(vt->x(), vt->y(), vt->z());
              SPoint3 p_ps = LINE.orthogonalProjection(ps);
              SPoint3 p_pt = LINE.orthogonalProjection(pt);
              SVector3 dist1 = ps - pt;
              SVector3 dist2 = p_ps - p_pt;
              if(dist1.norm() > CTX::instance()->geom.tolerance) {
                  if(dist2.norm() > 1.e-8 * dist1.norm()) { rotation = false; }
                  SVector3 t1 = ps - p_ps;
                  SVector3 t2 = pt - p_pt;
                  if(t1.norm() > 1.e-8 * dist1.norm()) {
                      if(count == 0)
                          ANGLE = myAngle(t1, t2, LINE.t);
                      else {
                          double ANGLE2 = myAngle(t1, t2, LINE.t);
                          if(fabs(ANGLE2 - ANGLE) > 1.e-8) { rotation = false; }
                      }
                      count++;
                  }
              }
          }
  
          if(rotation) {
              Msg::Info("Periodic mesh rotation found: axis (%g,%g,%g) point (%g %g "
                        "%g) angle %g",
                        LINE.t.x(), LINE.t.y(), LINE.t.z(), LINE.p.x(), LINE.p.y(),
                        LINE.p.z(), ANGLE * 180 / M_PI);
  
              double ux = LINE.t.x();
              double uy = LINE.t.y();
              double uz = LINE.t.z();
  
              tfo[0 * 4 + 0] = cos(ANGLE) + ux * ux * (1. - cos(ANGLE));
              tfo[0 * 4 + 1] = ux * uy * (1. - cos(ANGLE)) - uz * sin(ANGLE);
              tfo[0 * 4 + 2] = ux * uz * (1. - cos(ANGLE)) + uy * sin(ANGLE);
  
              tfo[1 * 4 + 0] = ux * uy * (1. - cos(ANGLE)) + uz * sin(ANGLE);
              tfo[1 * 4 + 1] = cos(ANGLE) + uy * uy * (1. - cos(ANGLE));
              tfo[1 * 4 + 2] = uy * uz * (1. - cos(ANGLE)) - ux * sin(ANGLE);
  
              tfo[2 * 4 + 0] = ux * uz * (1. - cos(ANGLE)) - uy * sin(ANGLE);
              tfo[2 * 4 + 1] = uy * uz * (1. - cos(ANGLE)) + ux * sin(ANGLE);
              tfo[2 * 4 + 2] = cos(ANGLE) + uz * uz * (1. - cos(ANGLE));
  
              double origin[3] = {LINE.p.x(), LINE.p.y(), LINE.p.z()};
  
              for(int i = 0; i < 3; i++) tfo[i * 4 + 3] = origin[i];
              for(int i = 0; i < 3; i++)
                  for(int j = 0; j < 3; j++) tfo[i * 4 + 3] -= tfo[i * 4 + j] * origin[j];
              for(int i = 0; i < 3; i++) tfo[12 + i] = 0;
              tfo[15] = 1;
          }
          else {
              Msg::Error("Only rotations or translations can currently be computed "
                         "automatically for periodic faces: face %d not meshed",
                         tag());
              return;
          }
      }
  
      // --- now check and encode the transformation
      // --- including for edges and vertices
  
      setMeshMaster(master, tfo);
  }
  
  void GFace::addElement(int type, MElement *e)
  {
      switch(type) {
      case TYPE_TRI: addTriangle(reinterpret_cast<MTriangle *>(e)); break;
      case TYPE_QUA: addQuadrangle(reinterpret_cast<MQuadrangle *>(e)); break;
      case TYPE_POLYG: addPolygon(reinterpret_cast<MPolygon *>(e)); break;
      default: Msg::Error("Trying to add unsupported element in face");
      }
  }
  
  void GFace::removeElement(int type, MElement *e)
  {
      switch(type) {
      case TYPE_TRI: {
          std::vector<MTriangle *>::iterator it = std::find(
                      triangles.begin(), triangles.end(), reinterpret_cast<MTriangle *>(e));
          if(it != triangles.end()) triangles.erase(it);
      } break;
      case TYPE_QUA: {
          std::vector<MQuadrangle *>::iterator it =
                  std::find(quadrangles.begin(), quadrangles.end(),
                            reinterpret_cast<MQuadrangle *>(e));
          if(it != quadrangles.end()) quadrangles.erase(it);
      } break;
      case TYPE_POLYG: {
          std::vector<MPolygon *>::iterator it = std::find(
                      polygons.begin(), polygons.end(), reinterpret_cast<MPolygon *>(e));
          if(it != polygons.end()) polygons.erase(it);
      } break;
      default: Msg::Error("Trying to remove unsupported element in face");
      }
  }
  
  bool GFace::reorder(const int elementType, const std::vector<std::size_t> &ordering)
  {
      if(triangles.size() != 0) {
          if(triangles.front()->getTypeForMSH() == elementType) {
              if(ordering.size() != triangles.size()) return false;
  
              for(std::vector<std::size_t>::const_iterator it = ordering.begin();
                  it != ordering.end(); ++it) {
                  if(*it >= triangles.size()) return false;
              }
  
              std::vector<MTriangle *> newTrianglesOrder(triangles.size());
              for(std::size_t i = 0; i < ordering.size(); i++) {
                  newTrianglesOrder[i] = triangles[ordering[i]];
              }
  #if __cplusplus >= 201103L
              triangles = std::move(newTrianglesOrder);
  #else
              triangles = newTrianglesOrder;
  #endif
  
              return true;
          }
      }
  
      if(quadrangles.size() != 0) {
          if(quadrangles.front()->getTypeForMSH() == elementType) {
              if(ordering.size() != quadrangles.size()) return false;
  
              for(std::vector<std::size_t>::const_iterator it = ordering.begin();
                  it != ordering.end(); ++it) {
                  if(*it >= quadrangles.size()) return false;
              }
  
              std::vector<MQuadrangle *> newQuadranglesOrder(quadrangles.size());
              for(std::size_t i = 0; i < ordering.size(); i++) {
                  newQuadranglesOrder[i] = quadrangles[ordering[i]];
              }
  #if __cplusplus >= 201103L
              quadrangles = std::move(newQuadranglesOrder);
  #else
              quadrangles = newQuadranglesOrder;
  #endif
  
              return true;
          }
      }
  
      if(polygons.size() != 0) {
          if(polygons.front()->getTypeForMSH() == elementType) {
              if(ordering.size() != polygons.size()) return false;
  
              for(std::vector<std::size_t>::const_iterator it = ordering.begin();
                  it != ordering.end(); ++it) {
                  if(*it >= polygons.size()) return false;
              }
  
              std::vector<MPolygon *> newPolygonsOrder(polygons.size());
              for(std::size_t i = 0; i < ordering.size(); i++) {
                  newPolygonsOrder[i] = polygons[ordering[i]];
              }
  #if __cplusplus >= 201103L
              polygons = std::move(newPolygonsOrder);
  #else
              polygons = newPolygonsOrder;
  #endif
  
              return true;
          }
      }
  
      return false;
  }
  
  void GFace::alignElementsWithMaster()
  {
      GEntity *master = getMeshMaster();
  
      if(master != this) {
          std::set<MFace, Less_Face> srcFaces;
  
          for(std::size_t i = 0; i < master->getNumMeshElements(); i++) {
              MElement *face = master->getMeshElement(i);
              std::vector<MVertex *> vtcs;
              vtcs.reserve(face->getNumVertices());
              for(std::size_t j = 0; j < face->getNumVertices(); j++) {
                  vtcs.push_back(face->getVertex(j));
              }
              srcFaces.insert(MFace(vtcs));
          }
  
          for(std::size_t i = 0; i < getNumMeshElements(); i++) {
              MElement *face = getMeshElement(i);
              std::vector<MVertex *> vtcs;
              for(std::size_t j = 0; j < face->getNumVertices(); j++) {
                  MVertex *tv = face->getVertex(j);
  
                  std::map<MVertex *, MVertex *>::iterator cIter =
                          correspondingVertices.find(tv);
                  if(cIter != correspondingVertices.end()) vtcs.push_back(cIter->second);
              }
  
              MFace mf(vtcs);
  
              std::set<MFace, Less_Face>::iterator sIter = srcFaces.find(mf);
  
              if(sIter == srcFaces.end()) continue;
  
              MFace of = *sIter;
  
              int orientation;
              bool swap;
              mf.computeCorrespondence(of, orientation, swap);
  
              switch(face->getNumVertices()) {
              case 3: {
                  MTriangle *tri = dynamic_cast<MTriangle *>(face);
                  if(tri) tri->reorient(orientation, swap);
                  break;
              }
              case 4: {
                  MQuadrangle *qua = dynamic_cast<MQuadrangle *>(face);
                  if(qua) qua->reorient(orientation, swap);
                  break;
              }
              }
          }
      }
  }
  