Mesh optimizer: changed input, improved robustness of barrier functions, added...

Mesh optimizer: changed input, improved robustness of barrier functions, added different scaling for node displacement, removed NCJ contribution, added IdealJac contribution.

CMakeLists.txt

@@ -154,8 +154,10 @@ set(GMSH_API
   contrib/MeshOptimizer/MeshOptPatch.h contrib/MeshOptimizer/MeshOpt.h
   contrib/MeshOptimizer/MeshOptCommon.h contrib/MeshOptimizer/MeshOptimizer.h
   contrib/MeshOptimizer/MeshOptObjContribFunc.h contrib/MeshOptimizer/MeshOptObjContrib.h
+  contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h
   contrib/MeshOptimizer/MeshOptObjectiveFunction.h contrib/MeshOptimizer/MeshOptVertexCoord.h 
-  contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h
+  contrib/MeshQualityOptimizer/MeshQualityObjContribIdealJac.h
+  contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
   contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
   contrib/MathEx/mathex.h)
 

contrib/HighOrderMeshOptimizer/OptHomRun.cpp

@@ -753,12 +753,14 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
 #include "MeshOptimizer.h"
 
 
-struct HOPatchDefParameters : public MeshOptParameters::PatchDefParameters
+struct HOPatchDefParameters : public MeshOptPatchDef
 {
   HOPatchDefParameters(const OptHomParameters &p);
   virtual ~HOPatchDefParameters() {}
-  virtual double elBadness(MElement *el);
-  virtual double maxDistance(MElement *el);
+  virtual double elBadness(MElement *el) const;
+  virtual double maxDistance(MElement *el) const;
+  virtual int inPatch(const SPoint3 &badBary,
+                      double limDist, MElement *el) const;
 private:
   double jacMin, jacMax;
   double distanceFactor;
@@ -790,7 +792,8 @@ HOPatchDefParameters::HOPatchDefParameters(const OptHomParameters &p)
 }
 
 
