diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index 1af5b3fd32e28bf1f876231e3360a2bcad4b6e7c..c4f71eb0620b6abce8c24e29c8b36c2683e52430 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -269,7 +269,7 @@ void Patch::initScaledJac()
if ((_dim == 2) && _scaledNormEl.empty()) {
_scaledNormEl.resize(nEl());
// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
- for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, true);
+ for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, NS_INVNORM);
}
else if (_invStraightJac.empty()) {
_invStraightJac.resize(nEl(),1.);
@@ -293,7 +293,7 @@ void Patch::initMetricMin()
// TODO: Re-introduce normal to geometry?
//void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl)
-void Patch::calcNormalEl2D(int iEl, bool scale)
+void Patch::calcNormalEl2D(int iEl, NORMALSCALING scaling)
{
const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
@@ -320,11 +320,12 @@ void Patch::calcNormalEl2D(int iEl, bool scale)
// }
fullMatrix<double> uElNorm(1, 3);
- fullMatrix<double> &elNorm = scale ? _scaledNormEl[iEl] : uElNorm;
+ fullMatrix<double> &elNorm = (scaling == NS_INVNORM) ? _scaledNormEl[iEl] :
+ ((scaling == NS_UNIT) ? uElNorm : _metricNormEl[iEl]);
elNorm.resize(1, 3);
const double norm = jac->getPrimNormal2D(primNodesXYZ, elNorm);
- if (scale) {
- double factor = 1./norm;
+ if (scaling != 0) {
+ double factor = (scaling == NS_SQRTNORM) ? 1./norm : sqrt(norm);
// if (hasGeoNorm) {
// const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2);
// if (scal < 0.) factor = -factor;
@@ -520,7 +521,7 @@ void Patch::initNCJ()
if ((_dim == 2) && _unitNormEl.empty()) {
_unitNormEl.resize(nEl());
// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
- for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, false);
+ for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, NS_UNIT);
}
}
@@ -591,3 +592,68 @@ void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<doubl
}
}
}
+
+
+void Patch::initInvCond()
+{
+ // Initialize _nBezEl
+ if (_nBezEl.empty()) {
+ _nBezEl.resize(nEl());
+ for (int iEl=0; iEl<nEl(); iEl++)
+ _nBezEl[iEl] = _el[iEl]->getJacobianFuncSpace()->getNumJacNodes();
+ }
+
+ // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
+ // Jacobians of 3D elements
+ if ((_dim == 2) && _metricNormEl.empty()) {
+ _metricNormEl.resize(nEl());
+// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
+ for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, NS_SQRTNORM);
+ }
+}
+
+
+void Patch::invCondAndGradients(int iEl, std::vector<double> &invCond,
+ std::vector<double> &gInvCond)
+{
+ const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
+ const int &numJacNodes = _nBezEl[iEl];
+ const int &numMapNodes = _nNodEl[iEl];
+ fullMatrix<double> IDI(numJacNodes,3*numMapNodes+1);
+
+ // Coordinates of nodes
+ fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3);
+ for (int i = 0; i < numMapNodes; i++) {
+ int &iVi = _el2V[iEl][i];
+ nodesXYZ(i,0) = _xyz[iVi].x();
+ nodesXYZ(i,1) = _xyz[iVi].y();
+ nodesXYZ(i,2) = _xyz[iVi].z();
+ }
+
+ // Calculate Jacobian and gradients, scale if 3D (already scaled by
+ // regularization normals in 2D)
+ jacBasis->getInvCondAndGradients(nodesXYZ, _metricNormEl[iEl], IDI);
+
+ // Inverse condition number
+ for (int l = 0; l < numJacNodes; l++) invCond[l] = IDI(l, 3*numMapNodes);
+
+ // Gradients of the inverse condition number
+ int iPC = 0;
+ std::vector<SPoint3> gXyzV(numJacNodes);
+ std::vector<SPoint3> gUvwV(numJacNodes);
+ for (int i = 0; i < numMapNodes; i++) {
+ int &iFVi = _el2FV[iEl][i];
+ if (iFVi >= 0) {
+ for (int l = 0; l < numJacNodes; l++)
+ gXyzV[l] = SPoint3(IDI(l, i), IDI(l, i+numMapNodes),
+ IDI(l, i+2*numMapNodes));
+ _coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV);
+ for (int l = 0; l < numJacNodes; l++) {
+ gInvCond[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
+ if (_nPCFV[iFVi] >= 2) gInvCond[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+ if (_nPCFV[iFVi] == 3) gInvCond[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
+ }
+ iPC += _nPCFV[iFVi];
+ }
+ }
+}
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index 9ff8e026e5df44a837691ec0f8a200cf1efdca77..4269ba833d02510802bea41b289f0f77d8f86a5b 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/contrib/MeshOptimizer/MeshOptPatch.h
@@ -98,9 +98,13 @@ public:
void initNCJ();
void NCJ(int iEl, std::vector<double> &NCJ);
void NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<double> &gNCJ);
+ void initInvCond();
+ void invCondAndGradients(int iEl, std::vector<double> &invCond, std::vector<double> &ginvCond);
private:
+ enum NORMALSCALING { NS_UNIT, NS_INVNORM, NS_SQRTNORM };
+
// Mesh entities and variables
int _dim;
int _nPC; // Total nb. of parametric coordinates
@@ -144,11 +148,12 @@ private:
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 calcNormalEl2D(int iEl, bool scale);
+ void calcNormalEl2D(int iEl, NORMALSCALING scaling);
// Mesh quality
std::vector<int> _nNCJEl; // Number of NCJ values for an el.
std::vector<SVector3> _unitNormEl; // Normals to 2D elements for NCJ regularization
+ std::vector<fullMatrix<double> > _metricNormEl; // Normals to 2D elements for inverse conditioning computation
};
diff --git a/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h b/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
new file mode 100644
index 0000000000000000000000000000000000000000..ccdac14fcf1f1566c0a5f46dbca406dbdb12c807
--- /dev/null
+++ b/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
@@ -0,0 +1,98 @@
+// TODO: Copyright
+
+#ifndef _MESHQUALITYOBJCONTRIBINVCOND_H_
+#define _MESHQUALITYOBJCONTRIBINVCOND_H_
+
+#include "MeshOptPatch.h"
+#include "MeshOptObjContrib.h"
+
+
+template<class FuncType>
+class ObjContribInvCond : public ObjContrib, public FuncType
+{
+public:
+ ObjContribInvCond(double weight);
+ virtual ~ObjContribInvCond() {}
+ virtual ObjContrib *copy() const;
+ virtual void initialize(Patch *mesh);
+ virtual bool fail() { return false; }
+ 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>
+ObjContribInvCond<FuncType>::ObjContribInvCond(double weight) :
+ ObjContrib("InvCond", FuncType::getNamePrefix()+"InvCond"),
+ _mesh(0), _weight(weight)
+{
+}
+
+
+template<class FuncType>
+ObjContrib *ObjContribInvCond<FuncType>::copy() const
+{
+ return new ObjContribInvCond<FuncType>(*this);
+}
+
+
+template<class FuncType>
+void ObjContribInvCond<FuncType>::initialize(Patch *mesh)
+{
+ _mesh = mesh;
+ _mesh->initInvCond();
+ updateMinMax();
+ FuncType::initialize(_min, _max);
+}
+
+
+template<class FuncType>
+bool ObjContribInvCond<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
+{
+ _min = BIGVAL;
+ _max = -BIGVAL;
+
+ for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ std::vector<double> invCond(_mesh->nBezEl(iEl)); // Min. of Metric
+ std::vector<double> gInvCond(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
+ _mesh->invCondAndGradients(iEl, invCond, gInvCond);
+ for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Add contribution for each Bezier coeff.
+ Obj += _weight * FuncType::compute(invCond[l]);
+ const double dfact = _weight * FuncType::computeDiff(invCond[l]);
+ for (int iPC = 0; iPC < _mesh->nPCEl(iEl); iPC++)
+ gradObj[_mesh->indPCEl(iEl,iPC)] += dfact * gInvCond[_mesh->indGSJ(iEl, l, iPC)];
+ _min = std::min(_min, invCond[l]);
+ _max = std::max(_max, invCond[l]);
+ }
+ }
+
+ return true;
+}
+
+
+template<class FuncType>
+void ObjContribInvCond<FuncType>::updateMinMax()
+{
+ _min = BIGVAL;
+ _max = -BIGVAL;
+
+ for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ std::vector<double> invCond(_mesh->nBezEl(iEl)); // Min. of Metric
+ std::vector<double> dumGInvCond(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
+ _mesh->invCondAndGradients(iEl, invCond, dumGInvCond);
+ for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Add contribution for each Bezier coeff.
+ _min = std::min(_min, invCond[l]);
+ _max = std::max(_max, invCond[l]);
+ }
+ }
+}
+
+
+#endif /* _MESHQUALITYOBJCONTRIBINVCOND_H_ */
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
index a50d9eed61b1dd925e5246a4da63679b9b8e3805..c50d513767ae7ea9c2726e4501a221036e21a95f 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
@@ -11,6 +11,7 @@
#include "MeshOptObjContrib.h"
#include "MeshOptObjContribScaledNodeDispSq.h"
#include "MeshQualityObjContribNCJ.h"
+#include "MeshQualityObjContribInvCond.h"
#include "MeshOptimizer.h"
#include "MeshQualityOptimizer.h"
@@ -22,15 +23,17 @@ struct QualPatchDefParameters : public MeshOptParameters::PatchDefParameters
virtual double elBadness(MElement *el);
virtual double maxDistance(MElement *el);
private:
- double NCJMin, NCJMax;
+// double NCJMin, NCJMax;
+ double invCondMin;
double distanceFactor;
};
QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
{
- NCJMin = p.minTargetNCJ;
- NCJMax = (p.maxTargetNCJ > 0.) ? p.maxTargetNCJ : 1.e300;
+// NCJMin = p.minTargetNCJ;
+// NCJMax = (p.maxTargetNCJ > 0.) ? p.maxTargetNCJ : 1.e300;
+ invCondMin = p.minTargetInvCond;
strategy = (p.strategy == 1) ? MeshOptParameters::STRAT_ONEBYONE :
MeshOptParameters::STRAT_CONNECTED;
minLayers = (p.dim == 3) ? 1 : 0;
@@ -47,18 +50,25 @@ QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
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.);
+// 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;
+// }
+// }
+ const JacobianBasis *jac = el->getJacobianFuncSpace();
+ fullMatrix<double> nodesXYZ(el->getNumShapeFunctions(), 3);
+ el->getNodesCoord(nodesXYZ);
+ fullVector<double> invCond(jac->getNumJacNodes());
+ jac->getInvCond(nodesXYZ, invCond);
+ const double valMin = *std::min_element(invCond.getDataPtr(),
+ invCond.getDataPtr()+jac->getNumJacNodes());
+ double badness = std::min(valMin-invCondMin, 0.);
return badness;
}
@@ -82,33 +92,38 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
par.verbose = 4;
ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree);
- ObjContribNCJ<ObjContribFuncBarrierMovMin> minNCJBarFunc(1.);
- minNCJBarFunc.setTarget(p.minTargetNCJ, 1.);
-// minNCJBarFunc.setTarget(p.minTargetNCJ, 0.866025404);
- ObjContribNCJ<ObjContribFuncBarrierFixMinMovMax> minMaxNCJBarFunc(1.);
- minMaxNCJBarFunc.setTarget(p.maxTargetNCJ, 1.);
+// ObjContribNCJ<ObjContribFuncBarrierMovMin> minNCJBarFunc(1.);
+// minNCJBarFunc.setTarget(p.minTargetNCJ, 1.);
+//// minNCJBarFunc.setTarget(p.minTargetNCJ, 0.866025404);
+// ObjContribNCJ<ObjContribFuncBarrierFixMinMovMax> minMaxNCJBarFunc(1.);
+// minMaxNCJBarFunc.setTarget(p.maxTargetNCJ, 1.);
+ ObjContribInvCond<ObjContribFuncBarrierMovMin> minInvCondBarFunc(1.);
+ minInvCondBarFunc.setTarget(p.minTargetInvCond, 1.);
MeshOptParameters::PassParameters minJacPass;
minJacPass.barrierIterMax = p.optPassMax;
minJacPass.optIterMax = p.itMax;
minJacPass.contrib.push_back(&nodeDistFunc);
- minJacPass.contrib.push_back(&minNCJBarFunc);
+// minJacPass.contrib.push_back(&minNCJBarFunc);
+ minJacPass.contrib.push_back(&minInvCondBarFunc);
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);
- }
+// 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();
+// p.minNCJ = minMaxNCJBarFunc.getMin();
+// p.maxNCJ = minMaxNCJBarFunc.getMax();
+ p.minInvCond = minInvCondBarFunc.getMin();
+ p.maxInvCond = minInvCondBarFunc.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
index 287b3a418559943d25e7099178616392a1877d39..d34d1083d33b01a16617e1d57d5402a836d6c753 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
@@ -33,7 +33,8 @@
class GModel;
struct MeshQualOptParameters {
- double minTargetNCJ, maxTargetNCJ;
+// double minTargetNCJ, maxTargetNCJ;
+ double minTargetInvCond;
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
@@ -51,11 +52,13 @@ struct MeshQualOptParameters {
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 minNCJ, maxNCJ; // after optimization, range of jacobians
+ double minInvCond, maxInvCond; // after optimization, range of jacobians
double CPU; // Time for optimization
MeshQualOptParameters ()
- : minTargetNCJ(0.5), maxTargetNCJ(1.5), weightFixed(1000.),
+// : minTargetNCJ(0.5), maxTargetNCJ(1.5), weightFixed(1000.),
+ : minTargetInvCond(0.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.)