diff --git a/CMakeLists.txt b/CMakeLists.txt
index bc25e6bb451aec142fff5fb10180ffa462ef1113..97b9e3b0ec634458b5f8646d7c39a09ed938a3b5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -151,6 +151,8 @@ set(GMSH_API
contrib/MeshOptimizer/MeshOptCommon.h contrib/MeshOptimizer/MeshOptimizer.h
contrib/MeshOptimizer/MeshOptObjContribFunc.h contrib/MeshOptimizer/MeshOptObjContrib.h
contrib/MeshOptimizer/MeshOptObjectiveFunction.h contrib/MeshOptimizer/MeshOptVertexCoord.h
+ contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h
+ contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
contrib/MathEx/mathex.h)
get_property(IAMCHILD DIRECTORY PROPERTY PARENT_DIRECTORY)
@@ -645,6 +647,8 @@ if(ENABLE_OPTHOM)
include_directories(contrib/HighOrderMeshOptimizer)
add_subdirectory(contrib/MeshOptimizer)
include_directories(contrib/MeshOptimizer)
+ add_subdirectory(contrib/MeshQualityOptimizer)
+ include_directories(contrib/MeshQualityOptimizer)
set_config_option(HAVE_OPTHOM "OptHom")
endif(ENABLE_OPTHOM)
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index ea9615bf1255aeeefb0af65b8eb4ff62c92f2d73..37bd90050bc09511bf2a5a66c9a4d7b46898f2fe 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -295,6 +295,16 @@ double qmTriangle::angles(MTriangle *e)
}
+void qmTriangle::NCJRange(const MTriangle *el, double &valMin, double &valMax)
+{
+ fullVector<double> ncj(3);
+ const MVertex *v0 = el->getVertex(0), *v1 = el->getVertex(1), *v2 = el->getVertex(2);
+ NCJ(v0->x(), v0->y(), v0->z(), v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z(), ncj);
+ valMin = *std::min_element(ncj.getDataPtr(), ncj.getDataPtr()+3);
+ valMax = *std::max_element(ncj.getDataPtr(), ncj.getDataPtr()+3);
+}
+
+
void qmTriangle::NCJ(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2,
@@ -466,6 +476,18 @@ double qmQuadrangle::angles(MQuadrangle *e)
}
+void qmQuadrangle::NCJRange(const MQuadrangle *el, double &valMin, double &valMax)
+{
+ fullVector<double> ncj(4);
+ const MVertex *v0 = el->getVertex(0), *v1 = el->getVertex(1),
+ *v2 = el->getVertex(2), *v3 = el->getVertex(3);
+ NCJ(v0->x(), v0->y(), v0->z(), v1->x(), v1->y(), v1->z(),
+ v2->x(), v2->y(), v2->z(), v3->x(), v3->y(), v3->z(), ncj);
+ valMin = *std::min_element(ncj.getDataPtr(), ncj.getDataPtr()+4);
+ valMax = *std::max_element(ncj.getDataPtr(), ncj.getDataPtr()+4);
+}
+
+
void qmQuadrangle::NCJ(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2,
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index 8086d20856051492719b78ec61de90f5402f64dc..324f9745a249c6afcbe7d8d08b7981b621a997b9 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -33,6 +33,8 @@ public:
static double eta(MTriangle *el);
static double angles(MTriangle *e);
static double minNCJ(const MTriangle *e);
+ static void NCJRange(const MTriangle *e, double &valMin, double &valMax);
+ static inline int numNCJVal() { return 3; }
static void NCJ(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2,
@@ -51,6 +53,8 @@ public:
static double eta(MQuadrangle *el);
static double angles(MQuadrangle *e);
static double minNCJ(const MQuadrangle *e);
+ static void NCJRange(const MQuadrangle *e, double &valMin, double &valMax);
+ static inline int numNCJVal() { return 4; }
static void NCJ(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2,
@@ -93,6 +97,8 @@ public:
const double &x4, const double &y4, const double &z4,
double *volume = 0);
static double minNCJ(const MTetrahedron *e);
+// static void NCJRange(const MTetrahedron *e, double &valMin, double &valMax);
+ static inline int numNCJVal() { return 4; }
static void NCJ(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2,
@@ -110,6 +116,7 @@ class qmPrism
{
public:
static double minNCJ(const MPrism *el);
+ static inline int numNCJVal() { return 6; }
// static void NCJ(const double &x0, const double &y0, const double &z0,
// const double &x1, const double &y1, const double &z1,
// const double &x2, const double &y2, const double &z2,
@@ -123,6 +130,7 @@ class qmHexahedron
{
public:
static double angles(MHexahedron *el);
+ static inline int numNCJVal() { return 8; }
// static void NCJ(const double &x0, const double &y0, const double &z0,
// const double &x1, const double &y1, const double &z1,
// const double &x2, const double &y2, const double &z2,
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 8cc3573451f308d0f85ed6afba3be70934b6cbeb..94b970e882ace3609b6a80c130ee0546d4cdf2c5 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -758,7 +758,7 @@ struct HOPatchDefParameters : public MeshOptParameters::PatchDefParameters
HOPatchDefParameters(const OptHomParameters &p);
virtual ~HOPatchDefParameters() {}
virtual double elBadness(MElement *el);
- virtual double maxDistance(const MElement *el);
+ virtual double maxDistance(MElement *el);
private:
double jacMin, jacMax;
double distanceFactor;
@@ -802,7 +802,7 @@ double HOPatchDefParameters::elBadness(MElement *el) {
}
-double HOPatchDefParameters::maxDistance(const MElement *el) {
+double HOPatchDefParameters::maxDistance(MElement *el) {
return distanceFactor * el->maxDistToStraight();
}
diff --git a/contrib/MeshOptimizer/MeshOptCommon.h b/contrib/MeshOptimizer/MeshOptCommon.h
index 80fd94edaeb51806ed740562133ea0633641116e..9f895e84f9c5ede16ae4f6847aa42dbcae194a1e 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.h
+++ b/contrib/MeshOptimizer/MeshOptCommon.h
@@ -56,7 +56,7 @@ struct MeshOptParameters { // Parameters controlling
bool weakMerge; // If connected strategy: weak or strong merging of patches
};
virtual double elBadness(MElement *el) = 0; // Pointer to function returning "badness" of element (for patch creation)
- virtual double maxDistance(const MElement *el) = 0; // Pointer to function checking the patch distance criterion for a given bad element
+ virtual double maxDistance(MElement *el) = 0; // Pointer to function checking the patch distance criterion for a given bad element
};
int dim ; // Which dimension to optimize
bool onlyVisible ; // Apply optimization to visible entities ONLY
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index 0577a2525783f9f93454297eed8bc8d6918a9a4e..cfed7a6823dc30c971bf9b53a1d46d1239c27294 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -35,6 +35,7 @@
#include "MTetrahedron.h"
#include "BasisFactory.h"
#include "OptHomIntegralBoundaryDist.h"
+#include "qualityMeasures.h"
#include "MeshOptPatch.h"
@@ -496,3 +497,98 @@ bool Patch::bndDistAndGradients(int iEl, double &f, std::vector<double> &gradF,
}
return edgeFound;
}
+
+
+void Patch::initNCJ()
+{
+ // Initialize _nBezEl
+ if (_nNCJEl.empty()) {
+ _nNCJEl.resize(nEl());
+ for (int iEl=0; iEl<nEl(); iEl++)
+ switch(_el[iEl]->getType()) { // TODO: Complete with other types?
+ case TYPE_TRI: _nNCJEl[iEl] = qmTriangle::numNCJVal(); break;
+ case TYPE_QUA: _nNCJEl[iEl] = qmQuadrangle::numNCJVal(); break;
+ }
+ }
+}
+
+
+void Patch::NCJ(int iEl, std::vector<double> &NCJ)
+{
+ const int &numNCJVal = _nNCJEl[iEl];
+ const int &numNodes = _nNodEl[iEl];
+
+ std::vector<int> &iVEl = _el2V[iEl];
+ fullVector<double> val(numNCJVal);
+ switch(_el[iEl]->getType()) { // TODO: Complete with other types?
+ case TYPE_TRI: {
+ const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], &iV2 = _el2V[iEl][2];
+ qmTriangle::NCJ(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(),
+ _xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(),
+ _xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(), val);
+ break;
+ }
+ case TYPE_QUA: {
+ const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1],
+ &iV2 = _el2V[iEl][2], &iV3 = _el2V[iEl][3];
+ qmQuadrangle::NCJ(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(),
+ _xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(),
+ _xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(),
+ _xyz[iV3].x(), _xyz[iV3].y(), _xyz[iV3].z(), val);
+ break;
+ }
+ }
+
+ for (int l = 0; l < numNCJVal; l++) NCJ[l] = val(l);
+}
+
+
+void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<double> &gNCJ)
+{
+ const int &numNCJVal = _nNCJEl[iEl];
+ const int &numNodes = _nNodEl[iEl];
+
+ std::vector<int> &iVEl = _el2V[iEl];
+ fullMatrix<double> gradVal(numNCJVal, 3*numNodes+1);
+ switch(_el[iEl]->getType()) { // TODO: Complete with other types?
+ case TYPE_TRI: {
+ const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], &iV2 = _el2V[iEl][2];
+ qmTriangle::NCJAndGradients(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(),
+ _xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(),
+ _xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(), gradVal);
+ break;
+ }
+ case TYPE_QUA: {
+ const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1],
+ &iV2 = _el2V[iEl][2], &iV3 = _el2V[iEl][3];
+ qmQuadrangle::NCJAndGradients(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(),
+ _xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(),
+ _xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(),
+ _xyz[iV3].x(), _xyz[iV3].y(), _xyz[iV3].z(), gradVal);
+ break;
+ }
+ }
+
+ // NCJ
+ for (int l = 0; l < numNCJVal; l++) NCJ[l] = gradVal(l,3*numNodes);
+
+ // Gradients of the NCJ
+ int iPC = 0;
+ std::vector<SPoint3> gXyzV(numNCJVal);
+ std::vector<SPoint3> gUvwV(numNCJVal);
+ for (int i = 0; i < numNodes; i++) {
+ int &iFVi = _el2FV[iEl][i];
+ if (iFVi >= 0) {
+ for (int l = 0; l < numNCJVal; l++)
+ gXyzV [l] = SPoint3(gradVal(l,i+0*numNodes), gradVal(l,i+1*numNodes),
+ gradVal(l,i+2*numNodes));
+ _coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV);
+ for (int l = 0; l < numNCJVal; l++) {
+ gNCJ[indGNCJ(iEl,l,iPC)] = gUvwV[l][0];
+ if (_nPCFV[iFVi] >= 2) gNCJ[indGNCJ(iEl,l,iPC+1)] = gUvwV[l][1];
+ if (_nPCFV[iFVi] == 3) gNCJ[indGNCJ(iEl,l,iPC+2)] = gUvwV[l][2];
+ }
+ iPC += _nPCFV[iFVi];
+ }
+ }
+}
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index 209f173c0a5b30d7df0625a3d08d82db79de7a9b..5a2f0613860caedf002d24a4a5a6301e38540a5d 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/contrib/MeshOptimizer/MeshOptPatch.h
@@ -92,6 +92,13 @@ public:
void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
bool bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps);
+ // Mesh quality
+ inline const int &nNCJEl(int iEl) { return _nNCJEl[iEl]; }
+ inline int indGNCJ(int iEl, int l, int iPC) { return iPC*_nNCJEl[iEl]+l; }
+ void initNCJ();
+ void NCJ(int iEl, std::vector<double> &NCJ);
+ void NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<double> &gNCJ);
+
private:
// Mesh entities and variables
@@ -133,11 +140,14 @@ private:
double _invLengthScaleSq; // Square inverse of a length for node displacement scaling
// High-order: scaled Jacobian and metric measures
- std::vector<int> _nBezEl; // Number of Bezier poly. and for an el.
+ std::vector<int> _nBezEl; // Number of Bezier poly. for an el.
std::vector<fullMatrix<double> > _scaledNormEl; // Normals to 2D elements for Jacobian regularization and scaling
std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements
// void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl);
void calcScaledNormalEl2D(int iEl);
+
+ // Mesh quality
+ std::vector<int> _nNCJEl; // Number of NCJ values for an el.
};
diff --git a/contrib/MeshQualityOptimizer/CMakeLists.txt b/contrib/MeshQualityOptimizer/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d958b468cce1bf9c6a38b7113413ef5adb099508
--- /dev/null
+++ b/contrib/MeshQualityOptimizer/CMakeLists.txt
@@ -0,0 +1,11 @@
+# Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+#
+# See the LICENSE.txt file for license information. Please report all
+# bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+set(SRC
+ MeshQualityOptimizer.cpp
+)
+
+file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.hpp)
+append_gmsh_src(contrib/MeshQualityOptimizer "${SRC};${HDR}")
diff --git a/contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h b/contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h
new file mode 100644
index 0000000000000000000000000000000000000000..369ff30b54a145b4e7febc3b1b4c9659b8796861
--- /dev/null
+++ b/contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h
@@ -0,0 +1,97 @@
+// TODO: Copyright
+
+#ifndef _MESHQUALITYCONTRIBNCJ_H_
+#define _MESHQUALITYCONTRIBNCJ_H_
+
+#include "MeshOptPatch.h"
+#include "MeshOptObjContrib.h"
+
+
+template<class FuncType>
+class ObjContribNCJ : public ObjContrib, public FuncType
+{
+public:
+ ObjContribNCJ(double weight);
+ virtual ~ObjContribNCJ() {}
+ virtual ObjContrib *copy() const;
+ virtual void initialize(Patch *mesh);
+ virtual bool fail() { return _min <= 0.; }
+ virtual bool addContrib(double &Obj, alglib::real_1d_array &gradObj);
+ virtual void updateParameters() { FuncType::updateParameters(_min, _max); }
+ virtual bool targetReached() { return FuncType::targetReached(_min, _max); }
+ virtual bool stagnated() { return FuncType::stagnated(_min, _max); }
+ virtual void updateMinMax();
+
+protected:
+ Patch *_mesh;
+ double _weight;
+};
+
+
+template<class FuncType>
+ObjContribNCJ<FuncType>::ObjContribNCJ(double weight) :
+ ObjContrib("NCJ", FuncType::getNamePrefix()+"NCJ"),
+ _mesh(0), _weight(weight)
+{
+}
+
+
+template<class FuncType>
+ObjContrib *ObjContribNCJ<FuncType>::copy() const
+{
+ return new ObjContribNCJ<FuncType>(*this);
+}
+
+
+template<class FuncType>
+void ObjContribNCJ<FuncType>::initialize(Patch *mesh)
+{
+ _mesh = mesh;
+ _mesh->initNCJ();
+ updateMinMax();
+ FuncType::initialize(_min, _max);
+}
+
+
+template<class FuncType>
+bool ObjContribNCJ<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
+{
+ _min = BIGVAL;
+ _max = -BIGVAL;
+
+ for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ std::vector<double> NCJ(_mesh->nNCJEl(iEl)); // Scaled Jacobians
+ std::vector<double> gNCJ(_mesh->nNCJEl(iEl)*_mesh->nPCEl(iEl)); // Gradients of scaled Jacobians
+ _mesh->NCJAndGradients(iEl, NCJ, gNCJ);
+ for (int l = 0; l < _mesh->nNCJEl(iEl); l++) { // Add contribution for each Bezier coeff.
+ Obj += _weight * FuncType::compute(NCJ[l]);
+ const double dfact = _weight * FuncType::computeDiff(NCJ[l]);
+ for (int iPC = 0; iPC < _mesh->nPCEl(iEl); iPC++)
+ gradObj[_mesh->indPCEl(iEl, iPC)] += dfact * gNCJ[_mesh->indGNCJ(iEl, l, iPC)];
+ _min = std::min(_min, NCJ[l]);
+ _max = std::max(_max, NCJ[l]);
+ }
+ }
+
+ return true;
+}
+
+
+template<class FuncType>
+void ObjContribNCJ<FuncType>::updateMinMax()
+{
+ _min = BIGVAL;
+ _max = -BIGVAL;
+
+ for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ std::vector<double> NCJ(_mesh->nNCJEl(iEl)); // Scaled Jacobians
+ _mesh->NCJ(iEl, NCJ);
+ for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Check each Bezier coeff.
+ _min = std::min(_min, NCJ[l]);
+ _max = std::max(_max, NCJ[l]);
+ }
+ }
+}
+
+
+#endif /* _MESHQUALITYCONTRIBNCJ_H_ */
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c0c502457c2f62b4674a13ffdbae02af87318110
--- /dev/null
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
@@ -0,0 +1,113 @@
+// TODO: Copyright
+
+//#include "GModel.h"
+#include "GEntity.h"
+#include "MElement.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "qualityMeasures.h"
+#include "MeshOptCommon.h"
+#include "MeshOptObjContribFunc.h"
+#include "MeshOptObjContrib.h"
+#include "MeshOptObjContribScaledNodeDispSq.h"
+#include "MeshQualityObjContribNCJ.h"
+#include "MeshOptimizer.h"
+#include "MeshQualityOptimizer.h"
+
+
+struct QualPatchDefParameters : public MeshOptParameters::PatchDefParameters
+{
+ QualPatchDefParameters(const MeshQualOptParameters &p);
+ virtual ~QualPatchDefParameters() {}
+ virtual double elBadness(MElement *el);
+ virtual double maxDistance(MElement *el);
+private:
+ double NCJMin, NCJMax;
+ double distanceFactor;
+};
+
+
+QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
+{
+ NCJMin = p.minTargetNCJ;
+ NCJMax = (p.maxTargetNCJ > 0.) ? p.maxTargetNCJ : 1.e300;
+ strategy = (p.strategy == 1) ? MeshOptParameters::STRAT_ONEBYONE :
+ MeshOptParameters::STRAT_CONNECTED;
+ minLayers = (p.dim == 3) ? 1 : 0;
+ maxLayers = p.nbLayers;
+ distanceFactor = p.distanceFactor;
+ if (strategy == MeshOptParameters::STRAT_CONNECTED)
+ weakMerge = (p.strategy == 2);
+ else {
+ maxAdaptPatch = p.maxAdaptBlob;
+ maxLayersAdaptFact = p.adaptBlobLayerFact;
+ distanceAdaptFact = p.adaptBlobDistFact;
+ }
+}
+
+
+double QualPatchDefParameters::elBadness(MElement *el) {
+ double valMin, valMax;
+ switch(el->getType()) { // TODO: Complete with other types?
+ case TYPE_TRI: {
+ qmTriangle::NCJRange(static_cast<MTriangle*>(el), valMin, valMax);
+ break;
+ }
+ case TYPE_QUA: {
+ qmQuadrangle::NCJRange(static_cast<MQuadrangle*>(el), valMin, valMax);
+ break;
+ }
+ }
+ double badness = std::min(valMin-NCJMin, 0.) + std::min(NCJMin-valMax, 0.);
+ return badness;
+}
+
+
+double QualPatchDefParameters::maxDistance(MElement *el) {
+ return distanceFactor * el->getOuterRadius();
+}
+
+
+void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
+{
+ Msg::StatusBar(true, "Optimizing mesh quality...");
+
+ MeshOptParameters par;
+ par.dim = p.dim;
+ par.onlyVisible = p.onlyVisible;
+ par.fixBndNodes = p.fixBndNodes;
+ QualPatchDefParameters patchDef(p);
+ par.patchDef = &patchDef;
+ par.optDisplay = 30;
+ par.verbose = 4;
+
+ ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree);
+ ObjContribNCJ<ObjContribFuncBarrierMovMin> minNCJBarFunc(1.);
+ minNCJBarFunc.setTarget(p.minTargetNCJ, 1.);
+ ObjContribNCJ<ObjContribFuncBarrierFixMinMovMax> minMaxNCJBarFunc(1.);
+ minMaxNCJBarFunc.setTarget(p.maxTargetNCJ, 1.);
+
+ MeshOptParameters::PassParameters minJacPass;
+ minJacPass.barrierIterMax = p.optPassMax;
+ minJacPass.optIterMax = p.itMax;
+ minJacPass.contrib.push_back(&nodeDistFunc);
+ minJacPass.contrib.push_back(&minNCJBarFunc);
+ par.pass.push_back(minJacPass);
+
+ if (p.maxTargetNCJ > 0.) {
+ MeshOptParameters::PassParameters minMaxJacPass;
+ minMaxJacPass.barrierIterMax = p.optPassMax;
+ minMaxJacPass.optIterMax = p.itMax;
+ minMaxJacPass.contrib.push_back(&nodeDistFunc);
+ minMaxJacPass.contrib.push_back(&minMaxNCJBarFunc);
+ par.pass.push_back(minMaxJacPass);
+ }
+
+ meshOptimizer(gm, par);
+
+ p.CPU = par.CPU;
+ p.minNCJ = minMaxNCJBarFunc.getMin();
+ p.maxNCJ = minMaxNCJBarFunc.getMax();
+
+ Msg::StatusBar(true, "Done optimizing high order mesh (%g s)", p.CPU);
+}
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
new file mode 100644
index 0000000000000000000000000000000000000000..287b3a418559943d25e7099178616392a1877d39
--- /dev/null
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
@@ -0,0 +1,68 @@
+// Copyright (C) 2013 ULg-UCL
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use, copy,
+// modify, merge, publish, distribute, and/or sell copies of the
+// Software, and to permit persons to whom the Software is furnished
+// to do so, provided that the above copyright notice(s) and this
+// permission notice appear in all copies of the Software and that
+// both the above copyright notice(s) and this permission notice
+// appear in supporting documentation.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR
+// ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY
+// DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
+// WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
+// ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
+// OF THIS SOFTWARE.
+//
+// Please report all bugs and problems to the public mailing list
+// <gmsh@geuz.org>.
+//
+// Contributors: Thomas Toulorge, Jonathan Lambrechts
+
+#ifndef _MESHQUALITYOPTIMIZER_H_
+#define _MESHQUALITYOPTIMIZER_H_
+
+class GModel;
+
+struct MeshQualOptParameters {
+ double minTargetNCJ, maxTargetNCJ;
+ 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
+ int optPassMax ; // max number of optimization passes
+ bool onlyVisible ; // apply optimization to visible entities ONLY
+ double distanceFactor; // filter elements such that no elements further away
+ // than DistanceFactor times the max distance to
+ // straight sided version of an element are optimized
+ bool fixBndNodes; // how jacobians are computed and if points can move on boundaries
+ int strategy; // 0 = connected blobs, 1 = adaptive one-by-one
+ int maxAdaptBlob; // Max. nb. of blob adaptation interations
+ int adaptBlobLayerFact; // Growth factor in number of layers for blob adaptation
+ double adaptBlobDistFact; // Growth factor in distance factor for blob adaptation
+
+ int SUCCESS ; // 0 --> success , 1 --> Not converged
+ double minNCJ, maxNCJ; // after optimization, range of jacobians
+ double CPU; // Time for optimization
+
+ MeshQualOptParameters ()
+ : minTargetNCJ(0.5), maxTargetNCJ(1.5), weightFixed(1000.),
+ weightFree (1.), nbLayers (6) , dim(3) , itMax(300), optPassMax(50),
+ onlyVisible(true), distanceFactor(12), fixBndNodes(false), strategy(0),
+ maxAdaptBlob(3), adaptBlobLayerFact(2.), adaptBlobDistFact(2.)
+ {
+ }
+};
+
+void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p);
+
+#endif
diff --git a/wrappers/gmshpy/CMakeLists.txt b/wrappers/gmshpy/CMakeLists.txt
index 7a34421ddafebd817a4722ad8bb1266c3a80fb55..687dc25866416fe2e87fab0362cc8380d0da7b77 100644
--- a/wrappers/gmshpy/CMakeLists.txt
+++ b/wrappers/gmshpy/CMakeLists.txt
@@ -141,6 +141,7 @@ set (GMSH_PYTHON_EXTRA_INCLUDE
Mesh/meshGFaceLloyd.h
contrib/HighOrderMeshOptimizer/OptHomRun.h
contrib/HighOrderMeshOptimizer/OptHomElastic.h
+ contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
Post/PViewAsSimpleFunction.h
Post/PViewDataList.h
Post/PViewFactory.h
diff --git a/wrappers/gmshpy/gmshMesh.i b/wrappers/gmshpy/gmshMesh.i
index 8974de1a4d9fba4b7c334ed3c5fe2eeba429b854..361cfed3e9c4ef1aa06fb98b0df7cb4bda1579fb 100644
--- a/wrappers/gmshpy/gmshMesh.i
+++ b/wrappers/gmshpy/gmshMesh.i
@@ -18,6 +18,7 @@
#include "OptHomRun.h"
#include "OptHomElastic.h"
#include "OptHomFastCurving.h"
+ #include "MeshQualityOptimizer.h"
#endif
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
#include "meshPartition.h"
@@ -57,6 +58,7 @@ namespace std {
%include "OptHomRun.h"
%include "OptHomElastic.h"
%include "OptHomFastCurving.h"
+%include "MeshQualityOptimizer.h"
#endif
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
%include "meshPartition.h"