-double HOPatchDefParameters::elBadness(MElement *el) {
+double HOPatchDefParameters::elBadness(MElement *el) const
+{
   double jmin, jmax;
   el->scaledJacRange(jmin, jmax);
   double badness = std::min(jmin-jacMin, 0.) + std::min(jacMax-jmax, 0.);
@@ -802,11 +805,19 @@ double HOPatchDefParameters::elBadness(MElement *el) {
 }
 
 
-double HOPatchDefParameters::maxDistance(MElement *el) {
+double HOPatchDefParameters::maxDistance(MElement *el) const
+{
   return distanceFactor * el->maxDistToStraight();
 }
 
 
+int HOPatchDefParameters::inPatch(const SPoint3 &badBary,
+                                  double limDist, MElement *el) const
+{
+  return testElInDist(badBary, limDist, el) ? 1 : 0;
+}
+
+
 void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
 {
   Msg::StatusBar(true, "Optimizing high order mesh...");
@@ -820,7 +831,8 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
   par.optDisplay = 30;
   par.verbose = 4;
 
-  ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree);
+  ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree,
+                                                                Patch::LS_MAXNODEDIST);
   ObjContribScaledJac<ObjContribFuncBarrierMovMin> minJacBarFunc(1.);
   minJacBarFunc.setTarget(p.BARRIER_MIN, 1.);
   ObjContribScaledJac<ObjContribFuncBarrierFixMinMovMax> minMaxJacBarFunc(1.);
@@ -828,7 +840,7 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
   ObjContribCADDist<ObjContribFuncSimpleTargetMax> CADDistFunc(p.optCADWeight, p.discrTolerance);
   CADDistFunc.setTarget(p.optCADDistMax);
 
-  MeshOptParameters::PassParameters minJacPass;
+  MeshOptPass minJacPass;
   minJacPass.barrierIterMax = p.optPassMax;
   minJacPass.optIterMax = p.itMax;
   minJacPass.contrib.push_back(&nodeDistFunc);
@@ -837,7 +849,7 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
   par.pass.push_back(minJacPass);
 
   if (p.BARRIER_MAX > 0.) {
-    MeshOptParameters::PassParameters minMaxJacPass;
+    MeshOptPass minMaxJacPass;
     minMaxJacPass.barrierIterMax = p.optPassMax;
     minMaxJacPass.optIterMax = p.itMax;
     minMaxJacPass.contrib.push_back(&nodeDistFunc);

contrib/MeshOptimizer/CMakeLists.txt

@@ -5,6 +5,7 @@
 
 set(SRC
   MeshOpt.cpp
+  MeshOptCommon.cpp
   MeshOptimizer.cpp
   MeshOptPatch.cpp
   MeshOptObjectiveFunction.cpp

contrib/MeshOptimizer/MeshOptCommon.cpp

+// 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
+
+#include "MElement.h"
+#include "MeshOptCommon.h"
+
+
+namespace {
+
+
+// Test intersection between sphere and segment
+bool testSegSphereIntersect(SPoint3 A, SPoint3 B, const SPoint3& P, const double rr)
+{
+  // Test if separating plane between sphere and segment vertices
+  // For each vertex, separation if vertex is outside sphere and P on opposite side
+  // to other seg. vertex w.r.t plane of normal (vertex-P) through vertex
+  const SVector3 PA(P, A), PB(P, B);
+  const double aa = dot(PA, PA), ab = dot(PA, PB);
+  if ((aa > rr) & (ab > aa)) return false;
+  const double  bb = dot(PB, PB);
+  if ((bb > rr) & (ab > bb)) return false;
+
+  // Test if separating plane between sphere and line
+  // For A, separation if projection Q of P on (AB) lies outside the sphere
+  const SVector3 AB(A, B);
+  const double d = ab - aa, e = dot(AB, AB);
+  const  SVector3 PQ = PA * e - d * AB;
+  if (dot(PQ, PQ) > rr * e * e) return false;
+
+  // Return true (intersection) if no separation at all
+  return true;
+}
+
+
+// Test intersection between sphere and triangle
+// Inspired by Christer Ericson, http://realtimecollisiondetection.net/blog/?p=103
+bool testTriSphereIntersect(SPoint3 A, SPoint3 B, SPoint3 C, const SPoint3& P, const double rr)
+{
+  // Test if separating plane between sphere and triangle plane
+  const SVector3 PA(P, A), AB(A, B), AC(A, C);
+  const SVector3 V = crossprod(AB, AC);                               // Normal to triangle plane
+  const double d = dot(PA, V);                                         // Dist. from P to triangle plane times norm of V
+  const double e = dot(V, V);                                         // Norm of V
+  if (d * d > rr * e) return false;                                   // Test if separating plane between sphere and triangle plane
+
+  // Test if separating plane between sphere and triangle vertices
+  const SVector3 PB(P, B), PC(P, B);
+  const double aa = dot(PA, PA), ab = dot(PA, PB), ac = dot(PA, PC);
+  const double bb = dot(PB, PB), bc = dot(PB, PC), cc = dot(PC, PC);
+  if ((aa > rr) & (ab > aa) & (ac > aa)) return false;                    // For each triangle vertex, separation if vertex is outside sphere
+  if ((bb > rr) & (ab > bb) & (bc > bb)) return false;                    // and P on opposite side to other two triangle vertices w.r.t
+  if ((cc > rr) & (ac > cc) & (bc > cc)) return false;                    // plane of normal (vertex-P) through vertex
+
+  // Test if separating plane between sphere and triangle edges
+  const SVector3 BC(B, C);
+  const double d1 = ab - aa, d2 = bc - bb, d3 = ac - cc;
+  const double e1 = dot(AB, AB), e2 = dot(BC, BC), e3 = dot(AC, AC);
+  const  SVector3 PQ1 = PA * e1 - d1 * AB;                                // Q1 projection of P on line (AB)
+  const  SVector3 PQ2 = PB * e2 - d2 * BC;                                // Q2 projection of P on line (BC)
+  const  SVector3 PQ3 = PC * e3 + d3 * AC;                                // Q3 projection of P on line (AC)
+  const  SVector3 PQC = PC * e1 - PQ1;
+  const  SVector3 PQA = PA * e2 - PQ2;
+  const  SVector3 PQB = PB * e3 - PQ3;
+  if ((dot(PQ1, PQ1) > rr * e1 * e1) & (dot(PQ1, PQC) > 0)) return false; // For A, separation if Q lies outside the sphere and if P and C
+  if ((dot(PQ2, PQ2) > rr * e2 * e2) & (dot(PQ2, PQA) > 0)) return false; // are on opposite sides of plane through AB with normal PQ
+  if ((dot(PQ3, PQ3) > rr * e3 * e3) & (dot(PQ3, PQB) > 0)) return false; // Same for other two vertices
+
+  // Return true (intersection) if no separation at all
+  return true;
+}
+
+
+}
+
+
+// Test of intersection element with circle/sphere
+bool MeshOptPatchDef::testElInDist(const SPoint3 &P, double limDist,
+                                   MElement *el) const
+{
+  const double sampleLen = 0.5*limDist;                                   // Distance between sample points
+  const double limDistSq = limDist*limDist;
+
+  if (el->getDim() == 2) {                                                // 2D?
+    for (int iEd = 0; iEd < el->getNumEdges(); iEd++) {                   // Loop over edges of element
+      MEdge ed = el->getEdge(iEd);
+      const SPoint3 A = ed.getVertex(0)->point();
+      const SPoint3 B = ed.getVertex(1)->point();
+      if (testSegSphereIntersect(A, B, P, limDistSq)) {
+        return true;
+      }
+    }
+  }
+  else {                                                                  // 3D
+    for (int iFace = 0; iFace < el->getNumFaces(); iFace++) {             // Loop over faces of element
+      MFace face = el->getFace(iFace);
+      const SPoint3 A = face.getVertex(0)->point();
+      const SPoint3 B = face.getVertex(1)->point();
+      const SPoint3 C = face.getVertex(2)->point();
+      if (face.getNumVertices() == 3)
+        return testTriSphereIntersect(A, B, C, P, limDistSq);
+      else {
+        const SPoint3 D = face.getVertex(3)->point();
+        return (testTriSphereIntersect(A, B, C, P, limDistSq) ||
+                    testTriSphereIntersect(A, C, D, P, limDistSq));
+      }
+    }
+  }
+
+  return false;
+}

contrib/MeshOptimizer/MeshOptCommon.h

@@ -34,17 +34,12 @@
 
 
 class MElement;
+class SPoint3;
 class ObjContrib;
 
 
-struct MeshOptParameters {                             // Parameters controlling the strategy
-  enum { STRAT_CONNECTED, STRAT_ONEBYONE };
-  struct PassParameters {                              // Parameters controlling the optimization procedure in each pass
-    std::vector<ObjContrib*> contrib;                  // Indices of contributions to objective function
-    int optIterMax;                                    // Max. number of opt. iterations each time the barrier is moved
-    int barrierIterMax;                                // Max. number of times the barrier is moved
-  };
-  struct PatchDefParameters {
+class MeshOptPatchDef {
+public:
   int strategy;                                         // Strategy: connected patches or adaptive one-by-one
   int minLayers, maxLayers;                             // Min. and max. nb. of layers around a bad element in patch
   union {
@@ -55,14 +50,32 @@ 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(MElement *el) = 0; // Pointer to function checking the patch distance criterion for a given bad element
+  virtual ~MeshOptPatchDef() {}
+  virtual double elBadness(MElement *el) const = 0;     // Determine "badness" of a given element (for patch creation)
+  virtual double maxDistance(MElement *el) const = 0;   // Compute max. distance to a given bad element for elements in patch
+  virtual int inPatch(const SPoint3 &badBary,           // Determine whether a given element should be included in the patch around a...
+                      double limDist,                   // ... given bad element barycenter, with a limit distance if needed. Output: ...
+                      MElement *el) const = 0;          // ... -1 = excluded, 0 = included only up to minLayers, 1 = included up to maxLayers
+protected:
+  bool testElInDist(const SPoint3 &P, double limDist,   // Test whether an element is within a certain distance from a point
+                    MElement *el) const;
+};
+
+
+struct MeshOptPass {                                    // Parameters controlling the optimization procedure in each pass
+  std::vector<ObjContrib*> contrib;                     // Indices of contributions to objective function
+  int optIterMax;                                       // Max. number of opt. iterations each time the barrier is moved
+  int barrierIterMax;                                   // Max. number of times the barrier is moved
 };
+
+
+struct MeshOptParameters {                              // Parameters controlling the strategy
+  enum { STRAT_CONNECTED, STRAT_ONEBYONE };
   int dim ;                                             // Which dimension to optimize
   bool onlyVisible ;                                    // Apply optimization to visible entities ONLY
   bool fixBndNodes;                                     // If points can move on boundaries
-  PatchDefParameters *patchDef;
-  std::vector<PassParameters> pass;
+  MeshOptPatchDef *patchDef;
+  std::vector<MeshOptPass> pass;
   int optDisplay;                                       // Sampling rate in opt. iterations for display
   int verbose;                                          // Level of information displayed and written to disk
   int success;                                          // Success flag: -1 = fail, 0 = partial fail (target not reached), 1 = success

contrib/MeshOptimizer/MeshOptObjContribFunc.cpp

 // TODO: Copyright
 
 #include <cmath>
+#include "GmshMessage.h"
 #include "MeshOptObjContribFunc.h"
 
 
@@ -31,28 +32,37 @@ bool ObjContribFuncSimpleTargetMax::stagnated(double vMin, double vMax)
 }
 
 
-const double ObjContribFuncBarrier::LOWMARGINMULT = 0.9;
-const double ObjContribFuncBarrier::UPMARGINMULT = 1.1;
+const double ObjContribFuncBarrier::MARGINCOEFF = 0.1;
 const double ObjContribFuncBarrier::STAGTHRESHOLD = 0.01;
 
 
 ObjContribFuncBarrier::ObjContribFuncBarrier() :
-    _target(0.), _barrier(0.), _init(0.), _opt(0.)
+    _target(0.), _barrier(0.), _init(0.), _opt(0.), _defaultMargin(0.)
 {
 }
 
 
-void ObjContribFuncBarrier::setTarget(double target, double opt)
+void ObjContribFuncBarrier::setTarget(double target, double opt, double defaultMargin)
 {
   _target = target;
   _opt = opt;
+  _defaultMargin = defaultMargin;
+  if (_defaultMargin == 0.) {
+    if (opt != 0.)
+      _defaultMargin = MARGINCOEFF*fabs(opt);
+    else if (target != 0.)
+      _defaultMargin = MARGINCOEFF*fabs(target);
+    else
+      Msg::Warning("Could not find value to define a scale for default barrier margin");
+  }
 }
 
 
 void ObjContribFuncBarrierMovMin::updateParameters(double vMin, double vMax)
 {
   _init = vMin;
-  _barrier = vMin > 0. ? LOWMARGINMULT*vMin : UPMARGINMULT*vMin;
+  const double margin = (vMin == 0.) ? _defaultMargin : MARGINCOEFF*fabs(vMin);
+  _barrier = vMin - margin;
 }
 
 
@@ -65,7 +75,8 @@ bool ObjContribFuncBarrierMovMin::stagnated(double vMin, double vMax)
 void ObjContribFuncBarrierMovMax::updateParameters(double vMin, double vMax)
 {
   _init = vMax;
-  _barrier = vMax > 0. ? UPMARGINMULT*vMax : LOWMARGINMULT*vMax;
+  const double margin = (vMax == 0.) ? _defaultMargin : MARGINCOEFF*fabs(vMax);
+  _barrier = vMax + margin;
 }
 
 
@@ -78,5 +89,6 @@ bool ObjContribFuncBarrierMovMax::stagnated(double vMin, double vMax)
 void ObjContribFuncBarrierFixMinMovMax::initialize(double vMin, double vMax)
 {
   ObjContribFuncBarrierMovMax::initialize(vMin, vMax);
-  _fixedMinBarrier = vMin > 0. ? LOWMARGINMULT*vMin : UPMARGINMULT*vMin;
+  const double margin = (vMin == 0.) ? _defaultMargin : MARGINCOEFF*fabs(vMin);
+  _fixedMinBarrier = vMin - margin;
 }

contrib/MeshOptimizer/MeshOptObjContribFunc.h

@@ -38,13 +38,16 @@ class ObjContribFuncBarrier
 {
 public:
   ObjContribFuncBarrier();
-  void setTarget(double target, double opt=1.);
+  void setTarget(double target, double opt=1.,
+                 double defaultMargin=0.);
 
 protected:
-  static const double LOWMARGINMULT, UPMARGINMULT;                                // Upper and lower margins w.r.t. min. & max. to set barrier
+  static const double MARGINCOEFF;                                                // Margin coefficient w.r.t. min. & max. to set barrier
   static const double STAGTHRESHOLD;                                              // Threshold to consider that measures stagnates
   double _opt;                                                                    // Optimal value of measure in barrier function
-  double _barrier, _target, _init;                                                // Current barrier, target and initial values of min./max. of measure
+  double _defaultMargin;                                                          // Default margin value to set barrier w.r.t. of min./max. of measure
+  double _barrier, _target;                                                       // Current barrier and target of min./max. of measure
+  double _init;                                                                   // Initial value of min./max. of measure
   static double logBarrier(double v, double barrier, double opt);
   static double diffLogBarrier(double v, double barrier, double opt);
 };

contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h

@@ -11,7 +11,8 @@ template<class FuncType>
 class ObjContribScaledNodeDispSq : public ObjContrib, public FuncType
 {
 public:
-  ObjContribScaledNodeDispSq(double weightFixed, double weightFree);
+  ObjContribScaledNodeDispSq(double weightFixed, double weightFree,
+                             Patch::LengthScaling scaling);
   virtual ~ObjContribScaledNodeDispSq() {}
   virtual ObjContrib *copy() const;
   virtual void initialize(Patch *mesh);
@@ -25,15 +26,16 @@ public:
 protected:
   Patch *_mesh;
   double _weightFixed, _weightFree;
-
+  Patch::LengthScaling _scaling;
 };
 
 
 template<class FuncType>
 ObjContribScaledNodeDispSq<FuncType>::ObjContribScaledNodeDispSq(double weightFixed,
-                                                                    double weightFree) :
+                                                                 double weightFree,
+                                                                 Patch::LengthScaling scaling) :
   ObjContrib("ScaledNodeDispSq", FuncType::getNamePrefix()+"ScaledNodeDispSq"),
-  _mesh(0), _weightFixed(weightFixed), _weightFree(weightFree)
+  _mesh(0), _weightFixed(weightFixed), _weightFree(weightFree), _scaling(scaling)
 {
 }
 
@@ -49,7 +51,7 @@ template<class FuncType>
 void ObjContribScaledNodeDispSq<FuncType>::initialize(Patch *mesh)
 {
   _mesh = mesh;
-  _mesh->initScaledNodeDispSq();
+  _mesh->initScaledNodeDispSq(_scaling);
   updateMinMax();
   FuncType::initialize(_min, _max);
 }

contrib/MeshOptimizer/MeshOptPatch.cpp

@@ -41,7 +41,8 @@
 
 Patch::Patch(const std::map<MElement*,GEntity*> &element2entity,
            const std::set<MElement*> &els, std::set<MVertex*> &toFix,
-           bool fixBndNodes)
+           bool fixBndNodes) :
+  _typeLengthScale(LS_NONE), _invLengthScaleSq(0.)
 {
 
   _dim = (*els.begin())->getDim();
@@ -212,24 +213,36 @@ void Patch::writeMSH(const char *filename)
   fclose(f);
 }
 
-// TODO: allow user to choose type of scaling length
-void Patch::initScaledNodeDispSq()
+
+void Patch::initScaledNodeDispSq(LengthScaling scaling)
 {
-  if (_invLengthScaleSq == 0.) {
+  if ((_invLengthScaleSq == 0.) || _typeLengthScale != scaling) {
+    _typeLengthScale = scaling;
     double maxDSq = 0.;
+    switch(scaling) {
+      case LS_MAXNODEDIST : {
         for (int iEl = 0; iEl < nEl(); iEl++) {
           const double d = el(iEl)->maxDistToStraight(), dd = d*d;
           if (dd > maxDSq) maxDSq = dd;
         }
-    if (maxDSq < 1.e-10) {
-      double maxSSq = 0.;
+        break;
+      }
+      case LS_MAXOUTERRADIUS : {
+        for (int iEl = 0; iEl < nEl(); iEl++) {
+          const double d = el(iEl)->getOuterRadius(), dd = d*d;
+          if (dd > maxDSq) maxDSq = dd;
+        }
+        break;
+      }
+      case LS_MINEDGELENGTH : {
         for (int iEl = 0; iEl < nEl(); iEl++) {
-        const double s = el(iEl)->getOuterRadius(), ss = s*s;
-        if (ss > maxSSq) maxSSq = ss;
+          const double d = el(iEl)->minEdge(), dd = d*d;
+          if (dd > maxDSq) maxDSq = dd;
+        }
+        break;
       }
-      _invLengthScaleSq = 1./maxSSq;
     }
-    else _invLengthScaleSq = 1./maxDSq;
+    _invLengthScaleSq = 1./maxDSq;
   }
 }
 
@@ -266,10 +279,11 @@ void Patch::initScaledJac()
 
   // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
   // Jacobians of 3D elements
-  if ((_dim == 2) && _scaledNormEl.empty()) {
-    _scaledNormEl.resize(nEl());
+  if ((_dim == 2) && _JacNormEl.empty()) {
+    _JacNormEl.resize(nEl());
 //    for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
-    for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, NS_INVNORM);
+    for (int iEl = 0; iEl < nEl(); iEl++)
+      calcNormalEl2D(iEl, NS_INVNORM, _JacNormEl[iEl], false);
   }
   else if (_invStraightJac.empty()) {
     _invStraightJac.resize(nEl(), 1.);
@@ -293,7 +307,8 @@ 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, NORMALSCALING scaling)
+void Patch::calcNormalEl2D(int iEl, NormalScaling scaling,
+                           fullMatrix<double> &elNorm, bool ideal)
 {
   const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
 
@@ -319,22 +334,26 @@ void Patch::calcNormalEl2D(int iEl, NORMALSCALING scaling)
 //    geoNorm = ((GFace*)ge)->normal(param);
 //  }
 
-  fullMatrix<double> uElNorm(1, 3);
-  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 (scaling != 0) {
-    double factor = (scaling == NS_SQRTNORM) ? 1./norm : sqrt(norm);
+  const double norm = jac->getPrimNormal2D(primNodesXYZ, elNorm, ideal);
+  double factor;
+  switch (scaling) {
+    case NS_UNIT:
+      factor = 1.;
+      break;
+    case NS_INVNORM:
+      factor = 1./norm;
+      break;
+    case NS_SQRTNORM:
+      factor = sqrt(norm);
+      break;
+  }
 //  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;
 //  }
   elNorm.scale(factor);   // Re-scaling normal here is faster than an extra scaling operation on the Jacobian
 }
-  else
-    _unitNormEl[iEl] = SVector3(uElNorm(0, 0), uElNorm(0, 1), uElNorm(0, 2));
-}
 
 
 void Patch::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
@@ -356,7 +375,7 @@ void Patch::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
 
   // Calculate Jacobian and gradients, scale if 3D (already scaled by
   // regularization normals in 2D)
-  jacBasis->getSignedJacAndGradients(nodesXYZ,_scaledNormEl[iEl],JDJ);
+  jacBasis->getSignedJacAndGradients(nodesXYZ,_JacNormEl[iEl],JDJ);
   if (_dim == 3) JDJ.scale(_invStraightJac[iEl]);
 
   // Transform Jacobian and gradients from Lagrangian to Bezier basis
@@ -504,89 +523,90 @@ bool Patch::bndDistAndGradients(int iEl, double &f, std::vector<double> &gradF,
 }
 
 
-void Patch::initNCJ()
+void Patch::initIdealJac()
 {
   // Initialize _nBezEl
-  if (_nNCJEl.empty()) {
-    _nNCJEl.resize(nEl());
+  if (_nIJacEl.empty()) {
+    _nIJacEl.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;
-      }
+      _nIJacEl[iEl] = _el[iEl]->getJacobianFuncSpace()->getNumJacNodes();
   }
 
   // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
   // Jacobians of 3D elements
-  if ((_dim == 2) && _unitNormEl.empty()) {
-    _unitNormEl.resize(nEl());
+  if ((_dim == 2) && _IJacNormEl.empty()) {
+    _IJacNormEl.resize(nEl());
 //    for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
-    for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, NS_UNIT);
-  }
-}
-
-
-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];
-  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], _xyz[iV1], _xyz[iV2], _unitNormEl[iEl], NCJ);
-    break;
+    for (int iEl = 0; iEl < nEl(); iEl++)
+      calcNormalEl2D(iEl, NS_INVNORM, _IJacNormEl[iEl], true);
   }
-  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], _xyz[iV1], _xyz[iV2], _xyz[iV3], _unitNormEl[iEl], NCJ);
-    break;
+  else if (_invStraightJac.empty()) {
+    _invIJac.resize(nEl(), 1.);
+    for (int iEl = 0; iEl < nEl(); iEl++) {
+      int nEd = _el[iEl]->getNumEdges();
+      double sumEdLength = 0.;
+      for(int iEd = 0; iEd < nEd; iEd++)
+        sumEdLength += _el[iEl]->getEdge(iEd).length();
+      const double invMeanEdLength = double(nEd)/sumEdLength;
+      _invIJac[iEl] = invMeanEdLength*invMeanEdLength*invMeanEdLength;
     }
   }
 }
 
 
