New high order tools optimizer (thanks to Thomas and Jon) !!

CMakeLists.txt

@@ -37,6 +37,7 @@ option(ENABLE_FL_TREE "Enable FLTK tree browser widget" ${DEFAULT})
 option(ENABLE_FOURIER_MODEL "Enable Fourier geometrical models" OFF)
 option(ENABLE_GMM "Enable GMM linear algebra solvers" ${DEFAULT})
 option(ENABLE_GRAPHICS "Compile-in OpenGL graphics even if there is no GUI" OFF)
+option(ENABLE_OPTHOM "Enable optimization tool for high order meshes" ON)
 option(ENABLE_KBIPACK "Enable Kbipack for homology solver" ${DEFAULT})
 option(ENABLE_MATHEX "Enable MathEx expression parser" ${DEFAULT})
 option(ENABLE_MED "Enable MED mesh and post-processing file formats" ${DEFAULT})
@@ -108,6 +109,9 @@ set(GMSH_API
   contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h 
     contrib/kbipack/gmp_blas.h contrib/kbipack/mpz.h
   contrib/DiscreteIntegration/Integration3D.h
+  contrib/HighOrderMeshOptimizer/OptHOM.h
+  contrib/HighOrderMeshOptimizer/OptHOM_Mesh.h
+  contrib/HighOrderMeshOptimizer/ParamCoord.h
   contrib/MathEx/mathex.h)
 
 execute_process(COMMAND date "+%Y%m%d" OUTPUT_VARIABLE DATE 
@@ -508,6 +512,13 @@ if(ENABLE_DINTEGRATION)
   set_config_option(HAVE_DINTEGRATION "DIntegration")
 endif(ENABLE_DINTEGRATION)
 
+if(ENABLE_OPTHOM)
+  add_subdirectory(contrib/HighOrderMeshOptimizer)
+  include_directories(contrib/HighOrderMeshOptimizer)
+  set_config_option(HAVE_OPTHOM "OptHom")
+endif(ENABLE_OPTHOM)
+
+
 if(ENABLE_KBIPACK)
   find_library(GMP_LIB gmp)
   if(GMP_LIB)

Fltk/highOrderToolsWindow.cpp

@@ -35,6 +35,7 @@ typedef unsigned long intptr_t;
 #include "Options.h"
 #include "Context.h"
 #include "HighOrder.h"
+#include "OptHomRun.h"
 
 #if defined(HAVE_PARSER)
 #include "Parser.h"
@@ -69,9 +70,17 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data)
   bool elastic = (bool)o->butt[3]->value(); 
   double threshold = o->value[1]->value(); 
   bool onlyVisible = (bool)o->butt[1]->value(); 
-  int nLayers = (int) o->value[2]->value(); 
+  int nbLayers = (int) o->value[2]->value(); 
 
   if(elastic)ElasticAnalogy(GModel::current(), threshold, onlyVisible);
+  else  {
+    OptHomParameters p;
+    p.nbLayers = nbLayers; 
+    p.BARRIER = threshold;
+    p.onlyVisible = onlyVisible;
+    p.dim  = GModel::current()->getNumRegions()  ? 3 : 2;
+    HighOrderMeshOptimizer (GModel::current(),p);
+  }
 
   CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
   drawContext::global()->draw();

contrib/HighOrderMeshOptimizer/CMakeLists.txt

+# Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+#
+# See the LICENSE.txt file for license information. Please report all
+# bugs and problems to <gmsh@geuz.org>.
+
+set(SRC
+  OptHomMesh.cpp 
+  OptHOM.cpp 
+  OptHomRun.cpp 
+  ParamCoord.cpp 
+)
+
+file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.hpp)
+append_gmsh_src(contrib/HighOrderMeshOptimizer "${SRC};${HDR}")

contrib/HighOrderMeshOptimizer/OptHOM.cpp

