diff --git a/CMakeLists.txt b/CMakeLists.txt
index de5f29e34a6fa86898c48d35dbe25fb3a5bfb7cb..af556300d5df1b17c8873506e6d4e1933c102f9a 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index 142b217e0175cab07112735f5af503953a146dc7..cf40d0e66b6cfabc65731c7ba67827ef7dbbf6ad 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/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();
diff --git a/contrib/HighOrderMeshOptimizer/CMakeLists.txt b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..31963a192265dadb4c534c86254aa28de27512b2
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
@@ -0,0 +1,14 @@
+# 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}")
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ce38c6b41dee62e436e49040a2a1f29eb1496519
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -0,0 +1,267 @@
+#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
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
new file mode 100644
index 0000000000000000000000000000000000000000..575a5eb3d3037c1050ee14c9ea7908c3433988a9
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -0,0 +1,85 @@
+#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
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a4249443864900a6c8daf6a7384bac04613b7754
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -0,0 +1,460 @@
+#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);
+
+}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..6a9067d98ed956fe63e3d053f3265f1fd870abf1
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -0,0 +1,113 @@
+#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
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b0d03813eb10b64aaee7c2a5c7309a60645efbd5
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -0,0 +1,228 @@
+//#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");
+ }
+ }
+}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
new file mode 100644
index 0000000000000000000000000000000000000000..10ad8a3f86390d0327571c97eedfde2ee8b4d399
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -0,0 +1,33 @@
+#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
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ccaba35d92abcfb79de17512b383fe578514cea3
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
@@ -0,0 +1,174 @@
+#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;
+ }
+ }
+
+}
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.h b/contrib/HighOrderMeshOptimizer/ParamCoord.h
new file mode 100644
index 0000000000000000000000000000000000000000..4f70573e9aeddeeb18e695d24f9cc9d48c158562
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.h
@@ -0,0 +1,101 @@
+#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