-void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<double> &gNCJ)
+void Patch::idealJacAndGradients(int iEl, std::vector<double> &iJ, std::vector<double> &gIJ)
 {
-  const int &numNCJVal = _nNCJEl[iEl];
-  const int &numNodes = _nNodEl[iEl];
-
-  std::vector<int> &iVEl = _el2V[iEl];
-  std::vector<double> dNCJ(numNCJVal*3*numNodes);
-  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], _xyz[iV1], _xyz[iV2],
-                                _unitNormEl[iEl], NCJ, dNCJ);
-    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], _xyz[iV1], _xyz[iV2], _xyz[iV3],
-                                  _unitNormEl[iEl], NCJ, dNCJ);
-    break;
-  }
+  const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
+  const int &numJacNodes = _nIJacEl[iEl];
+  const int &numMapNodes = _nNodEl[iEl];
+  fullMatrix<double> JDJ(numJacNodes,3*numMapNodes+1), BDB(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();
   }
 
-  // Gradients of the NCJ
+  // Calculate Jacobian and gradients, scale if 3D (already scaled by
+  // regularization normals in 2D)
+  jacBasis->getSignedIdealJacAndGradients(nodesXYZ,_IJacNormEl[iEl],JDJ);
+  if (_dim == 3) JDJ.scale(_invIJac[iEl]);
+//  if (_el[iEl]->getNum() == 90370) {
+//    std::cout << "DBGTT: bad el.: " << _el[iEl]->getNum() << "\n";
+//    for (int i = 0; i < numMapNodes; i++)
+//      std::cout << "DBGTT: {x,y,z}" << i << " = (" << nodesXYZ(i,0)
+//                << ", " << nodesXYZ(i,1) << ", " << nodesXYZ(i,2) << ")\n";
+//   for (int l = 0; l < numJacNodes; l++) {
+//      for (int i = 0; i < numMapNodes; i++)
+//        std::cout << "DBGTT: dJ" << l << "d{x,y,z}" << i << " = (" << JDJ(l, i)
+//                  << ", " << JDJ(l, i+numMapNodes) << ", " << JDJ(l, i+2*numMapNodes)<< ")\n";
+//      std::cout << "DBGTT: J" << l << " = " << JDJ(l, 3*numMapNodes)<< "\n";
+//    }
+//  }
+
+  // Transform Jacobian and gradients from Lagrangian to Bezier basis
+  jacBasis->lag2Bez(JDJ,BDB);
+
+  // Scaled jacobian
+  for (int l = 0; l < numJacNodes; l++) iJ [l] = BDB (l,3*numMapNodes);
+
+  // Gradients of the scaled jacobian
   int iPC = 0;