+#include <iostream>
+#include <algorithm>
+#include "OptHomMesh.h"
+#include "OptHOM.h"
+#include "GmshConfig.h"
+
+#ifdef HAVE_BFGS
+
+#include "ap.h"
+#include "alglibinternal.h"
+#include "alglibmisc.h"
+#include "linalg.h"
+#include "optimization.h"
+
+
+
+// Constructor
+OptHOM::OptHOM(GEntity *ge, std::set<MVertex*> & toFix, int method) :
+       mesh(ge, toFix, method)
+{
+};
+
+// Contribution of the element Jacobians to the objective function value and gradients (2D version)
+bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
+{
+
+  for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
+    std::vector<double> sJ(mesh.nBezEl(iEl));                   // Scaled Jacobians
+    mesh.scaledJac(iEl,sJ);
+    std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));  // Gradients of scaled Jacobians
+    mesh.gradScaledJac(iEl,gSJ);
+    for (int l = 0; l < mesh.nBezEl(iEl); l++) {
+      Obj += compute_f(sJ[l]);
+      const double f1 = compute_f1(sJ[l]);
+      for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++)
+        gradObj[mesh.indPCEl(iEl,iPC)] += f1*gSJ[mesh.indGSJ(iEl,l,iPC)];
+    }
+  }
+
+  return true;
+
+}
+
+
+
+// Contribution of the vertex distance to the objective function value and gradients (2D version)
+bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj)
+{
+
+  for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
+    const double Factor = invLengthScaleSq*(mesh.forced(iFV) ? Fact : Fact2);
+    Obj += Factor * mesh.distSq(iFV);
+    std::vector<double> gDSq(mesh.nPCFV(iFV));
+    mesh.gradDistSq(iFV,gDSq);
+    for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) gradObj[mesh.indPCFV(iFV,iPC)] += Factor*gDSq[iPC];
+  }
+
+  return true;
+
+}
+
+
+
+void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj)
+{
+
+  mesh.updateMesh(x.getcontent());
+
+  Obj = 0.;
+  for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
+
+  addJacObjGrad(Obj, gradObj);
+  addDistObjGrad(lambda, lambda2, Obj, gradObj);
+
+}
+
+
+
+void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj, void *HOInst)
+{
+  OptHOM &HO = *static_cast<OptHOM*> (HOInst);
+  HO.evalObjGrad(x, Obj, gradObj);
+  double distMaxBnd, distAvgBnd, minJac, maxJac;
+  HO.getDistances(distMaxBnd, distAvgBnd, minJac, maxJac);
+  if (minJac > HO.barrier) {
+    printf("STOP\n");
+    for (int i = 0; i < gradObj.length(); ++i) {
+      gradObj[i] = 0;
+    }
+  }
+}
+
+
+
+void OptHOM::getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac)
+{
+
+  distMaxBND = 0;
+  distAvgBND = 0;
+  int nbBnd = 0;
+
+  for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
+    if (mesh.forced(iFV)) {
+      double dSq = mesh.distSq(iFV);
+      distMaxBND = std::max(distMaxBND, sqrt(dSq));
+      distAvgBND += sqrt(dSq);
+      nbBnd++;
+    }
+  }
+  if (nbBnd != 0) distAvgBND /= nbBnd;
+
+  minJac = 1000;
+  maxJac = -1000;
+  for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
+    std::vector<double> sJ(mesh.nBezEl(iEl));                   // Scaled Jacobians
+    mesh.scaledJac(iEl,sJ);
+    for (int l = 0; l < mesh.nBezEl(iEl); l++) {
+      minJac = std::min(minJac, sJ[l]);
+      maxJac = std::max(maxJac, sJ[l]);
+    }
+  }
+
+}
+
+
+
+void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj)
+{
+
+  iter++;
+
+  if (iter % progressInterv == 0) {
+    double maxD, avgD, minJ, maxJ;
+    getDistances(maxD, avgD, minJ, maxJ);
+
+    printf("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E  maxJ = %12.5E Max D = %12.5E Avg D = %12.5E\n", iter, Obj, Obj/initObj, minJ, maxJ, maxD, avgD);
+  }
+
+}
+
+
+
+void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst)
+{
+  ((OptHOM*)HOInst)->printProgress(x,Obj);
+}
+
+
+
+void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax)
+{
+
+  static const double EPSG = 0.;
+  static const double EPSF = 1.e-8;
+//  static const double EPSF = 0.;
+  static const double EPSX = 0.;
+//  const double EPSX = x.length()*1.e-4/sqrt(invLengthScaleSq);
+//  std::cout << "DEBUG: EPSX = " << EPSX << ", EPSX/x.length() = " << EPSX/x.length() << std::endl;
+
+//  double iGONorm = 0;
+//  for (int i=0; i<initGradObj.length(); i++) iGONorm += initGradObj[i]*initGradObj[i];
+//  const double EPSG = 1.e-2*sqrt(iGONorm);
+
+  std::cout << "--- Optimization pass with jacBar = " << jacBar << ", lambda = " << lambda << ", lambda2 = " << lambda2 << std::endl;
+
+//  alglib::minlbfgsstate state;
+//  alglib::minlbfgsreport rep;
+  alglib::mincgstate state;
+  alglib::mincgreport rep;
+
+//  minlbfgscreate(3, x, state);
+//  minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
+//  minlbfgssetxrep(state, true);
+  mincgcreate(x, state);
+  mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
+  mincgsetxrep(state, true);
+
+  iter = 0;
+
+//  alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
+  alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
+
+//  minlbfgsresults(state, x, rep);
+  mincgresults(state, x, rep);
+
+  std::cout << "Optimization finalized after " << rep.iterationscount << " iterations (" << rep.nfev << " functions evaluations)";
+  switch(int(rep.terminationtype)) {
+//  case -2: std::cout << ", because rounding errors prevented further improvement"; break;
+//  case -1: std::cout << ", because incorrect parameters were specified"; break;
+//  case 1: std::cout << ", because relative function improvement is no more than EpsF"; break;
+//  case 2: std::cout << ", because relative step is no more than EpsX"; break;
+//  case 4: std::cout << ", because gradient norm is no more than EpsG"; break;
+//  case 5: std::cout << ", because the maximum number of steps was taken"; break;
+//  case 7: std::cout << ", because stopping conditions are too stringent, further improvement is impossible"; break;
+  default: std::cout << " with code " << int(rep.terminationtype); break;
+  }
+  std::cout << "." << std::endl;
+
+}
+
+
+
+int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int pInt, int itMax)
+{
+  barrier = barrier_;
+  progressInterv = pInt;
+//  powM = 4;
+//  powP = 3;
+
+  // Set weights & length scale for non-dimensionalization
+  lambda = weightFixed;
+  lambda2 = weightFree;
+  std::vector<double> dSq(mesh.nVert());
+  mesh.distSqToStraight(dSq);
+  const double maxDSq = *max_element(dSq.begin(),dSq.end());
+  invLengthScaleSq = 1./maxDSq;  // Length scale for non-dimensional distance
+
+  // Set initial guess
+  alglib::real_1d_array x;
+  x.setlength(mesh.nPC());
+  mesh.getUvw(x.getcontent());
+
+  // Calculate initial performance
+  double minJ, maxJ;
+  getDistances(initMaxD, initAvgD, minJ, maxJ);
+
+  const double jacBarStart = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
+  jacBar = jacBarStart;
+  setBarrierTerm(jacBarStart);
+  std::cout << "DEBUG: jacBarStart = " << jacBarStart << std::endl;
+
+  // Calculate initial objective function value and gradient
+  initObj = 0.;
+  alglib::real_1d_array gradObj;
+  gradObj.setlength(mesh.nPC());
+  for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
+  evalObjGrad(x, initObj, gradObj);
+
+  std::cout << "Initial mesh: Obj = " << initObj << ", minJ = " << minJ << ", maxJ = " << maxJ << ", maxD = " << initMaxD << ", avgD = " << initAvgD << std::endl;
+
+  std::cout << "Start optimizing with barrier = " << barrier << std::endl;
+
+  while (minJ < barrier) {
+    OptimPass(x, gradObj, itMax);
+    double dumMaxD, dumAvgD;
+    getDistances(dumMaxD, dumAvgD, minJ, maxJ);
+    jacBar = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
+    setBarrierTerm(jacBar);
+  }
+
+//  for (int i = 0; i<3; i++) {
+//    lambda *= 10;
+//    OptimPass(x, gradObj, itMax);
+//  }
+
+  double finalObj = 0., finalMaxD, finalAvgD;
+  evalObjGrad(x, finalObj, gradObj);
+  getDistances(finalMaxD, finalAvgD, minJ, maxJ);
+  std::cout << "Final mesh: Obj = " << finalObj << ", maxD = " << finalMaxD << ", avgD = " << finalAvgD << ", minJ = " << minJ << ", maxJ = " << maxJ << std::endl;
+
+  return 0;
+
+}
+
+
+
+#endif

contrib/HighOrderMeshOptimizer/OptHOM.h

+#ifndef _OPTHOM_H_
+#define _OPTHOM_H_
+
+#include <string>
+#include <math.h>
+
+
+
+#ifdef HAVE_BFGS
+
+#include "ap.h"
+
+
+
+class Mesh;
+
+class OptHOM
+{
+
+public:
+
+  Mesh mesh;
+
+  OptHOM(GEntity *gf, std::set<MVertex*> & toFix, int method);
+  void getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac);
+  int optimize(double lambda, double lambda2, double barrier, int pInt, int itMax);  // optimize one list of elements
+  void updateMesh(const alglib::real_1d_array &x);
+  void evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj);
+  void printProgress(const alglib::real_1d_array &x, double Obj);
+
+  double barrier;
+
+private:
+
+//  double lambda, lambda2, powM, powP, invLengthScaleSq;
+  double lambda, lambda2, jacBar, bTerm, invLengthScaleSq;
+  int iter, progressInterv;            // Current iteration, interval of iterations for reporting
+  double initObj, initMaxD, initAvgD;  // Values for reporting
+
+  inline double setBarrierTerm(double jacBarrier) { bTerm = jacBarrier/(1.-jacBarrier); }
+  inline double compute_f(double v);
+  inline double compute_f1(double v);
+  bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
+  bool addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj);
+  void OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax);
+
+};
+
+
+
+inline double OptHOM::compute_f(double v)
+{
+  if (v > jacBar) {
+    const double l = log((1 + bTerm) * v - bTerm);
+    const double m = (v - 1);
+    return l * l + m * m;
+  }
+  else return 1.e300;
+//  if (v < 1.) return pow(1.-v,powM);
+//  if (v < 1.) return exp((long double)pow(1.-v,3));
+//  else return pow(v-1.,powP);
+}
+
+
+
+inline double OptHOM::compute_f1(double v)
+{
+  if (v > jacBar) {
+    const double veps = (1 + bTerm) * v - bTerm;
+    const double m = 2 * (v - 1);
+    return m + 2 * log(veps) / veps * (1 + bTerm);
+  }
+  else return -1.e300;
+//  if (v < 1.) return -powM*pow(1.-v,powM-1.);
+//  if (v < 1.) return -3.*pow(1.-v,2)*exp((long double)pow(1.-v,3));
+//  else return powP*pow(v-1.,powP-1.);
+}
+
+
+
+#endif
+
+
+
+#endif

