Fixed evaluation and optimization of inverse condition number in 2D

Numeric/CondNumBasis.cpp

@@ -110,8 +110,6 @@ inline void calcGradInvCondNum2D(double dxdX, double dxdY,
                                  const fullMatrix<double> &gSMatY,
                                  fullMatrix<double> &IDI)
 {
-  static const double Eps = 1e-10;
-
   bool posJac = true;
   if (sign) {
     const double lnx = dydX*dzdY - dzdX*dydY, lny = dzdX*dxdY - dxdX*dzdY,            // Local normal from mapping gradients
@@ -138,7 +136,7 @@ inline void calcGradInvCondNum2D(double dxdX, double dxdY,
   const double S12 = Dx*Dy;
   const double S13 = Dx*Dx;
   const double S = 2.*(S2+S5+S12)+ 4.*(S7-S8+S9-S10) + 8.*(S3+S6) + S1+S4+S11+S13;
-  if (S < Eps) {                                                                      // S == 0. -> Ideal element
+  if (S == 0.) {                                                                        // S == 0. -> Ideal element
     for (int j = 0; j<3*numMapNodes; j++) IDI(i, j) = 0.;
     IDI(i, 3*numMapNodes) = posJac ? 1. : -1.;
     return;

contrib/MeshOptimizer/MeshOptPatch.cpp

@@ -626,14 +626,13 @@ void Patch::initInvCondNum()
     }
   }
 
-//  // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
-//  // Jacobians of 3D elements
-//  if ((_dim == 2) && _condNormEl.empty()) {
-//    _condNormEl.resize(nEl());
-////    for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
-//    for (int iEl = 0; iEl < nEl(); iEl++)
-//      calcNormalEl2D(iEl, NS_SQRTNORM, _condNormEl[iEl], true);
-//  }
+  // Set normals to 2D elements
+  if ((_dim == 2) && _condNormEl.empty()) {
+    _condNormEl.resize(nEl());
+//    for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
+    for (int iEl = 0; iEl < nEl(); iEl++)
+      calcNormalEl2D(iEl, NS_UNIT, _condNormEl[iEl], true);
+  }
 }
 
 
@@ -655,12 +654,7 @@ void Patch::invCondNumAndGradients(int iEl, std::vector<double> &condNum,
   }
 
   // Calculate ICN and gradients
-  // TODO: Use signed measure for 2D as well
-  if (_dim == 3) {
-    cnBasis->getSignedInvCondNumAndGradients(nodesXYZ, normals, IDI);
-  }
-  else
-    cnBasis->getInvCondNumAndGradients(nodesXYZ, IDI);
+  cnBasis->getSignedInvCondNumAndGradients(nodesXYZ, _condNormEl[iEl], IDI);
 
   // Inverse condition number
   for (int l = 0; l < numICN; l++) condNum[l] = IDI(l, 3*numMapNodes);

contrib/MeshOptimizer/MeshOptPatch.h

@@ -158,7 +158,7 @@ private:
   std::vector<fullMatrix<double> > _IJacNormEl;                 // Normals to 2D elements for Jacobian regularization and scaling
   std::vector<double> _invIJac;                                 // Initial Jacobians for 3D elements
   std::vector<int> _nICNEl;                                     // Number of inv. cond. number values for an el.
-//  std::vector<fullMatrix<double> > _condNormEl;                 // Normals to 2D elements for inverse conditioning computation
+  std::vector<fullMatrix<double> > _condNormEl;                 // Normals to 2D elements for inverse conditioning computation
 };
 
 

contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp

@@ -25,7 +25,7 @@ struct QualPatchDefParameters : public MeshOptPatchDef
   virtual int inPatch(const SPoint3 &badBary,
                       double limDist, MElement *el) const;
 private:
-  bool _excludeHex, _excludePrism;
+  bool _excludeQuad, _excludeHex, _excludePrism;
   double _idealJacMin, _invCondNumMin;
   double _distanceFactor;
 };
@@ -33,6 +33,7 @@ private:
 
 QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
 {
+  _excludeQuad = p.excludeQuad;
   _excludeHex = p.excludeHex;
   _excludePrism = p.excludePrism;
   _idealJacMin = p.minTargetIdealJac;
@@ -55,6 +56,7 @@ QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
 double QualPatchDefParameters::elBadness(MElement *el) const
 {
   const int typ = el->getType();
+  if (_excludeQuad && (typ == TYPE_QUA)) return 1.;
   if (_excludeHex && (typ == TYPE_HEX)) return 1.;
   if (_excludePrism && (typ == TYPE_PRI)) return 1.;
 //  double jMin, jMax;
@@ -76,7 +78,9 @@ int QualPatchDefParameters::inPatch(const SPoint3 &badBary,
                                     double limDist, MElement *el) const
 {
   const int typ = el->getType();
-  if ((typ == TYPE_HEX) || (typ == TYPE_PRI)) return -1;
+  if (_excludeQuad && (typ == TYPE_QUA)) return -1;
+  if (_excludeHex && (typ == TYPE_HEX)) return -1;
+  if (_excludePrism && (typ == TYPE_PRI)) return -1;
   return testElInDist(badBary, limDist, el) ? 1 : 0;
 }
 

contrib/MeshQualityOptimizer/MeshQualityOptimizer.h

@@ -33,7 +33,7 @@
 class GModel;
 
 struct MeshQualOptParameters {
-  bool excludeHex, excludePrism;
+  bool excludeQuad, excludeHex, excludePrism;
   double minTargetIdealJac;
   double minTargetInvCondNum;
   double weightFixed; // weight of the energy for fixed nodes
@@ -58,8 +58,8 @@ struct MeshQualOptParameters {
   double CPU; // Time for optimization
 
   MeshQualOptParameters ()
-    : excludeHex(false), excludePrism(false), minTargetIdealJac(0.1),
-      minTargetInvCondNum(0.1), weightFixed(1000.),
+    : excludeQuad(false), excludeHex(false), excludePrism(false),
+      minTargetIdealJac(0.1), minTargetInvCondNum(0.1), 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.), CPU(0.),