-  std::vector<SPoint3> gXyzV(numNCJVal);
-  std::vector<SPoint3> gUvwV(numNCJVal);
-  for (int i = 0; i < numNodes; i++) {
+  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 < numNCJVal; l++) {
-        int ind = (l*numNodes+i)*3;
-        gXyzV[l] = SPoint3(dNCJ[ind], dNCJ[ind+1], dNCJ[ind+2]);
-      }
+      for (int l = 0; l < numJacNodes; l++)
+        gXyzV [l] = SPoint3(BDB(l,i), BDB(l,i+numMapNodes), BDB(l,i+2*numMapNodes));
       _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];
+      for (int l = 0; l < numJacNodes; l++) {
+        gIJ[indGIJac(iEl,l,iPC)] = gUvwV[l][0];
+        if (_nPCFV[iFVi] >= 2) gIJ[indGIJac(iEl,l,iPC+1)] = gUvwV[l][1];
+        if (_nPCFV[iFVi] == 3) gIJ[indGIJac(iEl,l,iPC+2)] = gUvwV[l][2];
       }
       iPC += _nPCFV[iFVi];
     }
@@ -594,7 +614,7 @@ void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<doubl
 }
 
 
-void Patch::initInvCond()
+void Patch::initCondNum()
 {
   // Initialize _nBezEl
   if (_nBezEl.empty()) {
@@ -605,16 +625,17 @@ void Patch::initInvCond()
 
   // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
   // Jacobians of 3D elements
-  if ((_dim == 2) && _metricNormEl.empty()) {
-    _metricNormEl.resize(nEl());
+  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);
+    for (int iEl = 0; iEl < nEl(); iEl++)
+      calcNormalEl2D(iEl, NS_SQRTNORM, _condNormEl[iEl], true);
   }
 }
 
 