contrib/HighOrderMeshOptimizer/OptHomMesh.cpp

+#include "GRegion.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "ParamCoord.h"
+#include "OptHomMesh.h"
+#include "JacobianBasis.h"
+
+std::map<int, std::vector<double> > Mesh::_jacBez;
+
+static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
+{
+  int nbNodes = lagrange->points.size1();
+  
+  // bezier points are defined in the [0,1] x [0,1] quad
+  fullMatrix<double> bezierPoints = bezier->points;
+  if (lagrange->parentType == TYPE_QUA) {
+    for (int i = 0; i < bezierPoints.size1(); ++i) {
+      bezierPoints(i, 0) = -1 + 2 * bezierPoints(i, 0);
+      bezierPoints(i, 1) = -1 + 2 * bezierPoints(i, 1);
+    }
+  }
+
+  fullMatrix<double> allDPsi;
+  lagrange->df(bezierPoints, allDPsi);
+  int size = bezier->points.size1();
+  std::vector<double> JB;
+  for (int d = 0; d < lagrange->dimension; ++d) {
+    size *= nbNodes;
+  }
+  JB.resize(size, 0.);
+  for (int k = 0; k < bezier->points.size1(); ++k) {
+    fullMatrix<double> dPsi(allDPsi, k * 3, 3);
+    if (lagrange->dimension == 2) {
+      for (int i = 0; i < nbNodes; i++) {
+        for (int j = 0; j < nbNodes; j++) {
+          double Jij = dPsi(i, 0) * dPsi(j, 1) - dPsi(i, 1) * dPsi(j,0);
+          for (int l = 0; l < bezier->points.size1(); l++) {
+            JB[(l * nbNodes + i) * nbNodes + j] += bezier->matrixLag2Bez(l, k) * Jij;
+          }
+        }
+      }
+    }
+    if (lagrange->dimension == 3) {
+      for (int i = 0; i < nbNodes; i++) {
+        for (int j = 0; j < nbNodes; j++) {
+          for (int m = 0; m < nbNodes; m++) {
+            double Jijm = 
+                (dPsi(j, 1) * dPsi(m, 2) - dPsi(j, 2) * dPsi(m, 1)) * dPsi(i, 0)
+              + (dPsi(j, 2) * dPsi(m, 0) - dPsi(j, 0) * dPsi(m, 2)) * dPsi(i, 1)
+              + (dPsi(j, 0) * dPsi(m, 1) - dPsi(j, 1) * dPsi(m, 0)) * dPsi(i, 2);
+            for (int l = 0; l < bezier->points.size1(); l++) {
+              JB[(l * nbNodes + (i * nbNodes + j)) * nbNodes + m] += bezier->matrixLag2Bez(l, k) * Jijm;
+            }
+          }
+        }
+      }
+    }
+  }
+  return JB;
+}
+
+Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
+    _ge(ge)
+{
+  _dim = _ge->dim();
+
+  if (method & METHOD_PHYSCOORD) {
+    if (_dim == 2) _pc = new ParamCoordPhys2D;
+    else _pc = new ParamCoordPhys3D;
+    std::cout << "METHOD: Using physical coordinates" << std::endl;
+  }
+  else if (method & METHOD_SURFCOORD)  {
+    if (_dim == 2) {
+      _pc = new ParamCoordSurf(_ge);
+      std::cout << "METHOD: Using surface parametric coordinates" << std::endl;
+    }
+    else std::cout << "ERROR: Surface parametric coordinates only for 2D optimization" << std::endl;
+  }
+  else {
+    _pc = new ParamCoordParent;
+    std::cout << "METHOD: Using parent parametric coordinates" << std::endl;
+  }
+
+  if (method & METHOD_RELAXBND) std::cout << "METHOD: Relaxing boundary vertices" << std::endl;
+  else if (method & METHOD_FIXBND) std::cout << "METHOD: Fixing all boundary vertices" << std::endl;
+  else std::cout << "METHOD: Fixing vertices on geometric points and \"toFix\" boundary" << std::endl;
+
+  // Initialize elements, vertices, free vertices and element->vertices connectivity
+  _nPC = 0;
+  _el.resize(_ge->getNumMeshElements());
+  _el2FV.resize(_ge->getNumMeshElements());
+  _el2V.resize(_ge->getNumMeshElements());
+  _nBezEl.resize(_ge->getNumMeshElements());
+  _nNodEl.resize(_ge->getNumMeshElements());
+  _indPCEl.resize(_ge->getNumMeshElements());
+  for (int iEl = 0; iEl < nEl(); iEl++) {
+    MElement *el = _ge->getMeshElement(iEl);
+    _el[iEl] = el;
+    const polynomialBasis *lagrange = el->getFunctionSpace();
+    const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
+    if (_jacBez.find(lagrange->type) == _jacBez.end()) {
+      _jacBez[lagrange->type] = _computeJB(lagrange, bezier);
+    }
+    _nBezEl[iEl] = bezier->points.size1();
+    _nNodEl[iEl] = lagrange->points.size1();
+    for (int iVEl = 0; iVEl < lagrange->points.size1(); iVEl++) {
+      MVertex *vert = el->getVertex(iVEl);
+      int iV = addVert(vert);
+      _el2V[iEl].push_back(iV);
+      const int nPCV = _pc->nCoord(vert);
+      bool isFV = false;
+      if (method & METHOD_RELAXBND) isFV = true;
+      else if (method & METHOD_FIXBND) isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
+      else isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
+      if (isFV) {
+        int iFV = addFreeVert(vert,iV,nPCV,toFix);
+        _el2FV[iEl].push_back(iFV);
+        for (int i=_startPCFV[iFV]; i<_startPCFV[iFV]+nPCV; i++) _indPCEl[iEl].push_back(i);
+      }
+      else _el2FV[iEl].push_back(-1);
+    }
+  }
+
+  // Initial coordinates
+  _ixyz.resize(nVert());
+  for (int iV = 0; iV < nVert(); iV++) _ixyz[iV] = _vert[iV]->point();
+  _iuvw.resize(nFV());
+  for (int iFV = 0; iFV < nFV(); iFV++) _iuvw[iFV] = _pc->getUvw(_freeVert[iFV]);
+
+  // Set current coordinates
+  _xyz = _ixyz;
+  _uvw = _iuvw;
+
+  // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial Jacobians of 3D elements
+  if (_dim == 2) {
+    _normEl.resize(nEl());
+    for (int iEl = 0; iEl < nEl(); iEl++) _normEl[iEl] = getNormalEl(iEl);
+  }
+  else {
+    _invStraightJac.resize(nEl(),1.);
+    double dumJac[3][3];
+    for (int iEl = 0; iEl < nEl(); iEl++) _invStraightJac[iEl] = 1. / _el[iEl]->getPrimaryJacobian(0.,0.,0.,dumJac);
+  }
+  if ((_dim == 2) && (method & METHOD_PROJJAC)) {
+    projJac = true;
+    std::cout << "METHOD: Using projected Jacobians" << std::endl;
+  }
+  else {
+    projJac = false;
+    std::cout << "METHOD: Using usual Jacobians" << std::endl;
+  }
+}
+
+SVector3 Mesh::getNormalEl(int iEl)
+{
+
+  switch (_el[iEl]->getType()) {
+    case TYPE_TRI: {
+      const int iV0 = _el2V[iEl][0], iV1 = _el2V[iEl][1], iV2 = _el2V[iEl][2];
+      SVector3 v10 = _xyz[iV1]-_xyz[iV0], v20 = _xyz[iV2]-_xyz[iV0];
+      SVector3 n = crossprod(v10, v20);
+      double xxx = n.norm();
+      n *= 1./(xxx*xxx);
+      return n;
+      break;
+    }
+    case TYPE_QUA: {
+      const int iV0 = _el2V[iEl][0], iV1 = _el2V[iEl][1], iV3 = _el2V[iEl][3];
+      SVector3 v10 = _xyz[iV1]-_xyz[iV0], v30 = _xyz[iV3]-_xyz[iV0];
+      SVector3 n = crossprod(v10, v30);
+      double xxx = n.norm();
+      n *= 4./(xxx*xxx);
+      return n;
+      break;
+    }
+    case TYPE_TET: {
+      return SVector3(0.);
+      break;
+    }
+    default:
+      std::cout << "ERROR: getNormalEl: Unknown element type" << std::endl;
+      break;
+  }
+
+  return SVector3(0.);  // Just to avoid compilation warnings...
+
+}
+
+
+
+int Mesh::addVert(MVertex* vert)
+{
+
+  std::vector<MVertex*>::iterator itVert = find(_vert.begin(),_vert.end(),vert);
+  if (itVert == _vert.end()) {
+    _vert.push_back(vert);
+    return _vert.size()-1;
+  }
+  else return std::distance(_vert.begin(),itVert);
+
+}
+
+
+
+int Mesh::addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix)
+{
+
+  std::vector<MVertex*>::iterator itVert = find(_freeVert.begin(),_freeVert.end(),vert);
+  if (itVert == _freeVert.end()) {
+    const int iStart = (_startPCFV.size() == 0)? 0 : _startPCFV.back()+_nPCFV.back();
+    const bool forcedV = (vert->onWhat()->dim() < 2) || (toFix.find(vert) != toFix.end());
+    _freeVert.push_back(vert);
+    _fv2V.push_back(iV);
+    _startPCFV.push_back(iStart);
+    _nPCFV.push_back(nPCV);
+    _nPC += nPCV;
+    _forced.push_back(forcedV);
+    return _freeVert.size()-1;
+  }
+  else return std::distance(_freeVert.begin(),itVert);
+
+}
+
+
+
+void Mesh::getUvw(double *it)
+{
+
+//  std::vector<double>::iterator it = uvw.begin();
+  for (int iFV = 0; iFV < nFV(); iFV++) {
+    SPoint3 &uvwV = _uvw[iFV];
+    *it = uvwV[0]; it++;
+    if (_nPCFV[iFV] >= 2) { *it = uvwV[1]; it++; }
+    if (_nPCFV[iFV] == 3) { *it = uvwV[2]; it++; }
+  }
+
+}
+
+
+
+void Mesh::updateMesh(const double *it)
+{
+
+//  std::vector<double>::const_iterator it = uvw.begin();
+  for (int iFV = 0; iFV < nFV(); iFV++) {
+    int iV = _fv2V[iFV];
+    SPoint3 &uvwV = _uvw[iFV];
+    uvwV[0] = *it; it++;
+    if (_nPCFV[iFV] >= 2) { uvwV[1] = *it; it++; }
+    if (_nPCFV[iFV] == 3) { uvwV[2] = *it; it++; }
+    _xyz[iV] = _pc->uvw2Xyz(_freeVert[iFV],uvwV);
+  }
+
+}
+
+
+
+void Mesh::distSqToStraight(std::vector<double> &dSq)
+{
+  std::vector<SPoint3> sxyz(nVert());
+  for (int iEl = 0; iEl < nEl(); iEl++) {
+    MElement *el = _el[iEl];
+    const polynomialBasis *lagrange = el->getFunctionSpace();
+    const polynomialBasis *lagrange1 = el->getFunctionSpace(1);
+    int nV = lagrange->points.size1();
+    int nV1 = lagrange1->points.size1();
+    for (int i = 0; i < nV1; ++i) {
+      sxyz[_el2V[iEl][i]] = _vert[_el2V[iEl][i]]->point();
+    }
+    for (int i = nV1; i < nV; ++i) {
+      double f[256];
+      lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1), lagrange->points(i, 2), f);
+      for (int j = 0; j < nV1; ++j)
+        sxyz[_el2V[iEl][i]] += sxyz[_el2V[iEl][j]] * f[j];
+    }
+  }
+
+  for (int iV = 0; iV < nVert(); iV++) {
+    SPoint3 d = _xyz[iV]-sxyz[iV];
+    dSq[iV] = d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
+  }
+}
+
+
+
+void Mesh::updateGEntityPositions()
+{
+
+  for (int iV = 0; iV < nVert(); iV++) _vert[iV]->setXYZ(_xyz[iV].x(),_xyz[iV].y(),_xyz[iV].z());
+
+}
+
+
+
+void Mesh::scaledJac(int iEl, std::vector<double> &sJ)
+{
+  const std::vector<double> &jacBez = _jacBez[_el[iEl]->getTypeForMSH()];
+  if (_dim == 2) {
+    SVector3 &n = _normEl[iEl];
+    if (projJac) {
+      for (int l = 0; l < _nBezEl[iEl]; l++) {
+        sJ[l] = 0.;
+        for (int i = 0; i < _nNodEl[iEl]; i++) {
+          int &iVi = _el2V[iEl][i];
+          for (int j = 0; j < _nNodEl[iEl]; j++) {
+            int &iVj = _el2V[iEl][j];
+            sJ[l] += jacBez[indJB2D(iEl,l,i,j)]
+                  * (_xyz[iVi].x() * _xyz[iVj].y() * n.z() - _xyz[iVi].x() * _xyz[iVj].z() * n.y()
+                     + _xyz[iVi].y() * _xyz[iVj].z() * n.x());
+          }
+        }
+      }
+    }
+    else
+      for (int l = 0; l < _nBezEl[iEl]; l++) {
+        sJ[l] = 0.;
+        for (int i = 0; i < _nNodEl[iEl]; i++) {
+          int &iVi = _el2V[iEl][i];
+          for (int j = 0; j < _nNodEl[iEl]; j++) {
+            int &iVj = _el2V[iEl][j];
+            sJ[l] += jacBez[indJB2D(iEl,l,i,j)] * _xyz[iVi].x() * _xyz[iVj].y();
+          }
+        }
+        sJ[l] *= n.z();
+      }
+ }
+  else {
+    for (int l = 0; l < _nBezEl[iEl]; l++) {
+      sJ[l] = 0.;
+      for (int i = 0; i < _nNodEl[iEl]; i++) {
+        int &iVi = _el2V[iEl][i];
+        for (int j = 0; j < _nNodEl[iEl]; j++) {
+          int &iVj = _el2V[iEl][j];
+          for (int m = 0; m < _nNodEl[iEl]; m++) {
+            int &iVm = _el2V[iEl][m];
+            sJ[l] += jacBez[indJB3D(iEl,l,i,j,m)] * _xyz[iVi].x() * _xyz[iVj].y() * _xyz[iVm].z();
+          }
+        }
+      }
+      sJ[l] *= _invStraightJac[iEl];
+    }
+  }
+
+}
+
+
+
+void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
+{
+  const std::vector<double> &jacBez = _jacBez[_el[iEl]->getTypeForMSH()];
+  if (_dim == 2) {
+    int iPC = 0;
+    SVector3 n = _normEl[iEl];
+    if (projJac) {
+      for (int i = 0; i < _nNodEl[iEl]; i++) {
+        int &iFVi = _el2FV[iEl][i];
+        if (iFVi >= 0) {
+          std::vector<SPoint3> gXyzV(_nBezEl[iEl],SPoint3(0.,0.,0.));
+          std::vector<SPoint3> gUvwV(_nBezEl[iEl]);
+          for (int m = 0; m < _nNodEl[iEl]; m++) {
+            int &iVm = _el2V[iEl][m];
+            const double vpx = _xyz[iVm].y() * n.z() - _xyz[iVm].z() * n.y();
+            const double vpy = -_xyz[iVm].x() * n.z() + _xyz[iVm].z() * n.x();
+            const double vpz = _xyz[iVm].x() * n.y() - _xyz[iVm].y() * n.x();
+            for (int l = 0; l < _nBezEl[iEl]; l++) {
+              gXyzV[l][0] += jacBez[indJB2D(iEl,l,i,m)] * vpx;
+              gXyzV[l][1] += jacBez[indJB2D(iEl,l,i,m)] * vpy;
+              gXyzV[l][2] += jacBez[indJB2D(iEl,l,i,m)] * vpz;
+            }
+          }
+          _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
+          for (int l = 0; l < _nBezEl[iEl]; l++) {
+            gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
+            if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+          }
+          iPC += _nPCFV[iFVi];
+        }
+      }
+    }
+    else
+      for (int i = 0; i < _nNodEl[iEl]; i++) {
+        int &iFVi = _el2FV[iEl][i];
+        if (iFVi >= 0) {
+          std::vector<SPoint3> gXyzV(_nBezEl[iEl],SPoint3(0.,0.,0.));
+          std::vector<SPoint3> gUvwV(_nBezEl[iEl]);
+          for (int m = 0; m < _nNodEl[iEl]; m++) {
+            int &iVm = _el2V[iEl][m];
+            for (int l = 0; l < _nBezEl[iEl]; l++) {
+              gXyzV[l][0] += jacBez[indJB2D(iEl,l,i,m)] * _xyz[iVm].y() * n.z();
+              gXyzV[l][1] += jacBez[indJB2D(iEl,l,m,i)] * _xyz[iVm].x() * n.z();
+            }
+          }
+          _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
+          for (int l = 0; l < _nBezEl[iEl]; l++) {
+            gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
+            if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+          }
+          iPC += _nPCFV[iFVi];
+        }
+      }
+  }
+  else {
+    int iPC = 0;
+    for (int i = 0; i < _nNodEl[iEl]; i++) {
+      int &iFVi = _el2FV[iEl][i];
+      if (iFVi >= 0) {
+        std::vector<SPoint3> gXyzV(_nBezEl[iEl],SPoint3(0.,0.,0.));
+        std::vector<SPoint3> gUvwV(_nBezEl[iEl]);
+        for (int a = 0; a < _nNodEl[iEl]; a++) {
+          int &iVa = _el2V[iEl][a];
+          for (int b = 0; b < _nNodEl[iEl]; b++) {
+            int &iVb = _el2V[iEl][b];
+            for (int l = 0; l < _nBezEl[iEl]; l++) {
+              gXyzV[l][0] += jacBez[indJB3D(iEl,l,i,a,b)] * _xyz[iVa].y() * _xyz[iVb].z() * _invStraightJac[iEl];
+              gXyzV[l][1] += jacBez[indJB3D(iEl,l,a,i,b)] * _xyz[iVa].x() * _xyz[iVb].z() * _invStraightJac[iEl];
+              gXyzV[l][2] += jacBez[indJB3D(iEl,l,a,b,i)] * _xyz[iVa].x() * _xyz[iVb].y() * _invStraightJac[iEl];
+            }
+          }
+        }
+        _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
+        for (int l = 0; l < _nBezEl[iEl]; l++) {
+          gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
+          if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][1];
+          if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][2];
+        }
+        iPC += _nPCFV[iFVi];
+      }
+    }
+  }
+
+}
+
+void Mesh::writeMSH(const char *filename)
+{
+  FILE *f = fopen(filename, "w");
+
+  fprintf(f, "$MeshFormat\n");
+  fprintf(f, "2.2 0 8\n");
+  fprintf(f, "$EndMeshFormat\n");
+
+  fprintf(f, "$Nodes\n");
+  fprintf(f, "%d\n", nVert());
+  for (int i = 0; i < nVert(); i++)
+    fprintf(f, "%d %22.15E %22.15E %22.15E\n", i + 1, _xyz[i].x(), _xyz[i].y(), _xyz[i].z());
+  fprintf(f, "$EndNodes\n");
+
+  fprintf(f, "$Elements\n");
+  fprintf(f, "%d\n", nEl());
+  for (int iEl = 0; iEl < nEl(); iEl++) {
+//    MElement *MEl = _el[iEl];
+    fprintf(f, "%d %d 2 0 %d", iEl+1, _el[iEl]->getTypeForMSH(), _ge->tag());
+    for (int iVEl = 0; iVEl < _el2V[iEl].size(); iVEl++) fprintf(f, " %d", _el2V[iEl][iVEl] + 1);
+    fprintf(f, "\n");
+  }
+  fprintf(f, "$EndElements\n");
+
+  fclose(f);
+
+}

