From 577221ee8f42ce19c13ed1355366924f0cd73c57 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Mon, 28 Oct 2013 23:09:31 +0000
Subject: [PATCH] Added maxDistToStraight to MElement. High order mesh
optimization: fixed bug with length scale, worked on one-by-one strategy.
---
Fltk/highOrderToolsWindow.cpp | 4 +-
Geo/MElement.cpp | 23 ++
Geo/MElement.h | 3 +
contrib/HighOrderMeshOptimizer/OptHOM.cpp | 2 +-
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp | 27 +--
contrib/HighOrderMeshOptimizer/OptHomRun.cpp | 225 ++++++++----------
contrib/HighOrderMeshOptimizer/OptHomRun.h | 4 +-
7 files changed, 131 insertions(+), 157 deletions(-)
diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index d9410ed48e..93683f2427 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -343,11 +343,11 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
value[5] = new Fl_Value_Input
(x, y, IW/2, BH);
value[5]->align(FL_ALIGN_RIGHT);
- value[5]->value(1.e+5);
+ value[5]->value(1000.);
value[6] = new Fl_Value_Input
(x+IW/2,y, IW/2, BH, "W fixed / W free");
value[6]->align(FL_ALIGN_RIGHT);
- value[6]->value(1.e+2);
+ value[6]->value(1.);
y += BH;
value[3] = new Fl_Value_Input
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6799ec9e7e..7a4acf9cb8 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -112,6 +112,29 @@ double MElement::rhoShapeMeasure()
return 0.;
}
+double MElement::maxDistToStraight()
+{
+ const nodalBasis *lagBasis = getFunctionSpace();
+ const fullMatrix<double> &uvw = lagBasis->points;
+ const int &nV = uvw.size1();
+ const int &dim = uvw.size2();
+ const nodalBasis *lagBasis1 = getFunctionSpace(1);
+ const int &nV1 = lagBasis1->points.size1();
+ SPoint3 xyz1[nV1];
+ for (int iV = 0; iV < nV1; ++iV) xyz1[iV] = getVertex(iV)->point();
+ double maxdx = 0.;
+ for (int iV = nV1; iV < nV; ++iV) {
+ double f[256];
+ lagBasis1->f(uvw(iV, 0), (dim > 1) ? uvw(iV, 1) : 0., (dim > 2) ? uvw(iV, 2) : 0., f);
+ SPoint3 xyzS(0.,0.,0.);
+ for (int iSF = 0; iSF < nV1; ++iSF) xyzS += xyz1[iSF]*f[iSF];
+ SVector3 vec(xyzS,getVertex(iV)->point());
+ double dx = vec.norm();
+ if (dx > maxdx) maxdx = dx;
+ }
+ return maxdx;
+}
+
void MElement::scaledJacRange(double &jmin, double &jmax, GEntity *ge)
{
jmin = jmax = 1.0;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 67c280f409..12e362dc63 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -175,6 +175,9 @@ class MElement
virtual double maxEdge();
virtual double minEdge();
+ // Max. distance between curved and straight element among all high-order nodes
+ double maxDistToStraight();
+
// get the quality measures
virtual double rhoShapeMeasure();
virtual double gammaShapeMeasure(){ return 0.; }
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index c7fb2a88a3..6ecfb3bf9c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -316,7 +316,7 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min,
// Set weights & length scale for non-dimensionalization
lambda = weightFixed;
lambda2 = weightFree;
- std::vector<double> dSq(mesh.nVert());
+ std::vector<double> dSq(mesh.nEl());
mesh.distSqToStraight(dSq);
const double maxDSq = *max_element(dSq.begin(),dSq.end());
if (maxDSq < 1.e-10) { // Length scale for non-dim. distance
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index e24731670d..e8cdf9de65 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -207,31 +207,10 @@ void Mesh::updateMesh(const double *it)
}
-void Mesh::distSqToStraight(std::vector<double> &dSq)
-{
- std::vector<SPoint3> sxyz(nVert());
+void Mesh::distSqToStraight(std::vector<double> &dSq) {
for (int iEl = 0; iEl < nEl(); iEl++) {
- MElement *el = _el[iEl];
- const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
- const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
- int nV = lagrange->points.size1();
- int nV1 = lagrange1->points.size1();
- for (int i = 0; i < nV1; ++i) {
- sxyz[_el2V[iEl][i]] = _vert[_el2V[iEl][i]]->point();
- }
- int dim = lagrange->points.size2();
- for (int i = nV1; i < nV; ++i) {
- double f[256];
- lagrange1->f(lagrange->points(i, 0), dim > 1 ? lagrange->points(i, 1) : 0.,
- dim > 2 ? lagrange->points(i, 2) : 0., f);
- for (int j = 0; j < nV1; ++j)
- sxyz[_el2V[iEl][i]] += sxyz[_el2V[iEl][j]] * f[j];
- }
- }
-
- for (int iV = 0; iV < nVert(); iV++) {
- SPoint3 d = _xyz[iV]-sxyz[iV];
- dSq[iV] = d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
+ const double d = _el[iEl]->maxDistToStraight();
+ dSq[iEl] = d*d;
}
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 7a013dec96..b5f5d615ec 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -47,37 +47,6 @@
#if defined(HAVE_BFGS)
-double distMaxStraight(MElement *el)
-{
- const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
- const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
- int nV = lagrange->points.size1();
- int nV1 = lagrange1->points.size1();
- SPoint3 sxyz[256];
- for (int i = 0; i < nV1; ++i) {
- sxyz[i] = el->getVertex(i)->point();
- }
- for (int i = nV1; i < nV; ++i) {
- double f[256];
- double u = 0., v = 0., w = 0.;
- if(lagrange->points.size2() > 0) u = lagrange->points(i, 0);
- if(lagrange->points.size2() > 1) v = lagrange->points(i, 1);
- if(lagrange->points.size2() > 2) w = lagrange->points(i, 2);
- lagrange1->f(u, v, w, f);
- for (int j = 0; j < nV1; ++j)
- sxyz[i] += sxyz[j] * f[j];
- }
-
- double maxdx = 0.0;
- for (int iV = nV1; iV < nV; iV++) {
- SVector3 d = el->getVertex(iV)->point()-sxyz[iV];
- double dx = d.norm();
- if (dx > maxdx) maxdx = dx;
- }
-
- return maxdx;
-}
-
void exportMeshToDassault(GModel *gm, const std::string &fn, int dim)
{
FILE *f = fopen(fn.c_str(),"w");
@@ -167,7 +136,7 @@ static std::set<MElement*> getSurroundingBlob
{
const SPoint3 p = el->barycenter();
- const double dist = distMaxStraight(el);
+ const double dist = el->maxDistToStraight();
const double limDist = ((optPrimSurfMesh && (dist < 1.e-10)) ?
el->getOuterRadius() : dist) * distFactor;
@@ -395,7 +364,7 @@ static std::set<MElement*> getSurroundingBlob3D
const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
const double distFactor)
{
- const double limDist = distMaxStraight(el) * distFactor;
+ const double limDist = el->maxDistToStraight() * distFactor;
std::set<MElement*> blob;
std::list<MElement*> currentLayer, lastLayer;
@@ -499,111 +468,111 @@ 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, p.optPrimSurfMesh);
+// 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(element2entity, 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(element2entity, 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, p.optPrimSurfMesh);
- 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(element2entity, 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);
+ 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);
+ }
- // Second step: add external layer, check it and optimize it safely (all
- // bnd. vertices fixed) if new broken element
- if (success > 0) {
+ // 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::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(element2entity, 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);
- }
}
-// std::set<MElement*> toOptimizePrim = getSurroundingBlob
-// (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(),
-// std::inserter(toOptimize, toOptimize.end()));
-// Msg::Info("Optimizing blob %i (max. %i remaining) "
-// "composed of %4d elements", iBadEl, badElts.size(),
-// toOptimize.size());
-// fflush(stdout);
-// OptHOM *opt = new OptHOM(element2entity, toOptimize, toFix, 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);
-// }
-//
-// // 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) {
-// distanceFactor *= p.adaptBlobDistFact;
-// nbLayers *= p.adaptBlobLayerFact;
-// Msg::Info("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());
-//// }
-// }
-// else break;
+ //std::ostringstream ossI2;
+ //ossI2 << "final_ITER_" << iter << ".msh";
+ //temp.mesh.writeMSH(ossI2.str().c_str());
+ if (success <= 0) {
+ distanceFactor *= p.adaptBlobDistFact;
+ nbLayers *= p.adaptBlobLayerFact;
+ Msg::Info("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());
+// }
+ }
+ else break;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index f82751da05..e934f01ad3 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -61,8 +61,8 @@ struct OptHomParameters {
double CPU; // Time for optimization
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),
+ : BARRIER_MIN_METRIC(-1.), BARRIER_MIN(0.1), BARRIER_MAX(2.0), weightFixed(1000.),
+ weightFree (1.), nbLayers (6) , dim(3) , itMax(300), onlyVisible(true),
distanceFactor(12), fixBndNodes(false), strategy(0), maxAdaptBlob(3),
adaptBlobLayerFact(2.), adaptBlobDistFact(2.), optPrimSurfMesh(false)
{
--
GitLab