-void Patch::invCondAndGradients(int iEl, std::vector<double> &invCond,
-                                std::vector<double> &gInvCond)
+void Patch::condNumAndGradients(int iEl, std::vector<double> &condNum,
+                                std::vector<double> &gCondNum)
 {
   const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
   const int &numJacNodes = _nBezEl[iEl];
@@ -632,10 +653,10 @@ void Patch::invCondAndGradients(int iEl, std::vector<double> &invCond,
 
   // Calculate Jacobian and gradients, scale if 3D (already scaled by
   // regularization normals in 2D)
-  jacBasis->getInvCondAndGradients(nodesXYZ, _metricNormEl[iEl], IDI);
+  jacBasis->getCondNumAndGradients(nodesXYZ, _condNormEl[iEl], IDI);
 
   // Inverse condition number
-  for (int l = 0; l < numJacNodes; l++) invCond[l] = IDI(l, 3*numMapNodes);
+  for (int l = 0; l < numJacNodes; l++) condNum[l] = IDI(l, 3*numMapNodes);
 
   // Gradients of the inverse condition number
   int iPC = 0;
@@ -649,9 +670,9 @@ void Patch::invCondAndGradients(int iEl, std::vector<double> &invCond,
                             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];
+        gCondNum[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
+        if (_nPCFV[iFVi] >= 2) gCondNum[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+        if (_nPCFV[iFVi] == 3) gCondNum[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
       }
       iPC += _nPCFV[iFVi];
     }

contrib/MeshOptimizer/MeshOptPatch.h

@@ -78,7 +78,8 @@ public:
   void writeMSH(const char *filename);
 
   // Node distance measure
-  void initScaledNodeDispSq();
+  enum LengthScaling { LS_NONE, LS_MAXNODEDIST, LS_MAXOUTERRADIUS, LS_MINEDGELENGTH };
+  void initScaledNodeDispSq(LengthScaling scaling);
   inline double invLengthScaleSq() { return _invLengthScaleSq; }
   double scaledNodeDispSq(int iFV);
   void gradScaledNodeDispSq(int iFV, std::vector<double> &gDSq);
@@ -93,18 +94,15 @@ public:
   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);
-  void initInvCond();
-  void invCondAndGradients(int iEl, std::vector<double> &invCond, std::vector<double> &ginvCond);
+  inline const int &nIJacEl(int iEl) { return _nIJacEl[iEl]; }
+  inline int indGIJac(int iEl, int l, int iPC) { return iPC*_nIJacEl[iEl]+l; }
+  void initIdealJac();
+  void idealJacAndGradients(int iEl, std::vector<double> &NCJ, std::vector<double> &gNCJ);
+  void initCondNum();
+  void condNumAndGradients(int iEl, std::vector<double> &condNum, std::vector<double> &gCondNum);
 
 private:
 
-  enum NORMALSCALING { NS_UNIT, NS_INVNORM, NS_SQRTNORM };
-
   // Mesh entities and variables
   int _dim;
   int _nPC;                                         // Total nb. of parametric coordinates
@@ -141,19 +139,23 @@ private:
   }
 
   // Node displacement
+  LengthScaling _typeLengthScale;
   double _invLengthScaleSq;                                     // Square inverse of a length for node displacement scaling
 
   // High-order: scaled Jacobian and metric measures
+  enum NormalScaling { NS_UNIT, NS_INVNORM, NS_SQRTNORM };
   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<fullMatrix<double> > _JacNormEl;                  // 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, NORMALSCALING scaling);
+  void calcNormalEl2D(int iEl, NormalScaling scaling,
+                      fullMatrix<double> &elNorm, bool ideal);
 
   // 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
+  std::vector<int> _nIJacEl;                                    // Number of NCJ values for an el
+  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<fullMatrix<double> > _condNormEl;                 // Normals to 2D elements for inverse conditioning computation
 };
 
 