contrib/HighOrderMeshOptimizer/OptHomMesh.h

+#ifndef _OPT_HOM_MESH_H_
+#define _OPT_HOM_MESH_H_
+
+#include <vector>
+#include <map>
+#include <bitset>
+#include "GFace.h"
+#include "MVertex.h"
+#include "ParamCoord.h"
+
+class Mesh
+{
+
+public:
+
+  Mesh(GEntity *ge, std::set<MVertex*> & toFix, int method);
+
+  inline const int &nPC() { return _nPC; }
+  inline int nVert() { return _vert.size(); }
+  inline int nFV() { return _freeVert.size(); }
+  inline const bool forced(int iV) { return _forced[iV]; }
+  inline int nEl() { return _el.size(); }
+  inline const int &nPCFV(int iFV) { return _nPCFV[iFV]; }
+  inline int indPCFV(int iFV, int iPC) { return _startPCFV[iFV]+iPC; }
+  inline int nPCEl(int iEl) { return _indPCEl[iEl].size(); }
+  inline const int &indPCEl(int iEl, int iPC) { return _indPCEl[iEl][iPC]; }
+  inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
+
+  void scaledJac(int iEl, std::vector<double> &sJ);
+  void gradScaledJac(int iEl, std::vector<double> &gSJ);
+  inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; }
+
+  inline double distSq(int iFV);
+  inline void gradDistSq(int iV, std::vector<double> &gDSq);
+
+  void getUvw(double *it);
+  void updateMesh(const double *it);
+  void distSqToStraight(std::vector<double> &dSq);
+
+  void updateGEntityPositions();
+  void writeMSH(const char *filename);
+
+  enum { METHOD_RELAXBND = 1 << 0, METHOD_FIXBND = 1 << 1, METHOD_PHYSCOORD = 1 << 2,
+    METHOD_SURFCOORD = 1 << 3, METHOD_PROJJAC = 1 << 4 };
+
+//  inline double xyzDBG(int iV, int iCoord) { return _xyz[iV][iCoord]; }
+//  inline double ixyzDBG(int iV, int iCoord) { return _ixyz[iV][iCoord]; }
+//  inline double sxyzDBG(int iV, int iCoord) { return _sxyz[iV][iCoord]; }
+//  inline double uvwDBG(int iFV, int iCoord) { return _uvw[iFV][iCoord]; }
+//  inline int fv2VDBG(int iFV) { return _fv2V[iFV]; }
+//  inline bool forcedDBG(int iV) { return _forced[iV]; }
+
+private:
+
+  GEntity *_ge;
+  int _dim;
+  int _nPC;                                     // Total nb. of parametric coordinates
+
+  std::vector<MElement*> _el;                   // List of elements
+  std::vector<SVector3> _normEl;                // Normals to 2D elements
+  std::vector<double> _invStraightJac;          // Initial Jacobians for 3D elements
+  std::vector<MVertex*> _vert, _freeVert;       // List of vert., free vert.
+  std::vector<int> _fv2V;                       // Index of free vert. -> index of vert.
+  std::vector<bool> _forced;                    // Is vertex forced?
+  std::vector<SPoint3> _xyz, _ixyz;             // Physical coord. of ALL vertices (current, straight, init.)
+  std::vector<SPoint3> _uvw, _iuvw;             // Parametric coord. of FREE vertices (current, straight, init.)
+  std::vector<int> _startPCFV;                  // Start index of parametric coordinates for a free vertex
+  std::vector<int> _nPCFV;                      // Number of parametric coordinates for a free vertex
+  std::vector<std::vector<int> > _el2FV, _el2V; // Free vertices, vertices in element
+  std::vector<int> _nBezEl, _nNodEl;            // Number of Bezier poly. and nodes for an el.
+  std::vector<std::vector<int> > _indPCEl;      // Index of parametric coord. for an el.
+
+  ParamCoord *_pc;
+  bool projJac;                                 // Using "projected" Jacobians or not
+
+  static std::map<int, std::vector<double> > _jacBez;
+
+  int addVert(MVertex* vert);
+  int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
+  SVector3 getNormalEl(int iEl);
+  inline int indJB2D(int iEl, int l, int i, int j) { return (l*_nNodEl[iEl]+i)*_nNodEl[iEl]+j; }
+  inline int indJB3D(int iEl, int l, int i, int j, int m) { return ((l*_nNodEl[iEl]+i)*_nNodEl[iEl]+j)*_nNodEl[iEl]+m; }
+
+};
+
+
+
+double Mesh::distSq(int iFV)
+{
+
+  SPoint3 d = _xyz[_fv2V[iFV]]-_ixyz[_fv2V[iFV]];
+  return d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
+
+}
+
+
+
+void Mesh::gradDistSq(int iFV, std::vector<double> &gDSq)
+{
+
+  SPoint3 gXyz = (_xyz[_fv2V[iFV]]-_ixyz[_fv2V[iFV]])*2.;
+  SPoint3 gUvw;
+  _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],gXyz,gUvw);
+
+  gDSq[0] = gUvw[0];
+  if (_nPCFV[iFV] >= 2) gDSq[1] = gUvw[1];
+  if (_nPCFV[iFV] == 3) gDSq[2] = gUvw[2];
+
+}
+
+
+
+#endif

