Skip to content
Snippets Groups Projects
Commit 577221ee authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Added maxDistToStraight to MElement. High order mesh optimization: fixed bug... 

Added maxDistToStraight to MElement. High order mesh optimization: fixed bug with length scale, worked on one-by-one strategy.
parent 898108b7
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Fltk/highOrderToolsWindow.cpp
+ 2
− 2
View file @ 577221ee
......@@ -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
......
Geo/MElement.cpp
+ 23
− 0
View file @ 577221ee
......@@ -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;
......
Geo/MElement.h
+ 3
− 0
View file @ 577221ee
......@@ -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.; }
......
contrib/HighOrderMeshOptimizer/OptHOM.cpp
+ 1
− 1
View file @ 577221ee
......@@ -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
......
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+ 3
− 24
View file @ 577221ee
......@@ -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;
}
}
......
contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+ 97
− 128
View file @ 577221ee
......@@ -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,44 +468,75 @@ 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;
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;
(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(toOptimize2, toOptimize2.end()));
Msg::Info("Optimizing corrective blob %i (max. %i remaining) "
std::inserter(toOptimize, toOptimize.end()));
Msg::Info("Optimizing blob %i (max. %i remaining) "
"composed of %4d elements", iBadEl, badElts.size(),
toOptimize2.size());
toOptimize.size());
fflush(stdout);
opt = new OptHOM(element2entity, toOptimize2, toFix2, p.fixBndNodes);
OptHOM *opt = new OptHOM(element2entity, toOptimize, toFix, p.fixBndNodes);
//std::ostringstream ossI1;
//ossI1 << "initial_corrective_" << iter << ".msh";
//opt->mesh.writeMSH(ossI1.str().c_str());
......@@ -547,63 +547,32 @@ static void optimizeOneByOne
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);
}
// 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*> 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());
//// }
//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;
}
else break;
}
......
contrib/HighOrderMeshOptimizer/OptHomRun.h
+ 2
− 2
View file @ 577221ee
......@@ -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)
{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment