diff --git a/CMakeLists.txt b/CMakeLists.txt
index cb478406c11f61f9d6a1cb0d78591d908c9da4a9..ff4521c23cbe464a9306a9858cb313a46cf4dccd 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 94b970e882ace3609b6a80c130ee0546d4cdf2c5..0417c6ba336157cc33e970d9d8150d694b9cc3ad 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/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);
diff --git a/contrib/MeshOptimizer/CMakeLists.txt b/contrib/MeshOptimizer/CMakeLists.txt
index 525303ecd0a77e2f925e492e83c36a9f38cfa5db..1cee09e2ed669fee3713f2644b0c665c5aded625 100644
--- a/contrib/MeshOptimizer/CMakeLists.txt
+++ b/contrib/MeshOptimizer/CMakeLists.txt
@@ -5,6 +5,7 @@
set(SRC
MeshOpt.cpp
+ MeshOptCommon.cpp
MeshOptimizer.cpp
MeshOptPatch.cpp
MeshOptObjectiveFunction.cpp
diff --git a/contrib/MeshOptimizer/MeshOptCommon.cpp b/contrib/MeshOptimizer/MeshOptCommon.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d32a9528fb6e99e92d3822755f843c505256b598
--- /dev/null
+++ b/contrib/MeshOptimizer/MeshOptCommon.cpp
@@ -0,0 +1,136 @@
+// 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;
+}
diff --git a/contrib/MeshOptimizer/MeshOptCommon.h b/contrib/MeshOptimizer/MeshOptCommon.h
index 9f895e84f9c5ede16ae4f6847aa42dbcae194a1e..189c7689d5030a4d4c009f794c3ebe208a2ebfce 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.h
+++ b/contrib/MeshOptimizer/MeshOptCommon.h
@@ -34,35 +34,48 @@
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 {
- 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 {
- struct { // If adaptive strategy:
- int maxAdaptPatch; // Max. nb. of adaptation iterations
- int maxLayersAdaptFact; // Growth rate in number of layers around a bad element
- double distanceAdaptFact; // Growth rate in max. distance from bad element
- };
- bool weakMerge; // If connected strategy: weak or strong merging of patches
+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 {
+ struct { // If adaptive strategy:
+ int maxAdaptPatch; // Max. nb. of adaptation iterations
+ int maxLayersAdaptFact; // Growth rate in number of layers around a bad element
+ double distanceAdaptFact; // Growth rate in max. distance from bad element
};
- 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
+ bool weakMerge; // If connected strategy: weak or strong merging of patches
};
+ 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
diff --git a/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp b/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp
index 8be99a8567af2d97a16df53205f7f2fc493edd44..cd5757325c612de14504b348183341504b061d07 100644
--- a/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp
+++ b/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp
@@ -1,6 +1,7 @@
// 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;
}
diff --git a/contrib/MeshOptimizer/MeshOptObjContribFunc.h b/contrib/MeshOptimizer/MeshOptObjContribFunc.h
index 04c7aef64f5185f541b09f629c2aba6c60153188..ec12b3fde04ff8bf572bb51e208c19364f541c0f 100644
--- a/contrib/MeshOptimizer/MeshOptObjContribFunc.h
+++ b/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);
};
diff --git a/contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h b/contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h
index 070d744447e3df5d5f24a871a6d109aea12ac8af..6aa3c7ec05f6caacf82e4975732f073dcd2ea2f3 100644
--- a/contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h
+++ b/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);
}
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index c4f71eb0620b6abce8c24e29c8b36c2683e52430..0debf3111a64cd336ddb187dd8e0060f04a80bee 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/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.;
- 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.;
- for (int iEl = 0; iEl < nEl(); iEl++) {
- const double s = el(iEl)->getOuterRadius(), ss = s*s;
- if (ss > maxSSq) maxSSq = ss;
+ 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;
+ }
+ 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 d = el(iEl)->minEdge(), dd = d*d;
+ if (dd > maxDSq) maxDSq = dd;
+ }
+ break;
}
- _invLengthScaleSq = 1./maxSSq;
}
- else _invLengthScaleSq = 1./maxDSq;
+ _invLengthScaleSq = 1./maxDSq;
}
}
@@ -266,13 +279,14 @@ 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.);
+ _invStraightJac.resize(nEl(), 1.);
double dumJac[3][3];
for (int iEl = 0; iEl < nEl(); iEl++)
_invStraightJac[iEl] = 1. / fabs(_el[iEl]->getPrimaryJacobian(0.,0.,0.,dumJac));
@@ -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,21 +334,25 @@ 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);
-// 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));
+ 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
}
@@ -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);
+ for (int iEl = 0; iEl < nEl(); iEl++)
+ calcNormalEl2D(iEl, NS_INVNORM, _IJacNormEl[iEl], true);
+ }
+ 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::NCJ(int iEl, std::vector<double> &NCJ)
+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];
- 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;
- }
- 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;
- }
+ 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();
}
-}
+ // 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";
+// }
+// }
-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];
- 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;
- }
- }
+ // 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 NCJ
+ // 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]);
- }
- _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++)
+ 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 < 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];
}
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index 4269ba833d02510802bea41b289f0f77d8f86a5b..d9b477f27770d51e26f370937b0cfb4ba3c9ef7b 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/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
- double _invLengthScaleSq; // Square inverse of a length for node displacement scaling
+ LengthScaling _typeLengthScale;
+ 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. 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
+ enum NormalScaling { NS_UNIT, NS_INVNORM, NS_SQRTNORM };
+ std::vector<int> _nBezEl; // Number of Bezier poly. for an el.
+ 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
};
diff --git a/contrib/MeshOptimizer/MeshOptimizer.cpp b/contrib/MeshOptimizer/MeshOptimizer.cpp
index 4cf44f2d1acb7a1db9e317c82b7f8558d5e51a77..ac1a3e0f799797cc309c9f3681bac5ac4584ebca 100644
--- a/contrib/MeshOptimizer/MeshOptimizer.cpp
+++ b/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,9 +91,11 @@ 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.insert(*itN).second) currentLayer.push_back(*itN); // Assume that if an el is too far, its neighbours are too far as well
+ 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
+ }
}
}
}
@@ -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(),
diff --git a/contrib/MeshQualityOptimizer/MeshQualityObjContribIdealJac.h b/contrib/MeshQualityOptimizer/MeshQualityObjContribIdealJac.h
new file mode 100644
index 0000000000000000000000000000000000000000..3f2245f0813697dbeff70027aa1942fd1a077d47
--- /dev/null
+++ b/contrib/MeshQualityOptimizer/MeshQualityObjContribIdealJac.h
@@ -0,0 +1,98 @@
+// TODO: Copyright
+
+#ifndef _MESHQUALITYOBJCONTRIBIDEALJAC_H_
+#define _MESHQUALITYOBJCONTRIBIDEALJAC_H_
+
+#include "MeshOptPatch.h"
+#include "MeshOptObjContrib.h"
+
+
+template<class FuncType>
+class ObjContribIdealJac : public ObjContrib, public FuncType
+{
+public:
+ ObjContribIdealJac(double weight);
+ virtual ~ObjContribIdealJac() {}
+ 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>
+ObjContribIdealJac<FuncType>::ObjContribIdealJac(double weight) :
+ ObjContrib("IdealJac", FuncType::getNamePrefix()+"IdealJac"),
+ _mesh(0), _weight(weight)
+{
+}
+
+
+template<class FuncType>
+ObjContrib *ObjContribIdealJac<FuncType>::copy() const
+{
+ return new ObjContribIdealJac<FuncType>(*this);
+}
+
+
+template<class FuncType>
+void ObjContribIdealJac<FuncType>::initialize(Patch *mesh)
+{
+ _mesh = mesh;
+ _mesh->initIdealJac();
+ updateMinMax();
+ FuncType::initialize(_min, _max);
+}
+
+
+template<class FuncType>
+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> 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 * gIJ[_mesh->indGIJac(iEl, l, iPC)];
+ _min = std::min(_min, iJ[l]);
+ _max = std::max(_max, iJ[l]);
+ }
+ }
+
+ return true;
+}
+
+
+template<class FuncType>
+void ObjContribIdealJac<FuncType>::updateMinMax()
+{
+ _min = BIGVAL;
+ _max = -BIGVAL;
+
+ for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ 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 /* _MESHQUALITYOBJCONTRIBIDEALJAC_H_ */
diff --git a/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h b/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
index ccdac14fcf1f1566c0a5f46dbca406dbdb12c807..65b758fb3f0d91c9f8e5daa85bea529845d39a6d 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
+++ b/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]);
diff --git a/contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h b/contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h
deleted file mode 100644
index 0673a2812825d938ca7baa0adbd1cda6021cfd81..0000000000000000000000000000000000000000
--- a/contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h
+++ /dev/null
@@ -1,97 +0,0 @@
-// 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->nNCJEl(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
index c50d513767ae7ea9c2726e4501a221036e21a95f..db96312c383c825b9d9800d987aaa2737a6c4934 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
+++ b/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();
}
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
index d34d1083d33b01a16617e1d57d5402a836d6c753..36722f0a453d634a96f480992703533d7841cbe8 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
@@ -33,15 +33,15 @@
class GModel;
struct MeshQualOptParameters {
-// 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
- 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
+ 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
+ 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
@@ -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)
{
}
};