contrib/HighOrderMeshOptimizer/OptHomRun.cpp

+//#include <fenv.h>
+#include <stdio.h>
+#include <sstream>
+#include <string.h>
+#include "OptHomRun.h"
+#include "GModel.h"
+#include "Gmsh.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "highOrderTools.h"
+#include "OptHomMesh.h"
+#include "OptHOM.h"
+
+std::set<MVertex*> filter2D(GFace *gf, int nbLayers, double _qual, bool onlytheworst = false)
+{
+  std::set<MVertex*> touched;
+  std::set<MVertex*> not_touched;
+  std::set<MTriangle *> ts;
+  std::set<MQuadrangle *> qs;
+
+  if (onlytheworst) {
+    double minQ = 1.e22;
+    MQuadrangle *worst;
+    for (int i = 0; i < gf->quadrangles.size(); i++) {
+      MQuadrangle *q = gf->quadrangles[i];
+      double Q = q->distoShapeMeasure();
+      if (Q < minQ) {
+        worst = q;
+        minQ = Q;
+      }
+    }
+    qs.insert(worst);
+    for (int j = 0; j < worst->getNumVertices(); j++)
+      touched.insert(worst->getVertex(j));
+  }
+
+  else {
+    for (int i = 0; i < gf->triangles.size(); i++) {
+      MTriangle *t = gf->triangles[i];
+      double Q = t->distoShapeMeasure();
+      if (Q < _qual || Q > 1.01) {
+        ts.insert(t);
+        for (int j = 0; j < t->getNumVertices(); j++)
+          touched.insert(t->getVertex(j));
+      }
+    }
+    for (int i = 0; i < gf->quadrangles.size(); i++) {
+      MQuadrangle *q = gf->quadrangles[i];
+      double Q = q->distoShapeMeasure();
+      if (Q < _qual || Q > 1.0) {
+        qs.insert(q);
+        for (int j = 0; j < q->getNumVertices(); j++)
+          touched.insert(q->getVertex(j));
+      }
+    }
+  }
+  printf("FILTER2D : %d bad triangles found among %6d\n", ts.size(), gf->triangles.size());
+  printf("FILTER2D : %d bad quads     found among %6d\n", qs.size(), gf->quadrangles.size());
+
+  // add layers of elements around bad elements
+  for (int layer = 0; layer < nbLayers; layer++) {
+    for (int i = 0; i < gf->triangles.size(); i++) {
+      MTriangle *t = gf->triangles[i];
+      bool add_ = false;
+      for (int j = 0; j < t->getNumVertices(); j++) {
+        if (touched.find(t->getVertex(j)) != touched.end()) {
+          add_ = true;
+        }
+      }
+      if (add_) ts.insert(t);
+    }
+    for (int i = 0; i < gf->quadrangles.size(); i++) {
+      bool add_ = false;
+      MQuadrangle *q = gf->quadrangles[i];
+      for (int j = 0; j < q->getNumVertices(); j++) {
+        if (touched.find(q->getVertex(j)) != touched.end()) {
+          add_ = true;
+        }
+      }
+      if (add_) qs.insert(q);
+    }
+
+    for (std::set<MTriangle*>::iterator it = ts.begin(); it != ts.end(); ++it) {
+      for (int j = 0; j < (*it)->getNumVertices(); j++) {
+        touched.insert((*it)->getVertex(j));
+      }
+    }
+    for (std::set<MQuadrangle*>::iterator it = qs.begin(); it != qs.end(); ++it) {
+      for (int j = 0; j < (*it)->getNumVertices(); j++) {
+        touched.insert((*it)->getVertex(j));
+      }
+    }
+  }
+
+  for (int i = 0; i < gf->triangles.size(); i++) {
+    if (ts.find(gf->triangles[i]) == ts.end()) {
+      for (int j = 0; j < gf->triangles[i]->getNumVertices(); j++) {
+        if (touched.find(gf->triangles[i]->getVertex(j)) != touched.end()) not_touched.insert(gf->triangles[i]->getVertex(j));
+      }
+    }
+  }
+  for (int i = 0; i < gf->quadrangles.size(); i++) {
+    if (qs.find(gf->quadrangles[i]) == qs.end()) {
+      for (int j = 0; j < gf->quadrangles[i]->getNumVertices(); j++) {
+        if (touched.find(gf->quadrangles[i]->getVertex(j)) != touched.end()) not_touched.insert(gf->quadrangles[i]->getVertex(j));
+      }
+    }
+  }
+
+  gf->triangles.clear();
+  gf->quadrangles.clear();
+  gf->triangles.insert(gf->triangles.begin(), ts.begin(), ts.end());
+  gf->quadrangles.insert(gf->quadrangles.begin(), qs.begin(), qs.end());
+
+  printf("FILTER2D : AFTER FILTERING %d triangles %d quads remain %d nodes on the boundary\n", gf->triangles.size(), gf->quadrangles.size(), not_touched.size());
+  return not_touched;
+}
+
+std::set<MVertex*> filter3D(GRegion *gr, int nbLayers, double _qual)
+{
+  std::set<MVertex*> touched;
+  std::set<MVertex*> not_touched;
+  std::set<MTetrahedron *> ts;
+  for (int i = 0; i < gr->tetrahedra.size(); i++) {
+    MTetrahedron *t = gr->tetrahedra[i];
+    if (t->distoShapeMeasure() < _qual) {
+      ts.insert(t);
+      for (int j = 0; j < t->getNumVertices(); j++)
+        touched.insert(t->getVertex(j));
+    }
+    if (ts.size() == 1) break;
+  }
+
+  printf("FILTER3D : %d bad tets found among %6d\n", ts.size(), gr->tetrahedra.size());
+
+  // add layers of elements around bad elements
+  for (int layer = 0; layer < nbLayers; layer++) {
+    for (int i = 0; i < gr->tetrahedra.size(); i++) {
+      MTetrahedron *t = gr->tetrahedra[i];
+      bool add_ = false;
+      for (int j = 0; j < t->getNumVertices(); j++) {
+        if (touched.find(t->getVertex(j)) != touched.end()) {
+          add_ = true;
+        }
+      }
+      if (add_) ts.insert(t);
+    }
+
+    for (std::set<MTetrahedron*>::iterator it = ts.begin(); it != ts.end(); ++it) {
+      for (int j = 0; j < (*it)->getNumVertices(); j++) {
+        touched.insert((*it)->getVertex(j));
+      }
+    }
+  }
+
+  for (int i = 0; i < gr->tetrahedra.size(); i++) {
+    if (ts.find(gr->tetrahedra[i]) == ts.end()) {
+      for (int j = 0; j < gr->tetrahedra[i]->getNumVertices(); j++) {
+        if (touched.find(gr->tetrahedra[i]->getVertex(j)) != touched.end()) not_touched.insert(gr->tetrahedra[i]->getVertex(j));
+      }
+    }
+  }
+
+  gr->tetrahedra.clear();
+  gr->tetrahedra.insert(gr->tetrahedra.begin(), ts.begin(), ts.end());
+
+  printf("FILTER3D : AFTER FILTERING %d tets remain %d nodes on the boundary\n", gr->tetrahedra.size(), not_touched.size());
+  return not_touched;
+}
+
+void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
+{
+  int samples = 30;
+
+//  int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
+//  int method = Mesh::METHOD_FIXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
+//  int method = Mesh::METHOD_FIXBND | Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
+//  int method = Mesh::METHOD_FIXBND;
+//  int method = Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
+//  int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
+  int method = Mesh::METHOD_PROJJAC;
+
+  SVector3 BND(gm->bounds().max(), gm->bounds().min());
+  double SIZE = BND.norm();
+
+  if (p.dim == 2) {
+    for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
+
+      std::vector<MTriangle*> tris = (*itf)->triangles;
+      std::vector<MQuadrangle*> quas = (*itf)->quadrangles;
+      std::set<MVertex*> toFix = filter2D(*itf, p.nbLayers, p.BARRIER);
+      OptHOM *temp = new OptHOM(*itf, toFix, method);
+      (*itf)->triangles = tris;
+      (*itf)->quadrangles = quas;
+
+      double distMaxBND, distAvgBND, minJac, maxJac;
+      temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac);
+      std::ostringstream ossI;
+      ossI << "initial_" << (*itf)->tag() << ".msh";
+      temp->mesh.writeMSH(ossI.str().c_str());
+      printf("--> Mesh Loaded for GFace %d :  minJ %12.5E -- maxJ %12.5E\n", (*itf)->tag(), minJac, maxJac);
+      if (distMaxBND < 1.e-8 * SIZE && minJac > p.BARRIER) continue;
+      int result = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax);
+      temp->mesh.updateGEntityPositions();
+      std::ostringstream ossF;
+      ossF << "final_" << (*itf)->tag() << ".msh";
+      temp->mesh.writeMSH(ossF.str().c_str());
+    }
+  }
+  else if (p.dim == 3) {
+    for (GModel::riter itr = gm->firstRegion(); itr != gm->lastRegion(); ++itr) {
+      std::vector<MTetrahedron*> tets = (*itr)->tetrahedra;
+      std::set<MVertex*> toFix;
+      filter3D(*itr, p.nbLayers, p.BARRIER);
+      OptHOM *temp = new OptHOM(*itr, toFix, method);
+      double distMaxBND, distAvgBND, minJac, maxJac;
+      temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac);
+      temp->mesh.writeMSH("initial.msh");
+      printf("--> Mesh Loaded for GRegion %d :  minJ %12.5E -- maxJ %12.5E\n", (*itr)->tag(), minJac, maxJac);
+      if (distMaxBND < 1.e-8 * SIZE && minJac > p.BARRIER) continue;
+      int result = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax);
+      temp->mesh.updateGEntityPositions();
+      (*itr)->tetrahedra = tets;
+      temp->mesh.writeMSH("final.msh");
+    }
+  }  
+}

contrib/HighOrderMeshOptimizer/OptHomRun.h

+#ifndef _OPTHOMRUN_H_
+#define _OPTHOMRUN_H_
+
+class GModel;
+
+struct OptHomParameters {
+  // INPUT ------> 
+  double STOP ; // stop criterion
+  double BARRIER ; // minimum scaled jcaobian
+  double weightFixed ; // weight of the energy for fixed nodes
+  double weightFree ; // weight of the energy for free nodes
+  int nbLayers ; // number of layers taken around a bad element
+  int dim ; // which dimension to optimize
+  int itMax ; // max number of iterations in the optimization process
+  double TMAX ; // max CPU time allowed
+  bool onlyVisible ; // apply optimization to visible entities ONLY
+  // OUTPUT ------> 
+  int SUCCESS ; // 0 --> success , 1 --> Not converged
+  double minJac, maxJac; // after optimization, range of jacobians
+  double DT; // Time for optimization
+  
+  OptHomParameters () 
+  // default values    
+  : STOP (0.01) , BARRIER (0.1) , weightFixed (10.),  weightFree (1.e-3),
+    nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true)
+  {
+  }
+};
+
+void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p);
+
+
+#endif

contrib/HighOrderMeshOptimizer/ParamCoord.cpp

+#include <iostream>
+#include "GFace.h"
+#include "GEdge.h"
+#include "MVertex.h"
+#include "ParamCoord.h"
+
+
+
+ParamCoordSurf::ParamCoordSurf(GEntity *ge)
+{
+  if ((ge->dim() == 2) && (ge->geomType() != GEntity::DiscreteSurface)) _gf = static_cast<GFace*>(ge);
+  else std::cout << "ERROR: Using surface parametric coordinates with non-surface geometric entity" << std::endl;
+}
+
+
+
+SPoint3 ParamCoordSurf::getUvw(MVertex* vert)
+{
+  SPoint2 p;
+  reparamMeshVertexOnFace(vert,_gf,p);
+  return SPoint3(p[0],p[1],0.);
+}
+
+
+
+SPoint3 ParamCoordSurf::uvw2Xyz(MVertex* vert, const SPoint3 &uvw)
+{
+  GPoint gp = _gf->point(uvw[0],uvw[1]);
+  return SPoint3(gp.x(),gp.y(),gp.z());
+}
+
+
+
+void ParamCoordSurf::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw)
+{
+  Pair<SVector3,SVector3> der = _gf->firstDer(SPoint2(uvw[0],uvw[1]));
+  gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() + gXyz.z() * der.first().z();
+  gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() + gXyz.z() * der.second().z();
+}
+
+
+
+void ParamCoordSurf::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
+{
+  Pair<SVector3,SVector3> der = _gf->firstDer(SPoint2(uvw[0],uvw[1]));
+  std::vector<SPoint3>::iterator itUvw=gUvw.begin();
+  for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+    (*itUvw)[0] = itXyz->x() * der.first().x() + itXyz->y() * der.first().y() + itXyz->z() * der.first().z();
+    (*itUvw)[1] = itXyz->x() * der.second().x() + itXyz->y() * der.second().y() + itXyz->z() * der.second().z();
+    itUvw++;
+  }
+}
+
+
+
+SPoint3 ParamCoordParent::getUvw(MVertex* vert)
+{
+
+  GEntity *ge = vert->onWhat();
+  if ((ge->geomType() == GEntity::DiscreteCurve) || (ge->geomType() == GEntity::DiscreteSurface))
+    std::cout << "ERROR: using parent coordinates on discrete curve or surface" << std::endl;
+
+  switch (ge->dim()) {
+  case 1: {
+    SPoint3 p(0.,0.,0.);
+    reparamMeshVertexOnEdge(vert,static_cast<GEdge*>(ge),p[0]);
+    return p;
+    break;
+  }
+  case 2: {
+    SPoint2 p;
+    reparamMeshVertexOnFace(vert,static_cast<GFace*>(ge),p);
+    return SPoint3(p[0],p[1],0.);
+    break;
+  }
+  case 3: {
+    return vert->point();
+    break;
+  }
+  }
+
+}
+
+
+
+SPoint3 ParamCoordParent::uvw2Xyz(MVertex* vert, const SPoint3 &uvw)
+{
+
+  GEntity *ge = vert->onWhat();
+
+  switch (ge->dim()) {
+  case 1: {
+    GPoint gp = static_cast<GEdge*>(ge)->point(uvw[0]);
+    return SPoint3(gp.x(),gp.y(),gp.z());
+    break;
+  }
+  case 2: {
+    GPoint gp = static_cast<GFace*>(ge)->point(uvw[0],uvw[1]);
+    return SPoint3(gp.x(),gp.y(),gp.z());
+    break;
+  }
+  case 3: {
+    return uvw;
+    break;
+  }
+  }
+
+}
+
+
+
+void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw)
+{
+
+  GEntity *ge = vert->onWhat();
+
+  switch (ge->dim()) {
+  case 1: {
+    SVector3 der = static_cast<GEdge*>(ge)->firstDer(uvw[0]);
+    gUvw[0] = gXyz.x() * der.x() + gXyz.y() * der.y() + gXyz.z() * der.z();
+    break;
+  }
+  case 2: {
+    Pair<SVector3,SVector3> der = static_cast<GFace*>(ge)->firstDer(SPoint2(uvw[0],uvw[1]));
+    gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() + gXyz.z() * der.first().z();
+    gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() + gXyz.z() * der.second().z();
+    break;
+  }
+  case 3: {
+    gUvw = gXyz;
+    break;
+  }
+  }
+
+}
+
+
+
+void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
+{
+
+  GEntity *ge = vert->onWhat();
+
+  switch (ge->dim()) {
+  case 1: {
+    SVector3 der = static_cast<GEdge*>(ge)->firstDer(uvw[0]);
+    std::vector<SPoint3>::iterator itUvw=gUvw.begin();
+    for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+      (*itUvw)[0] = itXyz->x() * der.x() + itXyz->y() * der.y() + itXyz->z() * der.z();
+      itUvw++;
+    }
+    break;
+  }
+  case 2: {
+    Pair<SVector3,SVector3> der = static_cast<GFace*>(ge)->firstDer(SPoint2(uvw[0],uvw[1]));
+    std::vector<SPoint3>::iterator itUvw=gUvw.begin();
+    for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+      (*itUvw)[0] = itXyz->x() * der.first().x() + itXyz->y() * der.first().y() + itXyz->z() * der.first().z();
+      (*itUvw)[1] = itXyz->x() * der.second().x() + itXyz->y() * der.second().y() + itXyz->z() * der.second().z();
+      itUvw++;
+    }
+    break;
+  }
+  case 3: {
+    std::vector<SPoint3>::iterator itUvw=gUvw.begin();
+    for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+      *itUvw = *itXyz;
+      itUvw++;
+    }
+    break;
+  }
+  }
+
+}