contrib/MeshOptimizer/MeshOptimizer.cpp

@@ -71,109 +71,10 @@ static std::set<MVertex *> getAllBndVertices
 }
 
 
-// Test intersection between sphere and segment
-static bool testSegSphereIntersect(SPoint3 A, SPoint3 B, const SPoint3& P, const double rr)
-{
-  // Test if separating plane between sphere and segment vertices
-  // For each vertex, separation if vertex is outside sphere and P on opposite side
-  // to other seg. vertex w.r.t plane of normal (vertex-P) through vertex
-  const SVector3 PA(P, A), PB(P, B);
-  const double aa = dot(PA, PA), ab = dot(PA, PB);
-  if ((aa > rr) & (ab > aa)) return false;
-  const double  bb = dot(PB, PB);
-  if ((bb > rr) & (ab > bb)) return false;
-
-  // Test if separating plane between sphere and line
-  // For A, separation if projection Q of P on (AB) lies outside the sphere
-  const SVector3 AB(A, B);
-  const double d = ab - aa, e = dot(AB, AB);
-  const  SVector3 PQ = PA * e - d * AB;
-  if (dot(PQ, PQ) > rr * e * e) return false;
-
-  // Return true (intersection) if no separation at all
-  return true;
-}
-
-
-// Test intersection between sphere and triangle
-// Inspired by Christer Ericson, http://realtimecollisiondetection.net/blog/?p=103
-static bool testTriSphereIntersect(SPoint3 A, SPoint3 B, SPoint3 C, const SPoint3& P, const double rr)
-{
-  // Test if separating plane between sphere and triangle plane
-  const SVector3 PA(P, A), AB(A, B), AC(A, C);
-  const SVector3 V = crossprod(AB, AC);                               // Normal to triangle plane
-  const double d = dot(PA, V);                                         // Dist. from P to triangle plane times norm of V
-  const double e = dot(V, V);                                         // Norm of V
-  if (d * d > rr * e) return false;                                   // Test if separating plane between sphere and triangle plane
-
-  // Test if separating plane between sphere and triangle vertices
-  const SVector3 PB(P, B), PC(P, B);
-  const double aa = dot(PA, PA), ab = dot(PA, PB), ac = dot(PA, PC);
-  const double bb = dot(PB, PB), bc = dot(PB, PC), cc = dot(PC, PC);
-  if ((aa > rr) & (ab > aa) & (ac > aa)) return false;                    // For each triangle vertex, separation if vertex is outside sphere
-  if ((bb > rr) & (ab > bb) & (bc > bb)) return false;                    // and P on opposite side to other two triangle vertices w.r.t
-  if ((cc > rr) & (ac > cc) & (bc > cc)) return false;                    // plane of normal (vertex-P) through vertex
-
-  // Test if separating plane between sphere and triangle edges
-  const SVector3 BC(B, C);
-  const double d1 = ab - aa, d2 = bc - bb, d3 = ac - cc;
-  const double e1 = dot(AB, AB), e2 = dot(BC, BC), e3 = dot(AC, AC);
-  const  SVector3 PQ1 = PA * e1 - d1 * AB;                                // Q1 projection of P on line (AB)
-  const  SVector3 PQ2 = PB * e2 - d2 * BC;                                // Q2 projection of P on line (BC)
-  const  SVector3 PQ3 = PC * e3 + d3 * AC;                                // Q3 projection of P on line (AC)
-  const  SVector3 PQC = PC * e1 - PQ1;
-  const  SVector3 PQA = PA * e2 - PQ2;
-  const  SVector3 PQB = PB * e3 - PQ3;
-  if ((dot(PQ1, PQ1) > rr * e1 * e1) & (dot(PQ1, PQC) > 0)) return false; // For A, separation if Q lies outside the sphere and if P and C
-  if ((dot(PQ2, PQ2) > rr * e2 * e2) & (dot(PQ2, PQA) > 0)) return false; // are on opposite sides of plane through AB with normal PQ
-  if ((dot(PQ3, PQ3) > rr * e3 * e3) & (dot(PQ3, PQB) > 0)) return false; // Same for other two vertices
-
-  // Return true (intersection) if no separation at all
-  return true;
-}
-
-
-// Approximate test of intersection element with circle/sphere by sampling
-static bool testElInDist(const SPoint3 P, double limDist, MElement *el)
-{
-  const double sampleLen = 0.5*limDist;                                   // Distance between sample points
-  const double limDistSq = limDist*limDist;
-
-  if (el->getDim() == 2) {                                                // 2D?
-    for (int iEd = 0; iEd < el->getNumEdges(); iEd++) {                   // Loop over edges of element
-      MEdge ed = el->getEdge(iEd);
-      const SPoint3 A = ed.getVertex(0)->point();
-      const SPoint3 B = ed.getVertex(1)->point();
-      if (testSegSphereIntersect(A, B, P, limDistSq)) {
-        return true;
-      }
-    }
-  }
-  else {                                                                  // 3D
-    for (int iFace = 0; iFace < el->getNumFaces(); iFace++) {             // Loop over faces of element
-      MFace face = el->getFace(iFace);
-      const SPoint3 A = face.getVertex(0)->point();
-      const SPoint3 B = face.getVertex(1)->point();
-      const SPoint3 C = face.getVertex(2)->point();
-      if (face.getNumVertices() == 3)
-        return testTriSphereIntersect(A, B, C, P, limDistSq);
-      else {
-        const SPoint3 D = face.getVertex(3)->point();
-        return (testTriSphereIntersect(A, B, C, P, limDistSq) ||
-                    testTriSphereIntersect(A, C, D, P, limDistSq));
-      }
-    }
-  }
-
-  return false;
-}
-
-
-static std::set<MElement*> getSurroundingPatch(MElement *el, int minLayers,
-                            int maxLayers, double limDist,
+static std::set<MElement*> getSurroundingPatch(MElement *el, const MeshOptPatchDef *patchDef,
+                            double limDist, int maxLayers,
                             const std::map<MVertex*, std::vector<MElement*> > &vertex2elements)
 {
-
   const SPoint3 pnt = el->barycenter(true);
 
   std::set<MElement*> patch;
@@ -190,12 +91,14 @@ static std::set<MElement*> getSurroundingPatch(MElement *el, int minLayers,
           ((*it)->getVertex(i))->second;
         for (std::vector<MElement*>::const_iterator itN = neighbours.begin();
              itN != neighbours.end(); ++itN) {
-          if ((patch.find(*itN) == patch.end()) &&
-              ((d < minLayers) || testElInDist(pnt, limDist, *itN)))
+          if (patch.find(*itN) == patch.end()) {
+            const int elIn = patchDef->inPatch(pnt, limDist, *itN);
+            if ((elIn > 0) || ((d < patchDef->minLayers) && (elIn == 0)))
               if (patch.insert(*itN).second) currentLayer.push_back(*itN);      // Assume that if an el is too far, its neighbours are too far as well
           }
         }
       }
+    }
     lastLayer = currentLayer;
   }
 
