diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index 236d257873219ddbd1142411b601e44dfe5086d5..566a5b5e3616f69b9e6218aa460e7fbe320bb530 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -96,20 +96,33 @@ static void chooseopti_cb(Fl_Widget *w, void *data)
if (elastic == 1){
o->choice[0]->deactivate();
+ o->choice[3]->deactivate();
for (int i=3;i<=6;i++)
o->value[i]->deactivate();
+ for (int i=9;i<=11;i++) o->value[i]->deactivate();
// o->push[1]->deactivate();
}
else {
o->push[0]->activate();
o->choice[0]->activate();
+ o->choice[3]->activate();
for (int i=3;i<=6;i++)
o->value[i]->activate();
+ for (int i=9;i<=11;i++) o->value[i]->activate();
// o->push[1]->activate();
}
#endif
}
+static void chooseopti_strategy(Fl_Widget *w, void *data)
+{
+#if defined(HAVE_OPTHOM)
+ highOrderToolsWindow *o = FlGui::instance()->highordertools;
+ if (o->choice[3]->value() == 0) for (int i=9;i<=11;i++) o->value[i]->deactivate();
+ else for (int i=9;i<=11;i++) o->value[i]->activate();
+#endif
+}
+
static void highordertools_runelas_cb(Fl_Widget *w, void *data)
{
#if defined(HAVE_OPTHOM)
@@ -133,8 +146,12 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data)
p.optPassMax = (int) o->value[4]->value();
p.weightFixed = o->value[5]->value();
p.weightFree = o->value[6]->value();
- p.DistanceFactor = o->value[7]->value();
- p.method = o->CAD ? (int)o->choice[0]->value() : 2;
+ p.distanceFactor = o->value[7]->value();
+ p.fixBndNodes = (!o->CAD) || (o->choice[0]->value() == 0);
+ p.strategy = o->choice[3]->value();
+ p.maxAdaptBlob = o->value[9]->value();
+ p.adaptBlobLayerFact = (int) o->value[10]->value();
+ p.adaptBlobDistFact = o->value[11]->value();
HighOrderMeshOptimizer (GModel::current(),p);
printf("CPU TIME = %4f seconds\n",p.CPU);
}
@@ -149,7 +166,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
FL_NORMAL_SIZE -= deltaFontSize;
int width = 3 * IW + 4 * WB;
- int height = 25 * BH;
+ int height = 28 * BH;
win = new paletteWindow
(width, height, CTX::instance()->nonModalWindows ? true : false, "High Order Tools");
@@ -315,6 +332,48 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
value[4]->align(FL_ALIGN_RIGHT);
value[4]->value(50);
+ static Fl_Menu_Item menu_strategy[] = {
+ {"Connected blobs", 0, 0, 0},
+ {"Adaptive one-by-one", 0, 0, 0},
+ {0}
+ };
+
+ y += BH;
+ choice[3] = new Fl_Choice
+ (x,y, IW, BH, "Strategy");
+ choice[3]->menu(menu_strategy);
+ choice[3]->align(FL_ALIGN_RIGHT);
+ choice[3]->callback(chooseopti_strategy);
+
+ y += BH;
+ value[9] = new Fl_Value_Input
+ (x,y, IW, BH, "Max. number of blob adaptation iter.");
+ value[9]->minimum(1);
+ value[9]->maximum(100);
+ value[9]->step(1);
+ value[9]->align(FL_ALIGN_RIGHT);
+ value[9]->value(2);
+ value[9]->deactivate();
+
+ y += BH;
+ value[10] = new Fl_Value_Input
+ (x,y, IW, BH, "Num. layer adaptation factor");
+ value[10]->align(FL_ALIGN_RIGHT);
+ value[10]->minimum(1);
+ value[10]->maximum(100);
+ value[10]->step(1);
+ value[10]->value(2);
+ value[10]->deactivate();
+
+ y += BH;
+ value[11] = new Fl_Value_Input
+ (x,y, IW, BH, "Distance adaptation factor");
+ value[11]->align(FL_ALIGN_RIGHT);
+ value[11]->minimum(1.);
+ value[11]->maximum(100.);
+ value[11]->value(2.);
+ value[11]->deactivate();
+
y += 1.5*BH;
push[1] = new Fl_Button
(x,y, IW, BH, "Apply");
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index c2d7c049959265721f46ba13e79e3361f0f1b979..60fe3e071641403cd3b6d43f99bf9d7133d97efb 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -221,18 +221,17 @@ class MPyramidN : public MPyramid {
protected:
std::vector<MVertex*> _vs;
const char _order;
- double _disto;
public:
MPyramidN(MVertex* v0, MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
const std::vector<MVertex*> &v, char order, int num=0, int part=0)
- : MPyramid(v0, v1, v2, v3, v4, num, part), _vs(v), _order(order), _disto(-1e22)
+ : MPyramid(v0, v1, v2, v3, v4, num, part), _vs(v), _order(order)
{
for (unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
getFunctionSpace(order);
}
MPyramidN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
- : MPyramid(v[0], v[1], v[2], v[3], v[4], num, part), _order(order), _disto(-1e22)
+ : MPyramid(v[0], v[1], v[2], v[3], v[4], num, part), _order(order)
{
for (unsigned int i = 5; i < v.size(); i++ ) _vs.push_back(v[i]);
for (unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 41ecbf88fe8d8983e6e3107495cff1bc722a4d4a..e722f24a545dc4a9073f5c68bf826e5bc84417ef 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -42,11 +42,11 @@ static inline double compute_f1(double v, double barrier)
// Constructor
-OptHOM::OptHOM(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method) :
- mesh(ge, els, toFix, method)
+OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes) :
+ mesh(els, toFix, fixBndNodes)
{
_optimizeMetricMin = false;
-};
+}
@@ -241,6 +241,8 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
static int OPTMETHOD = 1;
Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar);
+ std::cout << "--- Optimization pass with initial jac. range ("
+ << minJac << "," << maxJac << "), jacBar = " << jacBar << "\n";
iter = 0;
@@ -321,7 +323,6 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
jacBar = jacBarStart;
- setBarrierTerm(jacBarStart);
_optimizeBarrierMax = false;
// Calculate initial objective function value and gradient
@@ -337,24 +338,37 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
<< " and max. barrier = " << barrier_max << std::endl;
int ITER = 0;
+ bool minJacOK = true;
while (minJac < barrier_min) {
+ const double startMinJac = minJac;
OptimPass(x, gradObj, itMax);
recalcJacDist();
jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
- setBarrierTerm(jacBar);
- if (ITER ++ > optPassMax) break;
+ if (ITER ++ > optPassMax) {
+ minJacOK = (minJac > barrier_min);
+ break;
+ }
+ if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
+ std::cout << "Stagnation in minJac detected, stopping optimization\n";
+ minJacOK = false;
+ break;
+ }
}
- if (!_optimizeMetricMin) {
+ ITER = 0;
+ if (minJacOK && (!_optimizeMetricMin)) {
_optimizeBarrierMax = true;
jacBar = 1.1 * maxJac;
- setBarrierTerm(jacBar);
while (maxJac > barrier_max ) {
+ const double startMaxJac = maxJac;
OptimPass(x, gradObj, itMax);
recalcJacDist();
jacBar = 1.1 * maxJac;
- setBarrierTerm(jacBar);
if (ITER ++ > optPassMax) break;
+ if (fabs((maxJac-startMaxJac)/startMaxJac) < 0.01) {
+ std::cout << "Stagnation in maxJac detected, stopping optimization\n";
+ break;
+ }
}
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index 2d6fcda080b0504318f649dc2fa5ad25209d8e8b..53e60216ec7084353b8fc2d0cb86f5412f90bde2 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -21,7 +21,7 @@ public:
Mesh mesh;
- OptHOM(GEntity *gf, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method);
+ OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes);
// 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, JMIN > 0) but JMIN < barrier_min || JMAX > barrier_max
// returns -1 if the mesh is invalid : some jacobians cannot be made positive
@@ -36,7 +36,6 @@ public:
private:
-// double lambda, lambda2, powM, powP, invLengthScaleSq;
double lambda, lambda2, jacBar, invLengthScaleSq;
int iter, progressInterv; // Current iteration, interval of iterations for reporting
bool _optimizeMetricMin;
@@ -45,7 +44,6 @@ private:
bool _optimizeBarrierMax; // false : only moving barrier min; true : fixed barrier min + moving barrier max
- inline void setBarrierTerm(double jacBarrier) {jacBar = jacBarrier;}
bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
bool addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj);
bool addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 08398951b89659e78805abe5d62e3a90f43b78e6..ca60318901b226100e15407f34c12c66b81a7318 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -8,32 +8,21 @@
-Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFix, int method) :
- _ge(ge)
+Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBndNodes)
{
- _dim = _ge->dim();
- if (method & METHOD_PHYSCOORD) {
+ _dim = (*els.begin())->getDim();
+
+ if (fixBndNodes) {
if (_dim == 2) _pc = new ParamCoordPhys2D;
else _pc = new ParamCoordPhys3D;
- Msg::Debug("METHOD: Using physical coordinates");
- }
- else if (method & METHOD_SURFCOORD) {
- if (_dim == 2) {
- _pc = new ParamCoordSurf(_ge);
- Msg::Debug("METHOD: Using surface parametric coordinates");
- }
- else Msg::Error("ERROR: Surface parametric coordinates only for 2D optimization");
+ Msg::Debug("METHOD: Fixing boundary nodes and using physical coordinates");
}
else {
_pc = new ParamCoordParent;
- Msg::Debug("METHOD: Using parent parametric coordinates");
+ Msg::Debug("METHOD: Freeing boundary nodes and using parent parametric coordinates");
}
- if (method & METHOD_RELAXBND)Msg::Debug("METHOD: Relaxing boundary vertices");
- else if (method & METHOD_FIXBND) Msg::Debug("METHOD: Fixing all boundary vertices");
- else Msg::Debug("METHOD: Fixing vertices on geometric points and \"toFix\" boundary");
-
// Initialize elements, vertices, free vertices and element->vertices connectivity
const int nElements = els.size();
_nPC = 0;
@@ -55,8 +44,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
_el2V[iEl].push_back(iV);
const int nPCV = _pc->nCoord(vert);
bool isFV = false;
- if (method & METHOD_RELAXBND) isFV = true;
- else if (method & METHOD_FIXBND) isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
+ if (fixBndNodes) isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
else isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
if (isFV) {
int iFV = addFreeVert(vert,iV,nPCV,toFix);
@@ -87,14 +75,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
double dumJac[3][3];
for (int iEl = 0; iEl < nEl(); iEl++) _invStraightJac[iEl] = 1. / _el[iEl]->getPrimaryJacobian(0.,0.,0.,dumJac);
}
- if ((_dim == 2) && (method & METHOD_PROJJAC)) {
- projJac = true;
- Msg::Debug("METHOD: Using projected Jacobians");
- }
- else {
- projJac = false;
- Msg::Debug("METHOD: Using usual Jacobians");
- }
+
}
@@ -371,7 +352,7 @@ void Mesh::writeMSH(const char *filename)
fprintf(f, "%d\n", nEl());
for (int iEl = 0; iEl < nEl(); iEl++) {
// MElement *MEl = _el[iEl];
- fprintf(f, "%d %d 2 0 %d", iEl+1, _el[iEl]->getTypeForMSH(), _ge->tag());
+ fprintf(f, "%d %d 2 0 0", iEl+1, _el[iEl]->getTypeForMSH());
for (size_t iVEl = 0; iVEl < _el2V[iEl].size(); iVEl++) fprintf(f, " %d", _el2V[iEl][iVEl] + 1);
fprintf(f, "\n");
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 0a716ae5a2760fc45126aaf2e3b4a21289d7e26a..481d6a2b8e292c1e0268626dea52378135efbd44 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -17,7 +17,7 @@ class Mesh
public:
- Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method);
+ Mesh(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes);
inline const int &nPC() { return _nPC; }
inline int nVert() { return _vert.size(); }
@@ -46,9 +46,6 @@ public:
void updateGEntityPositions();
void writeMSH(const char *filename);
- enum { METHOD_RELAXBND = 1 << 0, METHOD_FIXBND = 1 << 1, METHOD_PHYSCOORD = 1 << 2,
- METHOD_SURFCOORD = 1 << 3, METHOD_PROJJAC = 1 << 4 };
-
// inline double xyzDBG(int iV, int iCoord) { return _xyz[iV][iCoord]; }
// inline double ixyzDBG(int iV, int iCoord) { return _ixyz[iV][iCoord]; }
// inline double sxyzDBG(int iV, int iCoord) { return _sxyz[iV][iCoord]; }
@@ -58,7 +55,6 @@ public:
private:
- GEntity *_ge;
int _dim;
int _nPC; // Total nb. of parametric coordinates
@@ -77,7 +73,6 @@ private:
std::vector<std::vector<int> > _indPCEl; // Index of parametric coord. for an el.
ParamCoord *_pc;
- bool projJac; // Using "projected" Jacobians or not
int addVert(MVertex* vert);
int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index b713254b1e1bda8282e80ea24da9ee8e1f0d97c1..0cc9e087fa984893d8c50a81d651205d42fab05b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -138,12 +138,13 @@ void exportMeshToDassault (GModel *gm, const std::string &fn, int dim){
-static std::set<MVertex *> getBndVertices(std::set<MElement*> &elements, std::map<MVertex *, std::vector<MElement*> > &vertex2elements)
+static std::set<MVertex *> getAllBndVertices(std::set<MElement*> &elements,
+ const std::map<MVertex*, std::vector<MElement*> > &vertex2elements)
{
std::set<MVertex*> bnd;
for (std::set<MElement*>::iterator itE = elements.begin(); itE != elements.end(); ++itE) {
for (int i = 0; i < (*itE)->getNumPrimaryVertices(); ++i) {
- const std::vector<MElement*> &neighbours = vertex2elements[(*itE)->getVertex(i)];
+ const std::vector<MElement*> &neighbours = vertex2elements.find((*itE)->getVertex(i))->second;
for (size_t k = 0; k < neighbours.size(); ++k) {
if (elements.find(neighbours[k]) == elements.end()) {
for (int j = 0; j < neighbours[k]->getNumVertices(); ++j) {
@@ -158,7 +159,9 @@ static std::set<MVertex *> getBndVertices(std::set<MElement*> &elements, std::ma
-static std::set<MElement*> getSurroundingBlob(MElement *el, int depth, std::map<MVertex *, std::vector<MElement*> > &vertex2elements, const double distFactor)
+static std::set<MElement*> getSurroundingBlob(MElement *el, int depth,
+ const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
+ const double distFactor, int forceDepth = 0)
{
const SPoint3 p = el->barycenter_infty();
@@ -173,9 +176,9 @@ static std::set<MElement*> getSurroundingBlob(MElement *el, int depth, std::map<
currentLayer.clear();
for (std::list<MElement*>::iterator it = lastLayer.begin(); it != lastLayer.end(); ++it) {
for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
- const std::vector<MElement*> &neighbours = vertex2elements[(*it)->getVertex(i)];
+ const std::vector<MElement*> &neighbours = vertex2elements.find((*it)->getVertex(i))->second;
for (std::vector<MElement*>::const_iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN)
- if (p.distance((*itN)->barycenter_infty()) < limDist)
+ if ((d < forceDepth) || (p.distance((*itN)->barycenter_infty()) < limDist))
if (blob.insert(*itN).second) currentLayer.push_back(*itN); // Assume that if an el is too far, its neighbours are too far as well
}
}
@@ -203,20 +206,24 @@ static void addBlobChaintoGroup(std::set<int> &group, const std::vector<std::set
-static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getConnectedBlobs(GEntity &entity, const std::set<MElement*> &badElements, int depth, const double distFactor)
+static void calcVertex2Elements(int dim, GEntity *entity, std::map<MVertex*, std::vector<MElement *> > &vertex2elements)
{
+ for (size_t i = 0; i < entity->getNumMeshElements(); ++i) {
+ MElement *element = entity->getMeshElement(i);
+ if (element->getDim() == dim)
+ for (int j = 0; j < element->getNumPrimaryVertices(); ++j)
+ vertex2elements[element->getVertex(j)].push_back(element);
+ }
+}
- OptHomMessage("Starting blob generation from %i bad elements...",badElements.size());
- // Compute vertex -> element connectivity
- std::cout << "Computing vertex -> element connectivity...\n";
- std::map<MVertex*, std::vector<MElement *> > vertex2elements;
- for (size_t i = 0; i < entity.getNumMeshElements(); ++i) {
- MElement &element = *entity.getMeshElement(i);
- for (int j = 0; j < element.getNumPrimaryVertices(); ++j) {
- vertex2elements[element.getVertex(j)].push_back(&element);
- }
- }
+
+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)
+{
+
+ OptHomMessage("Starting blob generation from %i bad elements...",badElements.size());
// Contruct primary blobs
std::cout << "Constructing " << badElements.size() << " primary blobs...\n";
@@ -268,7 +275,7 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon
std::cout << "Computing boundaries for " << blobs.size() << " blobs...\n";
std::vector<std::pair<std::set<MElement *>, std::set<MVertex*> > > result;
for (std::list<std::set<MElement*> >::iterator itB = blobs.begin(); itB != blobs.end(); ++itB)
- result.push_back(std::pair<std::set<MElement*>, std::set<MVertex*> >(*itB, getBndVertices(*itB, vertex2elements)));
+ result.push_back(std::pair<std::set<MElement*>, std::set<MVertex*> >(*itB, getAllBndVertices(*itB, vertex2elements)));
OptHomMessage("Generated %i blobs",blobs.size());
@@ -278,90 +285,506 @@ static std::vector<std::pair<std::set<MElement*> , std::set<MVertex*> > > getCon
-void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
+static void optimizeConnectedBlobs(const std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+ std::set<MElement*> &badasses,
+ OptHomParameters &p, int samples)
{
- double t1 = Cpu();
+ p.SUCCESS = 1;
+ p.minJac = 1.e100;
+ p.maxJac = -1.e100;
+
+ std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > toOptimize =
+ getConnectedBlobs(vertex2elements, badasses, p.nbLayers, p.distanceFactor);
+
+ //#pragma omp parallel for schedule(dynamic, 1)
+ for (int i = 0; i < toOptimize.size(); ++i) {
+ OptHomMessage("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);
+//std::ostringstream ossI1;
+//ossI1 << "initial_ITER_" << i << ".msh";
+//temp.mesh.writeMSH(ossI1.str().c_str());
+ int success = temp.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) {
+ OptHomMessage("Jacobian optimization succeed, starting svd optimization");
+ success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
+ true, samples, p.itMax, p.optPassMax);
+ }
+ double minJac, maxJac, distMaxBND, distAvgBND;
+ temp.recalcJacDist();
+ temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+ p.minJac = std::min(p.minJac,minJac);
+ p.maxJac = std::max(p.maxJac,maxJac);
+ temp.mesh.updateGEntityPositions();
+ if (success <= 0) {
+ std::ostringstream ossI2;
+ ossI2 << "final_ITER_" << i << ".msh";
+ temp.mesh.writeMSH(ossI2.str().c_str());
+ }
+ //#pragma omp critical
+ p.SUCCESS = std::min(p.SUCCESS, success);
+ }
- int samples = 30;
+}
- int method = 0;
- if (p.method == 0)
- method = Mesh::METHOD_FIXBND | Mesh::METHOD_PROJJAC;
- else if (p.method == 1)
- method = Mesh::METHOD_PROJJAC;
- else if (p.method == 2)
- method = Mesh::METHOD_FIXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
- else if(p.method < 0)
- method = -p.method;
- SVector3 BND(gm->bounds().max(), gm->bounds().min());
- OptHomMessage("High order mesh optimizer starts");
- std::vector<GEntity*> entities;
- gm->getEntities(entities);
+static MElement *getWorstElement(const std::set<MElement*> &badasses, OptHomParameters &p){
+
+ double worst = 1.e300;
+ MElement *worstEl = 0;
- for (int i = 0; i < entities.size(); ++i) {
+ for (std::set<MElement*>::iterator it=badasses.begin(); it!=badasses.end(); it++) {
+ double jmin, jmax, val;
+ (*it)->scaledJacRange(jmin,jmax);
+ val = jmin-p.BARRIER_MIN + p.BARRIER_MAX-jmax;
+ if (val < worst) {
+ worst = val;
+ worstEl = *it;
+ }
+ }
- double tf1 = Cpu();
- GEntity &entity = *entities[i];
+ return worstEl;
- if (entity.dim() != p.dim || (p.onlyVisible && !entity.getVisibility())) continue;
+}
- OptHomMessage("Optimizing Model entity %4d...", entity.tag());
- std::set<MElement*> badasses;
- for (int i = 0; i < entity.getNumMeshElements();i++) {
- double jmin, jmax;
- entity.getMeshElement(i)->scaledJacRange(jmin, jmax);
- if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
- if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX)
- badasses.insert(entity.getMeshElement(i));
+
+
+static std::set<MVertex *> getPrimBndVertices(std::set<MElement*> &elements,
+ const std::map<MVertex*, std::vector<MElement*> > &vertex2elements)
+{
+ std::set<MVertex*> bnd;
+ for (std::set<MElement*>::iterator itE = elements.begin(); itE != elements.end(); ++itE) {
+ for (int i = 0; i < (*itE)->getNumPrimaryVertices(); ++i) {
+ const std::vector<MElement*> &neighbours = vertex2elements.find((*itE)->getVertex(i))->second;
+ for (size_t k = 0; k < neighbours.size(); ++k) {
+ if (elements.find(neighbours[k]) == elements.end()) {
+ bnd.insert((*itE)->getVertex(i));
+ }
+ }
+ }
+ }
+ return bnd;
+}
+
+
+
+static std::set<MElement*> getSurroundingBlob3D(MElement *el, int depth,
+ const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
+ const double distFactor)
+{
+
+ const double limDist = distMaxStraight(el)*distFactor;
+
+ std::set<MElement*> blob;
+ std::list<MElement*> currentLayer, lastLayer;
+
+ std::list<SPoint3> seedPts;
+
+ blob.insert(el);
+ lastLayer.push_back(el);
+ for (int d = 0; d < depth; ++d) {
+ currentLayer.clear();
+ for (std::list<MElement*>::iterator it = lastLayer.begin(); it != lastLayer.end(); ++it) {
+ for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
+ const std::vector<MElement*> &neighbours = vertex2elements.find((*it)->getVertex(i))->second;
+ for (std::vector<MElement*>::const_iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN) {
+ SPoint3 pt = (*itN)->barycenter_infty(); // Check distance from all seed points
+ bool nearSeed = false;
+ for (std::list<SPoint3>::const_iterator itS = seedPts.begin(); itS != seedPts.end(); ++itS)
+ if (itS->distance(pt) < limDist) {
+ nearSeed = true;
+ break;
+ }
+ if ((d == 0) || nearSeed)
+ if (blob.insert(*itN).second) currentLayer.push_back(*itN); // Assume that if an el is too far, its neighbours are too far as well
+ }
+ }
}
- if (badasses.empty()) continue;
+ if (d == 0) // Elts of 1st layer are seed points
+ for (std::list<MElement*>::iterator itL = currentLayer.begin(); itL != currentLayer.end(); ++itL)
+ seedPts.push_back((*itL)->barycenter_infty());
+ lastLayer = currentLayer;
+ }
+
+ return blob;
+
+}
- std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > toOptimize = getConnectedBlobs(entity, badasses, p.nbLayers, p.DistanceFactor);
- //#pragma omp parallel for schedule(dynamic, 1)
- p.SUCCESS = 1;
- p.minJac = 1.e100;
- p.maxJac = -1.e100;
- for (int i = 0; i < toOptimize.size(); ++i) {
- OptHomMessage("Optimizing a blob %i/%i composed of %4d elements", i+1, toOptimize.size(), toOptimize[i].first.size());
+
+static std::set<MElement*> addBlobLayer(std::set<MElement*> &blob,
+ const std::map<MVertex*, std::vector<MElement*> > &vertex2elements)
+{
+
+ std::set<MElement*> layer;
+ const std::set<MElement*> initBlob = blob;
+
+ for (std::set<MElement*>::const_iterator it = initBlob.begin(); it != initBlob.end(); ++it)
+ for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
+ const std::vector<MElement*> &neighbours = vertex2elements.find((*it)->getVertex(i))->second;
+ for (std::vector<MElement*>::const_iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN)
+ if (blob.insert(*itN).second) layer.insert(*itN);
+ }
+
+ return layer;
+
+}
+
+
+
+static bool detectNewBrokenElement(const std::set<MElement*> &layer,
+ const std::set<MElement*> &badasses, OptHomParameters &p){
+
+ for (std::set<MElement*>::iterator it=layer.begin(); it!=layer.end(); it++)
+ if (badasses.find(*it) == badasses.end()) {
+ double jmin, jmax, val;
+ (*it)->scaledJacRange(jmin,jmax);
+ if ((jmin < p.BARRIER_MIN) || (jmax > p.BARRIER_MAX)) return true;
+ }
+
+ return false;
+
+}
+
+
+
+static void optimizeOneByOne(const std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+ std::set<MElement*> badElts, OptHomParameters &p, int samples)
+{
+
+ p.SUCCESS = 1;
+ p.minJac = 1.e100;
+ p.maxJac = -1.e100;
+
+ const int initNumBadElts = badElts.size();
+ std::cout << "DBGTT: " << initNumBadElts << " badasses, starting to iterate...\n";
+
+ for (int iBadEl=0; iBadEl<initNumBadElts; iBadEl++) {
+
+ if (badElts.empty()) break;
+
+ // Create blob around worst element and remove it from badElts
+ MElement *worstEl = getWorstElement(badElts,p);
+ badElts.erase(worstEl);
+
+ int nbLayers = p.nbLayers;
+ double distanceFactor = p.distanceFactor;
+
+ int success;
+
+ for (int iterBlob=0; iterBlob<p.maxAdaptBlob; iterBlob++) {
+
+ 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()));
+ OptHomMessage("Optimizing primary blob %i (max. %i remaining) composed of %4d elements", iBadEl,
+ badElts.size(), toOptimize1.size());
fflush(stdout);
- OptHOM temp(&entity, toOptimize[i].first, toOptimize[i].second, method);
-std::ostringstream ossI1;
-ossI1 << "initial_" << entity.tag() << "ITER_" << i << ".msh";
-temp.mesh.writeMSH(ossI1.str().c_str());
- int success = temp.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) {
- OptHomMessage("Jacobian optimization succeed, starting svd optimization");
- success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
+ 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()));
+ OptHomMessage("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);
+ }
+ else {
+ OptHomMessage("Primary blob %i did not break new elements, no correction needed", iBadEl);
+ fflush(stdout);
+ }
}
- double minJac, maxJac, distMaxBND, distAvgBND;
- temp.recalcJacDist();
- temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
- p.minJac = std::min(p.minJac,minJac);
- p.maxJac = std::max(p.maxJac,maxJac);
- temp.mesh.updateGEntityPositions();
+
+ // Measure min and max Jac., update mesh
+ if ((success > 0) || (iterBlob == p.maxAdaptBlob-1)) {
+ double minJac, maxJac, distMaxBND, distAvgBND;
+ opt->recalcJacDist();
+ opt->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+ p.minJac = std::min(p.minJac,minJac);
+ p.maxJac = std::max(p.maxJac,maxJac);
+ opt->mesh.updateGEntityPositions();
+ }
+
+//std::ostringstream ossI2;
+//ossI2 << "final_ITER_" << iter << ".msh";
+//temp.mesh.writeMSH(ossI2.str().c_str());
if (success <= 0) {
- std::ostringstream ossI2;
- ossI2 << "final_" << entity.tag() << "ITER_" << i << ".msh";
- temp.mesh.writeMSH(ossI2.str().c_str());
+ distanceFactor *= p.adaptBlobDistFact;
+ nbLayers *= p.adaptBlobLayerFact;
+ OptHomMessage("Blob %i failed (adapt #%i), adapting with increased size", iBadEl, iterBlob);
+ if (iterBlob == p.maxAdaptBlob-1) {
+ std::ostringstream ossI2;
+ ossI2 << "final_" << iBadEl << ".msh";
+ opt->mesh.writeMSH(ossI2.str().c_str());
+ }
}
- //#pragma omp critical
- p.SUCCESS = std::min(p.SUCCESS, success);
+ else break;
+
}
- double DTF = Cpu()-tf1;
- if (p.SUCCESS == 1)
- OptHomMessage("Optimization succeeded for entity %d (CPU %g sec)", entity.tag(), DTF);
- else if (p.SUCCESS == 0)
- OptHomMessage("Warning : Model entity %4d has all jacobians positive but not all in the range (CPU %g sec)", entity.tag(), DTF);
- else if (p.SUCCESS == -1)
- OptHomMessage("Error : Model entity %4d has still negative jacobians (CPU %g sec)", entity.tag(), DTF);
+ //#pragma omp critical
+ p.SUCCESS = std::min(p.SUCCESS, success);
+
+ }
+
+}
+
+
+
+//static void optimizeOneByOneTest(const std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+// std::set<MElement*> badasses, OptHomParameters &p, int method, int samples)
+//{
+//
+// p.SUCCESS = 1;
+// p.minJac = 1.e100;
+// p.maxJac = -1.e100;
+//
+// const int initNumBadAsses = badasses.size();
+// std::cout << "DBGTT: " << initNumBadAsses << " badasses, starting to iterate...\n";
+//
+// for (int iter=0; iter<initNumBadAsses; iter++) {
+//
+// if (badasses.empty()) break;
+//
+// // Create blob around worst element
+// MElement *worstEl = getWorstElement(badasses,p);
+//
+// const int maxRepeatBlob = 1;
+// const int repeatBlobLayerFact = 2;
+// const int repeatBlobDistFact = 2.;
+//
+// int nbLayers = p.nbLayers, minNbLayers = 1;
+// double distanceFactor = p.distanceFactor;
+//
+// int success;
+//
+// for (int iterBlob=0; iterBlob<maxRepeatBlob; iterBlob++) {
+//
+// // TODO: First opt. with only 1st-order bnd. vertices fixed,
+// // then add external layer and opt. with all bnd. vert. fixed
+//// std::set<MElement*> toOptimizePrim = getSurroundingBlob(worstEl, nbLayers, vertex2elements, distanceFactor, 0);
+//// std::set<MVertex*> toFix = getBndVertices3D(toOptimizePrim, vertex2elements);
+// std::set<MElement*> toOptimizePrim = getSurroundingBlob3D(worstEl, nbLayers, vertex2elements, distanceFactor);
+// std::set<MVertex*> toFix = getAllBndVertices(toOptimizePrim, vertex2elements);
+// badasses.erase(worstEl);
+// std::set<MElement*> toOptimize;
+// std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),badasses.begin(),badasses.end(),
+// std::inserter(toOptimize, toOptimize.end()));
+//
+// // Optimize blob
+// OptHomMessage("Optimizing blob %i (max. %i remaining) composed of %4d elements", iter,
+// badasses.size(), toOptimize.size());
+// fflush(stdout);
+// OptHOM temp(toOptimize, toFix, method);
+////std::cout << "DBGTT: toOptimize:\n";
+////for (std::set<MElement*>::iterator it=toOptimize.begin(); it!=toOptimize.end(); it++ ) {
+//// SPoint3 bar = (*it)->barycenter();
+//// std::cout << "DBGTT: -> (" << bar.x() << "," << bar.y() << "," << bar.z() << ")\n";
+//// std::cout << "DBGTT: -> (" << (*it)->getVertex(0)->onWhat()->dim();
+//// for (int j=1; j<(*it)->getNumVertices(); j++) std::cout << "," << (*it)->getVertex(j)->onWhat()->dim();
+//// std::cout << ")\n";
+////}
+////std::cout << "DBGTT: toFix:\n";
+////for (std::set<MVertex*>::iterator it=toFix.begin(); it!=toFix.end(); it++ )
+//// std::cout << "DBGTT: -> (" << (*it)->x() << "," << (*it)->y() << "," << (*it)->z() << ")\n";
+//std::ostringstream ossI1;
+//ossI1 << "initial_ITER_" << iter << ".msh";
+//temp.mesh.writeMSH(ossI1.str().c_str());
+// success = temp.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) {
+// OptHomMessage("Jacobian optimization succeed, starting svd optimization");
+// success = temp.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
+// double minJac, maxJac, distMaxBND, distAvgBND;
+// temp.recalcJacDist();
+// temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+// p.minJac = std::min(p.minJac,minJac);
+// p.maxJac = std::max(p.maxJac,maxJac);
+// temp.mesh.updateGEntityPositions();
+//
+////std::ostringstream ossI2;
+////ossI2 << "final_ITER_" << iter << ".msh";
+////temp.mesh.writeMSH(ossI2.str().c_str());
+// if (success <= 0) {
+// distanceFactor *= repeatBlobDistFact;
+// nbLayers *= repeatBlobLayerFact;
+// OptHomMessage("Blob %i failed (repeat %i), repeating with increased size", iter, iterBlob);
+// if (iterBlob == maxRepeatBlob-1) {
+// std::ostringstream ossI2;
+// ossI2 << "final_ITER_" << iter << ".msh";
+// temp.mesh.writeMSH(ossI2.str().c_str());
+// }
+// }
+// else break;
+//
+// }
+//
+// //#pragma omp critical
+// p.SUCCESS = std::min(p.SUCCESS, success);
+//
+// }
+//
+//}
+//
+//
+//
+//static void optimizePropagate(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
+// std::set<MElement*> badasses, OptHomParameters &p, int method, int samples)
+//{
+//
+// p.SUCCESS = 1;
+// p.minJac = 1.e100;
+// p.maxJac = -1.e100;
+//
+// const int initNumBadAsses = badasses.size();
+// std::cout << "DBGTT: " << initNumBadAsses << " badasses, starting to iterate...\n";
+//
+// std::set<MVertex*> toFix;
+// std::set<MElement*> done;
+//
+// while (!badasses.empty()) {
+//
+// OptHomMessage("Optimizing first of %i bad elements, %i already done", badasses.size(), done.size());
+// fflush(stdout);
+////char dum;
+////std::cin >> dum;
+//
+// MElement *el = *badasses.begin();
+//
+// // Optimize element
+// std::set<MElement*> blob;
+// blob.insert(el);
+// OptHOM temp(blob, toFix, method);
+////std::ostringstream ossI1;
+////ossI1 << "initial.msh";
+////temp.mesh.writeMSH(ossI1.str().c_str());
+// int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
+// false, samples, p.itMax, p.optPassMax);
+// if (success <= 0) {
+// OptHomMessage("Error, element optimization failed");
+// fflush(stdout);
+// std::ostringstream ossI1;
+// ossI1 << "failed_" << el->getNum() << ".msh";
+// temp.mesh.writeMSH(ossI1.str().c_str());
+// }
+////std::ostringstream ossI2;
+////ossI2 << "final.msh";
+////temp.mesh.writeMSH(ossI2.str().c_str());
+//
+// // Move element from list of bad elements to liste of elements done
+// badasses.erase(el);
+// done.insert(el);
+//std::cout << "DBGTT: Removing el. " << el->getNum() << "\n";
+//
+// // Measure success, min and max Jac and update mesh
+// double minJac, maxJac, distMaxBND, distAvgBND;
+// temp.recalcJacDist();
+// temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+// p.minJac = std::min(p.minJac,minJac);
+// p.maxJac = std::max(p.maxJac,maxJac);
+// p.SUCCESS = std::min(p.SUCCESS, success);
+// temp.mesh.updateGEntityPositions();
+//
+// // Fix elements nodes
+// for (int i = 0; i < el->getNumVertices(); ++i) toFix.insert(el->getVertex(i));
+//
+// // Add elements that have been broken by optimization to list of bad elements
+// std::set<MElement*> layer;
+// for (int i = 0; i < el->getNumPrimaryVertices(); ++i) {
+// const std::vector<MElement*> &neighbours = vertex2elements.find(el->getVertex(i))->second;
+// for (size_t k = 0; k < neighbours.size(); ++k) layer.insert(neighbours[k]);
+// }
+// for (std::set<MElement*>::const_iterator it = layer.begin(); it != layer.end(); ++it)
+// if (done.find(*it) == done.end()) { // Add only if not done
+// double jmin, jmax;
+// (*it)->scaledJacRange(jmin, jmax);
+// if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
+// if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) {
+//std::cout << "DBGTT: Adding el. " << (*it)->getNum() << "\n";
+// badasses.insert(*it);
+// }
+// }
+//
+// }
+//
+//}
+
+
+
+void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
+{
+
+ double t1 = Cpu();
+
+ int samples = 30;
+
+ double tf1 = Cpu();
+
+ OptHomMessage("High order mesh optimizer starts");
+ std::vector<GEntity*> entities;
+ gm->getEntities(entities);
+ std::map<MVertex*, std::vector<MElement *> > vertex2elements;
+ 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;
+ std::cout << "Computing connectivity and bad elements for entity " << entity->tag() << " ...\n";
+ calcVertex2Elements(p.dim,entity,vertex2elements); // Compute vert. -> elt. connectivity
+//std::cout << " -> " << entity->getNumMeshElements()
+// << " elements, vertex2elements.size() = " << vertex2elements.size() << "\n";
+ 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);
+ if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
+ if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) badasses.insert(el);
+ }
+ }
}
+ if (p.strategy == 0) optimizeConnectedBlobs(vertex2elements, badasses, p, samples);
+ else optimizeOneByOne(vertex2elements, badasses, p, samples);
+
+ double DTF = Cpu()-tf1;
+ if (p.SUCCESS == 1)
+ OptHomMessage("Optimization succeeded (CPU %g sec)", DTF);
+ else if (p.SUCCESS == 0)
+ OptHomMessage("Warning : All jacobians positive but not all in the range (CPU %g sec)", DTF);
+ else if (p.SUCCESS == -1)
+ OptHomMessage("Error : Still negative jacobians (CPU %g sec)", DTF);
+
double t2 = Cpu();
p.CPU = t2-t1;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index 1e1c3a082c12031bd23014888a5a1c081dcb39e5..2382923312a25b0ddce40c373af01dc368c93057 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -16,9 +16,14 @@ struct OptHomParameters {
int optPassMax ; // max number of optimization passes
double TMAX ; // max CPU time allowed
bool onlyVisible ; // apply optimization to visible entities ONLY
- double DistanceFactor; // filter elements such that no elements further away than
+ 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
- int method ; // how jacobians are computed and if points can move on boundaries
+ bool fixBndNodes; // how jacobians are computed and if points can move on boundaries
+ int strategy; // 0 = connected blobs, 1 = adaptive one-by-one
+ 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
+
// OUTPUT ------>
int SUCCESS ; // 0 --> success , 1 --> Not converged
double minJac, maxJac; // after optimization, range of jacobians
@@ -29,8 +34,9 @@ struct OptHomParameters {
OptHomParameters ()
// default values
- : BARRIER_MIN_METRIC (-1.), BARRIER_MIN (0.1), BARRIER_MAX (2.0) , weightFixed (1.e6), weightFree (1.e2),
- nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true), DistanceFactor(12), method(1)
+ : 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.)
{
}
};
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
index cdb3e3c8db07d182b61f63b46123bd7cd430a5d0..b8f2518e194f6418097daa71a99ccd85b03540aa 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
@@ -6,53 +6,6 @@
-ParamCoordSurf::ParamCoordSurf(GEntity *ge)
-{
- if ((ge->dim() == 2) && (ge->geomType() != GEntity::DiscreteSurface)) _gf = static_cast<GFace*>(ge);
- else std::cout << "ERROR: Using surface parametric coordinates with non-surface geometric entity" << std::endl;
-}
-
-
-
-SPoint3 ParamCoordSurf::getUvw(MVertex* vert)
-{
- SPoint2 p;
- reparamMeshVertexOnFace(vert,_gf,p);
- return SPoint3(p[0],p[1],0.);
-}
-
-
-
-SPoint3 ParamCoordSurf::uvw2Xyz(MVertex* vert, const SPoint3 &uvw)
-{
- GPoint gp = _gf->point(uvw[0],uvw[1]);
- return SPoint3(gp.x(),gp.y(),gp.z());
-}
-
-
-
-void ParamCoordSurf::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw)
-{
- Pair<SVector3,SVector3> der = _gf->firstDer(SPoint2(uvw[0],uvw[1]));
- gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() + gXyz.z() * der.first().z();
- gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() + gXyz.z() * der.second().z();
-}
-
-
-
-void ParamCoordSurf::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
-{
- Pair<SVector3,SVector3> der = _gf->firstDer(SPoint2(uvw[0],uvw[1]));
- std::vector<SPoint3>::iterator itUvw=gUvw.begin();
- for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
- (*itUvw)[0] = itXyz->x() * der.first().x() + itXyz->y() * der.first().y() + itXyz->z() * der.first().z();
- (*itUvw)[1] = itXyz->x() * der.second().x() + itXyz->y() * der.second().y() + itXyz->z() * der.second().z();
- itUvw++;
- }
-}
-
-
-
SPoint3 ParamCoordParent::getUvw(MVertex* vert)
{
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.h b/contrib/HighOrderMeshOptimizer/ParamCoord.h
index cb3db207db7cd24f2381a02e607b33c24b25db7b..1927505508f26b6f4ec9729176314ebbfd46cc91 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.h
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.h
@@ -64,26 +64,6 @@ public:
-class ParamCoordSurf : public ParamCoord
-{
-
-public:
-
- ParamCoordSurf(GEntity *ge);
- int nCoord(MVertex* vert) { return 2; }
- virtual SPoint3 getUvw(MVertex* vert);
- virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw);
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw);
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw);
-
-private:
-
- GFace *_gf;
-
-};
-
-
-
class ParamCoordParent : public ParamCoord
{