contrib/HighOrderMeshOptimizer/ParamCoord.h

+#ifndef _PARAMCOORD_H_
+#define _PARAMCOORD_H_
+
+
+
+class ParamCoord
+{
+
+public:
+
+  virtual int nCoord(MVertex* vert) = 0;                                                 // Number of parametric coordinates for vertex
+  virtual SPoint3 getUvw(MVertex* vert) = 0;                                             // Get parametric coordinates of vertex
+  virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw) = 0;                        // Calculate physical coordinates from parametric coordinates of vertex
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) = 0; // Calculate derivatives w.r.t parametric coordinates
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw) = 0; // Calculate derivatives w.r.t parametric coordinates
+  virtual ~ParamCoord() {}
+
+};
+
+
+
+class ParamCoordPhys3D : public ParamCoord
+{
+
+public:
+
+  int nCoord(MVertex* vert) { return 3; }
+  virtual SPoint3 getUvw(MVertex* vert) { return vert->point(); }
+  virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw) { return uvw; }
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) { gUvw = gXyz; }
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
+  {
+    std::vector<SPoint3>::iterator itUvw=gUvw.begin();
+    for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+      *itUvw = *itXyz;
+      itUvw++;
+    }
+  }
+
+};
+
+
+
+class ParamCoordPhys2D : public ParamCoord
+{
+
+public:
+
+  int nCoord(MVertex* vert) { return 2; }
+  virtual SPoint3 getUvw(MVertex* vert) { return vert->point(); }
+  virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw) { return uvw; }
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) { gUvw = gXyz; }
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
+  {
+    std::vector<SPoint3>::iterator itUvw=gUvw.begin();
+    for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+      *itUvw = *itXyz;
+      itUvw++;
+    }
+  }
+
+};
+
+
+
+class ParamCoordSurf : public ParamCoord
+{
+
+public:
+
+  ParamCoordSurf(GEntity *ge);
+  int nCoord(MVertex* vert) { return 2; }
+  virtual SPoint3 getUvw(MVertex* vert);
+  virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw);
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw);
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw);
+
+private:
+
+  GFace *_gf;
+
+};
+
+
+
+class ParamCoordParent : public ParamCoord
+{
+
+public:
+
+  int nCoord(MVertex* vert) { return vert->onWhat()->dim(); }
+  virtual SPoint3 getUvw(MVertex* vert);
+  virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw);
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw);
+  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw);
+
+};
+
+
+
+#endif