@@ -251,8 +154,8 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
   primPatches.reserve(badElements.size());
   for (std::set<MElement*>::const_iterator it = badElements.begin(); it != badElements.end(); ++it) {
     const double limDist = par.patchDef->maxDistance(*it);
-    primPatches.push_back(getSurroundingPatch(*it, par.patchDef->minLayers,
-                            par.patchDef->maxLayers, limDist, vertex2elements));
+    primPatches.push_back(getSurroundingPatch(*it, par.patchDef, limDist,
+                                par.patchDef->maxLayers, vertex2elements));
   }
 
   // Compute patch connectivity
@@ -404,8 +307,8 @@ static void optimizeOneByOne
 
       // Set up patch
       const double limDist = par.patchDef->maxDistance(worstEl);
-      std::set<MElement*> toOptimizePrim = getSurroundingPatch(worstEl, par.patchDef->minLayers,
-                                                          maxLayers, limDist, vertex2elements);
+      std::set<MElement*> toOptimizePrim = getSurroundingPatch(worstEl, par.patchDef, limDist,
+                                                                    maxLayers, vertex2elements);
       std::set<MVertex*> toFix = getAllBndVertices(toOptimizePrim, vertex2elements);
       std::set<MElement*> toOptimize;
       std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),

contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h → contrib/MeshQualityOptimizer/MeshQualityObjContribIdealJac.h

 // TODO: Copyright
 
-#ifndef _MESHQUALITYCONTRIBNCJ_H_
-#define _MESHQUALITYCONTRIBNCJ_H_
+#ifndef _MESHQUALITYOBJCONTRIBIDEALJAC_H_
+#define _MESHQUALITYOBJCONTRIBIDEALJAC_H_
 
 #include "MeshOptPatch.h"
 #include "MeshOptObjContrib.h"
 
 
 template<class FuncType>
-class ObjContribNCJ : public ObjContrib, public FuncType
+class ObjContribIdealJac : public ObjContrib, public FuncType
 {
 public:
-  ObjContribNCJ(double weight);
-  virtual ~ObjContribNCJ() {}
+  ObjContribIdealJac(double weight);
+  virtual ~ObjContribIdealJac() {}
   virtual ObjContrib *copy() const;
   virtual void initialize(Patch *mesh);
   virtual bool fail() { return _min <= 0.; }
@@ -29,47 +29,47 @@ protected:
 
 
 template<class FuncType>
-ObjContribNCJ<FuncType>::ObjContribNCJ(double weight) :
-  ObjContrib("NCJ", FuncType::getNamePrefix()+"NCJ"),
+ObjContribIdealJac<FuncType>::ObjContribIdealJac(double weight) :
+  ObjContrib("IdealJac", FuncType::getNamePrefix()+"IdealJac"),
   _mesh(0), _weight(weight)
 {
 }
 
 
 template<class FuncType>
-ObjContrib *ObjContribNCJ<FuncType>::copy() const
+ObjContrib *ObjContribIdealJac<FuncType>::copy() const
 {
-  return new ObjContribNCJ<FuncType>(*this);
+  return new ObjContribIdealJac<FuncType>(*this);
 }
 
 
 template<class FuncType>
-void ObjContribNCJ<FuncType>::initialize(Patch *mesh)
+void ObjContribIdealJac<FuncType>::initialize(Patch *mesh)
 {
   _mesh = mesh;
-  _mesh->initNCJ();
+  _mesh->initIdealJac();
   updateMinMax();
   FuncType::initialize(_min, _max);
 }
 
 
 template<class FuncType>
-bool ObjContribNCJ<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
+bool ObjContribIdealJac<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]);
+    std::vector<double> iJ(_mesh->nIJacEl(iEl));                                       // Scaled Jacobians
+    std::vector<double> gIJ(_mesh->nIJacEl(iEl)*_mesh->nPCEl(iEl));                    // Gradients of scaled Jacobians
+    _mesh->idealJacAndGradients(iEl, iJ, gIJ);
+    for (int l = 0; l < _mesh->nIJacEl(iEl); l++) {                                    // Add contribution for each Bezier coeff.
+      Obj += _weight * FuncType::compute(iJ[l]);
+      const double dfact = _weight * FuncType::computeDiff(iJ[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]);
+        gradObj[_mesh->indPCEl(iEl, iPC)] += dfact * gIJ[_mesh->indGIJac(iEl, l, iPC)];
+      _min = std::min(_min, iJ[l]);
+      _max = std::max(_max, iJ[l]);
     }
   }
 
@@ -78,20 +78,21 @@ bool ObjContribNCJ<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gra
 
 
 template<class FuncType>
-void ObjContribNCJ<FuncType>::updateMinMax()
+void ObjContribIdealJac<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->nNCJEl(iEl); l++) {                      // Check each Bezier coeff.
-      _min = std::min(_min, NCJ[l]);
-      _max = std::max(_max, NCJ[l]);
+    std::vector<double> iJ(_mesh->nIJacEl(iEl));                         // Scaled Jacobians
+    std::vector<double> dumGIJ(_mesh->nIJacEl(iEl)*_mesh->nPCEl(iEl));   // Gradients of scaled Jacobians
+    _mesh->idealJacAndGradients(iEl, iJ, dumGIJ);
+    for (int l = 0; l < _mesh->nIJacEl(iEl); l++) {                      // Check each Bezier coeff.
+      _min = std::min(_min, iJ[l]);
+      _max = std::max(_max, iJ[l]);
     }
   }
 }
 
 
