From acb56d9865c6f34f25288cbf8abd793764c00332 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Thu, 3 Oct 2013 23:22:37 +0000
Subject: [PATCH] Added capability to fix invalid p1 surface meshes to OptHOM.
Changed 1-by-1 optimization strategy (to be cleaned up).
---
Fltk/highOrderToolsWindow.cpp | 1 +
Geo/MElement.cpp | 23 ++-
Geo/MElement.h | 2 +-
contrib/HighOrderMeshOptimizer/OptHOM.cpp | 13 +-
contrib/HighOrderMeshOptimizer/OptHOM.h | 3 +-
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp | 44 ++++-
contrib/HighOrderMeshOptimizer/OptHomMesh.h | 6 +-
contrib/HighOrderMeshOptimizer/OptHomRun.cpp | 169 +++++++++++-------
contrib/HighOrderMeshOptimizer/OptHomRun.h | 3 +-
9 files changed, 182 insertions(+), 82 deletions(-)
diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index 43c6093393..63b4feb541 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -151,6 +151,7 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
p.maxAdaptBlob = o->value[9]->value();
p.adaptBlobLayerFact = (int) o->value[10]->value();
p.adaptBlobDistFact = o->value[11]->value();
+ p.optPrimSurfMesh = false;
HighOrderMeshOptimizer(GModel::current(), p);
break;
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 2194c5a07c..17c2385f5a 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -112,7 +112,7 @@ double MElement::rhoShapeMeasure()
return 0.;
}
-void MElement::scaledJacRange(double &jmin, double &jmax)
+void MElement::scaledJacRange(double &jmin, double &jmax, GEntity *ge)
{
jmin = jmax = 1.0;
#if defined(HAVE_MESH)
@@ -122,6 +122,27 @@ void MElement::scaledJacRange(double &jmin, double &jmax)
getNodesCoord(nodesXYZ);
fullVector<double> SJi(numJacNodes), Bi(numJacNodes);
jac->getScaledJacobian(nodesXYZ,SJi);
+ if (ge && (ge->dim() == 2) && ge->haveParametrization()) { // If parametrized surface entity provided...
+ SVector3 geoNorm(0.,0.,0.); // ... correct Jacobian sign with geometrical normal
+ for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
+ MVertex *vert = getVertex(i);
+ GEntity *vge = vert->onWhat();
+ if (vert->onWhat() == ge) {
+ double u, v;
+ vert->getParameter(0,u);
+ vert->getParameter(1,v);
+ geoNorm += ((GFace*)ge)->normal(SPoint2(u,v));
+ }
+ }
+ if (geoNorm.normSq() == 0.) { // If no vertex on surface or average is zero, take normal at barycenter
+ SPoint2 param = ((GFace*)ge)->parFromPoint(barycenter(true),false);
+ geoNorm = ((GFace*)ge)->normal(param);
+ }
+ fullMatrix<double> elNorm(1,3);
+ jac->getPrimNormal2D(nodesXYZ,elNorm);
+ const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2);
+ if (scal < 0.) SJi.scale(-1.);
+ }
jac->lag2Bez(SJi,Bi);
jmin = *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
jmax = *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 8955f835f4..5e50155ca2 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -188,7 +188,7 @@ class MElement
return jmin;
}
virtual double angleShapeMeasure() { return 1.0; }
- virtual void scaledJacRange(double &jmin, double &jmax);
+ virtual void scaledJacRange(double &jmin, double &jmax, GEntity *ge = 0);
// get the radius of the inscribed circle/sphere if it exists,
// otherwise get the minimum radius of all the circles/spheres
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 45db5d249c..c7fb2a88a3 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -67,9 +67,10 @@ static inline double compute_f1(double v, double barrier)
// else return powP*pow(v-1.,powP-1.);
}
-OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix,
+OptHOM::OptHOM(const std::map<MElement*,GEntity*> &element2entity,
+ const std::set<MElement*> &els, std::set<MVertex*> & toFix,
bool fixBndNodes, bool fastJacEval) :
- mesh(els, toFix, fixBndNodes, fastJacEval)
+ mesh(element2entity, els, toFix, fixBndNodes, fastJacEval)
{
_optimizeMetricMin = false;
}
@@ -318,7 +319,13 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min,
std::vector<double> dSq(mesh.nVert());
mesh.distSqToStraight(dSq);
const double maxDSq = *max_element(dSq.begin(),dSq.end());
- invLengthScaleSq = 1./maxDSq; // Length scale for non-dimensional distance
+ if (maxDSq < 1.e-10) { // Length scale for non-dim. distance
+ std::vector<double> sSq(mesh.nEl());
+ mesh.elSizeSq(sSq);
+ const double maxSSq = *max_element(sSq.begin(),sSq.end());
+ invLengthScaleSq = 1./maxSSq;
+ }
+ else invLengthScaleSq = 1./maxDSq;
// Set initial guess
alglib::real_1d_array x;
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index debbc6a856..098d829132 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -43,7 +43,8 @@ class OptHOM
{
public:
Mesh mesh;
- OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix,
+ OptHOM(const std::map<MElement*,GEntity*> &element2entity,
+ const std::set<MElement*> &els, std::set<MVertex*> & toFix,
bool fixBndNodes, bool fastJacEval = false);
// returns 1 if the mesh has been optimized with success i.e. all jacobians
// are in the range; returns 0 if the mesh is valid (all jacobians positive,
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 727d645a4b..4b19087dde 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -35,7 +35,8 @@
#include "ParamCoord.h"
#include "OptHomMesh.h"
-Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix,
+Mesh::Mesh(const std::map<MElement*,GEntity*> &element2entity,
+ const std::set<MElement*> &els, std::set<MVertex*> &toFix,
bool fixBndNodes, bool fastJacEval) :
_fastJacEval(fastJacEval)
{
@@ -105,7 +106,7 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix,
// Jacobians of 3D elements
if (_dim == 2) {
_scaledNormEl.resize(nEl());
- for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(iEl);
+ for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
}
else {
_invStraightJac.resize(nEl(),1.);
@@ -116,24 +117,41 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix,
}
-void Mesh::calcScaledNormalEl2D(int iEl)
+void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl)
{
const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
fullMatrix<double> primNodesXYZ(jac->getNumPrimMapNodes(),3);
+ SVector3 geoNorm(0.,0.,0.);
+ std::map<MElement*,GEntity*>::const_iterator itEl2ent = element2entity.find(_el[iEl]);
+ GEntity *ge = (itEl2ent == element2entity.end()) ? 0 : itEl2ent->second;
+ const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
const int &iV = _el2V[iEl][i];
primNodesXYZ(i,0) = _xyz[iV].x();
primNodesXYZ(i,1) = _xyz[iV].y();
primNodesXYZ(i,2) = _xyz[iV].z();
+ if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) {
+ double u, v;
+ _vert[iV]->getParameter(0,u);
+ _vert[iV]->getParameter(1,v);
+ geoNorm += ((GFace*)ge)->normal(SPoint2(u,v));
+ }
+ }
+ if (hasGeoNorm && (geoNorm.normSq() == 0.)) {
+ SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false);
+ geoNorm = ((GFace*)ge)->normal(param);
}
- _scaledNormEl[iEl].resize(1,3);
- const double norm = jac->getPrimNormal2D(primNodesXYZ,_scaledNormEl[iEl]);
-
- _scaledNormEl[iEl](0,0) /= norm; // Re-scaling normal here is faster than an
- _scaledNormEl[iEl](0,1) /= norm; // extra scaling operation on the Jacobian
- _scaledNormEl[iEl](0,2) /= norm;
+ fullMatrix<double> &elNorm = _scaledNormEl[iEl];
+ elNorm.resize(1,3);
+ const double norm = jac->getPrimNormal2D(primNodesXYZ,elNorm);
+ double factor = 1./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
}
@@ -220,6 +238,14 @@ void Mesh::distSqToStraight(std::vector<double> &dSq)
}
}
+void Mesh::elSizeSq(std::vector<double> &sSq)
+{
+ for (int iEl = 0; iEl < nEl(); iEl++) {
+ const double s = _el[iEl]->getOuterRadius();
+ sSq[iEl] = s*s;
+ }
+}
+
void Mesh::updateGEntityPositions()
{
for (int iV = 0; iV < nVert(); iV++)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 06876da687..481cc1f96e 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -41,7 +41,8 @@
class Mesh
{
public:
- Mesh(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes, bool fastJacEval);
+ Mesh(const std::map<MElement*,GEntity*> &element2entity,
+ const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes, bool fastJacEval);
inline const int &nPC() { return _nPC; }
inline int nVert() { return _vert.size(); }
@@ -67,6 +68,7 @@ public:
void getUvw(double *it);
void updateMesh(const double *it);
void distSqToStraight(std::vector<double> &dSq);
+ void elSizeSq(std::vector<double> &sSq);
void updateGEntityPositions();
void writeMSH(const char *filename);
@@ -109,7 +111,7 @@ private:
int addVert(MVertex* vert);
int addFreeVert(MVertex* vert, const int iV, const int nPCV,
std::set<MVertex*> &toFix);
- void calcScaledNormalEl2D(int iEl);
+ void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl);
static inline int indJB2DBase(int nNod, int l, int i, int j)
{
return (l*nNod+i)*nNod+j;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 9edc393daf..8c0b20135c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -163,11 +163,13 @@ static std::set<MVertex *> getAllBndVertices
static std::set<MElement*> getSurroundingBlob
(MElement *el, int depth,
const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
- const double distFactor, int forceDepth = 0)
+ const double distFactor, int forceDepth, bool optPrimSurfMesh)
{
- const SPoint3 p = el->barycenter_infty();
- const double limDist = distMaxStraight(el) * distFactor;
+ const SPoint3 p = el->barycenter();
+ const double dist = distMaxStraight(el);
+ const double limDist = ((optPrimSurfMesh && (dist < 1.e-10)) ?
+ el->getOuterRadius() : dist) * distFactor;
std::set<MElement*> blob;
std::list<MElement*> currentLayer, lastLayer;
@@ -220,10 +222,18 @@ static void calcVertex2Elements(int dim, GEntity *entity,
}
}
+static void calcElement2Entity(GEntity *entity, std::map<MElement*,GEntity*> &element2entity)
+{
+ for (size_t i = 0; i < entity->getNumMeshElements(); ++i) {
+ MElement *element = entity->getMeshElement(i);
+ element2entity.insert(std::pair<MElement*,GEntity*>(element,entity));
+ }
+}
+
static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConnectedBlobs
(const std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
const std::set<MElement*> &badElements, int depth, const double distFactor,
- bool weakMerge = false)
+ bool weakMerge, bool optPrimSurfMesh)
{
Msg::Info("Starting blob generation from %i bad elements...", badElements.size());
@@ -232,9 +242,11 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
Msg::Info("Constructing %i primary blobs", badElements.size());
std::vector<std::set<MElement*> > primBlobs;
primBlobs.reserve(badElements.size());
- for (std::set<MElement*>::const_iterator it = badElements.begin();
- it != badElements.end(); ++it)
- primBlobs.push_back(getSurroundingBlob(*it, depth, vertex2elements, distFactor));
+ for (std::set<MElement*>::const_iterator it = badElements.begin(); it != badElements.end(); ++it) {
+ const int minLayers = ((*it)->getDim() == 3) ? 1 : 0;
+ primBlobs.push_back(getSurroundingBlob(*it, depth, vertex2elements,
+ distFactor, minLayers, optPrimSurfMesh));
+ }
// Compute blob connectivity
Msg::Info("Computing blob connectivity...");
@@ -292,6 +304,7 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
static void optimizeConnectedBlobs
(const std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+ const std::map<MElement*,GEntity*> &element2entity,
std::set<MElement*> &badasses,
OptHomParameters &p, int samples, bool weakMerge = false)
{
@@ -300,14 +313,15 @@ static void optimizeConnectedBlobs
p.maxJac = -1.e100;
std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > toOptimize =
- getConnectedBlobs(vertex2elements, badasses, p.nbLayers, p.distanceFactor, weakMerge);
+ getConnectedBlobs(vertex2elements, badasses, p.nbLayers,
+ p.distanceFactor, weakMerge, p.optPrimSurfMesh);
//#pragma omp parallel for schedule(dynamic, 1)
for (int i = 0; i < toOptimize.size(); ++i) {
Msg::Info("Optimizing a blob %i/%i composed of %4d elements", i+1,
toOptimize.size(), toOptimize[i].first.size());
fflush(stdout);
- OptHOM temp(toOptimize[i].first, toOptimize[i].second, p.fixBndNodes);
+ OptHOM temp(element2entity, toOptimize[i].first, toOptimize[i].second, p.fixBndNodes);
//std::ostringstream ossI1;
//ossI1 << "initial_ITER_" << i << ".msh";
//temp.mesh.writeMSH(ossI1.str().c_str());
@@ -460,6 +474,7 @@ static bool detectNewBrokenElement(std::set<MElement*> &layer,
static void optimizeOneByOne
(const std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+ const std::map<MElement*,GEntity*> &element2entity,
std::set<MElement*> badElts, OptHomParameters &p, int samples)
{
p.SUCCESS = 1;
@@ -484,60 +499,84 @@ static void optimizeOneByOne
for (int iterBlob=0; iterBlob<p.maxAdaptBlob; iterBlob++) {
- OptHOM *opt;
+// OptHOM *opt;
+//
+// // First step: small blob with unsafe optimization (only 1st-order
+// // bnd. vertices fixed)
+// std::set<MElement*> toOptimizePrim = getSurroundingBlob
+// (worstEl, nbLayers, vertex2elements, distanceFactor, 0);
+// std::set<MVertex*> toFix1 = getPrimBndVertices(toOptimizePrim, vertex2elements);
+// std::set<MElement*> toOptimize1;
+// std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
+// badElts.begin(),badElts.end(), // Do not optimize badElts
+// std::inserter(toOptimize1, toOptimize1.end()));
+// Msg::Info("Optimizing primary blob %i (max. %i remaining) composed of"
+// " %4d elements", iBadEl, badElts.size(), toOptimize1.size());
+// fflush(stdout);
+// opt = new OptHOM(toOptimize1, toFix1, p.fixBndNodes);
+// //std::ostringstream ossI1;
+// //ossI1 << "initial_primary_" << iter << ".msh";
+// //opt->mesh.writeMSH(ossI1.str().c_str());
+// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
+// false, samples, p.itMax, p.optPassMax);
+//
+// // Second step: add external layer, check it and optimize it safely (all
+// // bnd. vertices fixed) if new broken element
+// if (success > 0) {
+// opt->mesh.updateGEntityPositions();
+// std::set<MElement*> layer = addBlobLayer(toOptimizePrim, vertex2elements);
+// if (detectNewBrokenElement(layer, badElts, p)) {
+// delete opt;
+// std::set<MVertex*> toFix2 = getAllBndVertices(toOptimizePrim, vertex2elements);
+// std::set<MElement*> toOptimize2;
+// std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
+// badElts.begin(),badElts.end(),
+// std::inserter(toOptimize2, toOptimize2.end()));
+// Msg::Info("Optimizing corrective blob %i (max. %i remaining) "
+// "composed of %4d elements", iBadEl, badElts.size(),
+// toOptimize2.size());
+// fflush(stdout);
+// opt = new OptHOM(toOptimize2, toFix2, p.fixBndNodes);
+// //std::ostringstream ossI1;
+// //ossI1 << "initial_corrective_" << iter << ".msh";
+// //opt->mesh.writeMSH(ossI1.str().c_str());
+// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
+// p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
+// if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
+// Msg::Info("Jacobian optimization succeed, starting svd optimization");
+// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
+// p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
+// }
+// }
+// else {
+// Msg::Info("Primary blob %i did not break new elements, "
+// "no correction needed", iBadEl);
+// fflush(stdout);
+// }
+// }
- // First step: small blob with unsafe optimization (only 1st-order
- // bnd. vertices fixed)
std::set<MElement*> toOptimizePrim = getSurroundingBlob
- (worstEl, nbLayers, vertex2elements, distanceFactor, 0);
- std::set<MVertex*> toFix1 = getPrimBndVertices(toOptimizePrim, vertex2elements);
- std::set<MElement*> toOptimize1;
+ (worstEl, nbLayers, vertex2elements, distanceFactor, 1, p.optPrimSurfMesh);
+// std::set<MElement*> layer = addBlobLayer(toOptimizePrim, vertex2elements);
+ std::set<MVertex*> toFix = getAllBndVertices(toOptimizePrim, vertex2elements);
+ std::set<MElement*> toOptimize;
std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
- badElts.begin(),badElts.end(), // Do not optimize badElts
- std::inserter(toOptimize1, toOptimize1.end()));
- Msg::Info("Optimizing primary blob %i (max. %i remaining) composed of"
- " %4d elements", iBadEl, badElts.size(), toOptimize1.size());
+ badElts.begin(),badElts.end(),
+ std::inserter(toOptimize, toOptimize.end()));
+ Msg::Info("Optimizing blob %i (max. %i remaining) "
+ "composed of %4d elements", iBadEl, badElts.size(),
+ toOptimize.size());
fflush(stdout);
- opt = new OptHOM(toOptimize1, toFix1, p.fixBndNodes);
+ OptHOM *opt = new OptHOM(element2entity, toOptimize, toFix, p.fixBndNodes);
//std::ostringstream ossI1;
- //ossI1 << "initial_primary_" << iter << ".msh";
+ //ossI1 << "initial_corrective_" << iter << ".msh";
//opt->mesh.writeMSH(ossI1.str().c_str());
- success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
- false, samples, p.itMax, p.optPassMax);
-
- // Second step: add external layer, check it and optimize it safely (all
- // bnd. vertices fixed) if new broken element
- if (success > 0) {
- opt->mesh.updateGEntityPositions();
- std::set<MElement*> layer = addBlobLayer(toOptimizePrim, vertex2elements);
- if (detectNewBrokenElement(layer, badElts, p)) {
- delete opt;
- std::set<MVertex*> toFix2 = getAllBndVertices(toOptimizePrim, vertex2elements);
- std::set<MElement*> toOptimize2;
- std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
- badElts.begin(),badElts.end(),
- std::inserter(toOptimize2, toOptimize2.end()));
- Msg::Info("Optimizing corrective blob %i (max. %i remaining) "
- "composed of %4d elements", iBadEl, badElts.size(),
- toOptimize2.size());
- fflush(stdout);
- opt = new OptHOM(toOptimize2, toFix2, p.fixBndNodes);
- //std::ostringstream ossI1;
- //ossI1 << "initial_corrective_" << iter << ".msh";
- //opt->mesh.writeMSH(ossI1.str().c_str());
- success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
- p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
- if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
- Msg::Info("Jacobian optimization succeed, starting svd optimization");
- success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
- p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
- }
- }
- else {
- Msg::Info("Primary blob %i did not break new elements, "
- "no correction needed", iBadEl);
- fflush(stdout);
- }
+ success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
+ p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
+ if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
+ Msg::Info("Jacobian optimization succeed, starting svd optimization");
+ success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
+ p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
}
// Measure min and max Jac., update mesh
@@ -558,11 +597,11 @@ static void optimizeOneByOne
nbLayers *= p.adaptBlobLayerFact;
Msg::Info("Blob %i failed (adapt #%i), adapting with increased size",
iBadEl, iterBlob);
- if (iterBlob == p.maxAdaptBlob-1) {
+// if (iterBlob == p.maxAdaptBlob-1) {
std::ostringstream ossI2;
ossI2 << "final_" << iBadEl << ".msh";
opt->mesh.writeMSH(ossI2.str().c_str());
- }
+// }
}
else break;
}
@@ -589,18 +628,20 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
gm->getEntities(entities);
std::map<MVertex*, std::vector<MElement *> > vertex2elements;
+ std::map<MElement*,GEntity*> element2entity;
std::set<MElement*> badasses;
for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
GEntity* &entity = entities[iEnt];
if (entity->dim() != p.dim || (p.onlyVisible && !entity->getVisibility())) continue;
Msg::Info("Computing connectivity and bad elements for entity %d...",
entity->tag());
- calcVertex2Elements(p.dim,entity,vertex2elements); // Compute vert. -> elt. connectivity
+ calcVertex2Elements(p.dim,entity,vertex2elements);
+ if (p.optPrimSurfMesh) calcElement2Entity(entity,element2entity);
for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
double jmin, jmax;
MElement *el = entity->getMeshElement(iEl);
if (el->getDim() == p.dim) {
- el->scaledJacRange(jmin, jmax);
+ el->scaledJacRange(jmin, jmax, p.optPrimSurfMesh ? entity : 0);
if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) badasses.insert(el);
}
@@ -608,11 +649,11 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
}
if (p.strategy == 0)
- optimizeConnectedBlobs(vertex2elements, badasses, p, samples, false);
+ optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, false);
else if (p.strategy == 2)
- optimizeConnectedBlobs(vertex2elements, badasses, p, samples, true);
+ optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, true);
else
- optimizeOneByOne(vertex2elements, badasses, p, samples);
+ optimizeOneByOne(vertex2elements, element2entity, badasses, p, samples);
if (p.SUCCESS == 1)
Msg::Info("Optimization succeeded");
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index 865301c078..f82751da05 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -53,6 +53,7 @@ struct OptHomParameters {
int maxAdaptBlob; // Max. nb. of blob adaptation interations
int adaptBlobLayerFact; // Growth factor in number of layers for blob adaptation
double adaptBlobDistFact; // Growth factor in distance factor for blob adaptation
+ bool optPrimSurfMesh; // Enable optimisation of p1 surface meshes
// OUTPUT ------>
int SUCCESS ; // 0 --> success , 1 --> Not converged
@@ -63,7 +64,7 @@ struct OptHomParameters {
: BARRIER_MIN_METRIC(-1.), BARRIER_MIN(0.1), BARRIER_MAX(2.0), weightFixed(1.e6),
weightFree (1.e2), nbLayers (6) , dim(3) , itMax(300), onlyVisible(true),
distanceFactor(12), fixBndNodes(false), strategy(0), maxAdaptBlob(3),
- adaptBlobLayerFact(2.), adaptBlobDistFact(2.)
+ adaptBlobLayerFact(2.), adaptBlobDistFact(2.), optPrimSurfMesh(false)
{
}
};
--
GitLab