-#endif /* _MESHQUALITYCONTRIBNCJ_H_ */
+#endif /* _MESHQUALITYOBJCONTRIBIDEALJAC_H_ */

contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h

@@ -47,7 +47,7 @@ template<class FuncType>
 void ObjContribInvCond<FuncType>::initialize(Patch *mesh)
 {
   _mesh = mesh;
-  _mesh->initInvCond();
+  _mesh->initCondNum();
   updateMinMax();
   FuncType::initialize(_min, _max);
 }
@@ -62,7 +62,7 @@ bool ObjContribInvCond<FuncType>::addContrib(double &Obj, alglib::real_1d_array
   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);
+    _mesh->condNumAndGradients(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]);
@@ -86,7 +86,7 @@ void ObjContribInvCond<FuncType>::updateMinMax()
   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);
+    _mesh->condNumAndGradients(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]);

contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp

@@ -10,35 +10,37 @@
 #include "MeshOptObjContribFunc.h"
 #include "MeshOptObjContrib.h"
 #include "MeshOptObjContribScaledNodeDispSq.h"
-#include "MeshQualityObjContribNCJ.h"
+#include "MeshQualityObjContribIdealJac.h"
 #include "MeshQualityObjContribInvCond.h"
 #include "MeshOptimizer.h"
 #include "MeshQualityOptimizer.h"
 
 
-struct QualPatchDefParameters : public MeshOptParameters::PatchDefParameters
+struct QualPatchDefParameters : public MeshOptPatchDef
 {
   QualPatchDefParameters(const MeshQualOptParameters &p);
   virtual ~QualPatchDefParameters() {}
-  virtual double elBadness(MElement *el);
-  virtual double maxDistance(MElement *el);
+  virtual double elBadness(MElement *el) const;
+  virtual double maxDistance(MElement *el) const;
+  virtual int inPatch(const SPoint3 &badBary,
+                      double limDist, MElement *el) const;
 private:
-//  double NCJMin, NCJMax;
-  double invCondMin;
-  double distanceFactor;
+  bool _excludeHex, _excludePrism;
+  double _idealJacMin;
+  double _distanceFactor;
 };
 
 
 QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
 {
-//  NCJMin = p.minTargetNCJ;
-//  NCJMax = (p.maxTargetNCJ > 0.) ? p.maxTargetNCJ : 1.e300;
-  invCondMin = p.minTargetInvCond;
+  _excludeHex = p.excludeHex;
+  _excludePrism = p.excludePrism;
+  _idealJacMin = p.minTargetIdealJac;
   strategy = (p.strategy == 1) ? MeshOptParameters::STRAT_ONEBYONE :
                                         MeshOptParameters::STRAT_CONNECTED;
   minLayers = (p.dim == 3) ? 1 : 0;
   maxLayers = p.nbLayers;
-  distanceFactor = p.distanceFactor;
+  _distanceFactor = p.distanceFactor;
   if (strategy == MeshOptParameters::STRAT_CONNECTED)
     weakMerge = (p.strategy == 2);
   else {
@@ -49,32 +51,29 @@ 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;
-//  }
-//  }
-  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;
+double QualPatchDefParameters::elBadness(MElement *el) const
+{
+  const int typ = el->getType();
+  if (_excludeHex && (typ == TYPE_HEX)) return 1.;
+  if (_excludePrism && (typ == TYPE_PRI)) return 1.;
+  double jMin, jMax;
+  el->idealJacRange(jMin, jMax);
+  return jMin-_idealJacMin;
+}
+
+
+double QualPatchDefParameters::maxDistance(MElement *el) const
+{
+  return _distanceFactor * el->maxEdge();
 }
 
 
-double QualPatchDefParameters::maxDistance(MElement *el) {
-  return distanceFactor * el->getOuterRadius();
+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;
+  return testElInDist(badBary, limDist, el) ? 1 : 0;
 }
 
 
@@ -88,42 +87,24 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
   par.fixBndNodes = p.fixBndNodes;
   QualPatchDefParameters patchDef(p);
   par.patchDef = &patchDef;
-  par.optDisplay = 30;
+  par.optDisplay = 20;
   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.);
-  ObjContribInvCond<ObjContribFuncBarrierMovMin> minInvCondBarFunc(1.);
-  minInvCondBarFunc.setTarget(p.minTargetInvCond, 1.);
+  ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree,
+                                                                Patch::LS_MINEDGELENGTH);
+  ObjContribIdealJac<ObjContribFuncBarrierMovMin> minIdealJacBarFunc(1.);
+  minIdealJacBarFunc.setTarget(p.minTargetIdealJac, 1.);
 
-  MeshOptParameters::PassParameters minJacPass;
+  MeshOptPass minJacPass;
   minJacPass.barrierIterMax = p.optPassMax;
   minJacPass.optIterMax = p.itMax;
   minJacPass.contrib.push_back(&nodeDistFunc);
-//  minJacPass.contrib.push_back(&minNCJBarFunc);
-  minJacPass.contrib.push_back(&minInvCondBarFunc);
+  minJacPass.contrib.push_back(&minIdealJacBarFunc);
   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();
-  p.minInvCond = minInvCondBarFunc.getMin();
-  p.maxInvCond = minInvCondBarFunc.getMax();
-
-  Msg::StatusBar(true, "Done optimizing high order mesh (%g s)", p.CPU);
+  p.minIdealJac = minIdealJacBarFunc.getMin();
+  p.maxIdealJac = minIdealJacBarFunc.getMax();
 }

contrib/MeshQualityOptimizer/MeshQualityOptimizer.h

@@ -33,8 +33,8 @@
 class GModel;
 
 struct MeshQualOptParameters {
-//  double minTargetNCJ, maxTargetNCJ;
-  double minTargetInvCond;
+  bool excludeHex, excludePrism;
+  double minTargetIdealJac;
   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
@@ -52,16 +52,15 @@ 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 minInvCond, maxInvCond; // after optimization, range of jacobians
+  double minIdealJac, maxIdealJac; // after optimization, range of jacobians
   double CPU; // Time for optimization
 
   MeshQualOptParameters ()
-//    : minTargetNCJ(0.5), maxTargetNCJ(1.5), weightFixed(1000.),
-    : minTargetInvCond(0.5), weightFixed(1000.),
+    : excludeHex(false), excludePrism(false), minTargetIdealJac(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.)
+      maxAdaptBlob(3), adaptBlobLayerFact(2.), adaptBlobDistFact(2.), CPU(0.),
+      minIdealJac(0.), maxIdealJac(0.), SUCCESS(-1)
   {
   }
 };