From af362f914a82d50ff5f64df9ce6f8f455b3ef1d7 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 24 Apr 2012 14:04:19 +0000
Subject: [PATCH] high order tools can now deal with BL meshes !!!
---
Fltk/highOrderToolsWindow.cpp | 152 ++++--
Geo/GEdge.cpp | 48 +-
Geo/GEdge.h | 2 +
Geo/MElement.cpp | 21 +
Geo/MElement.h | 2 +
Mesh/HighOrder.cpp | 13 +
Mesh/HighOrder.h | 6 +
Mesh/meshGFaceDelaunayInsertion.cpp | 4 +-
contrib/HighOrderMeshOptimizer/OptHOM.cpp | 214 ++++-----
contrib/HighOrderMeshOptimizer/OptHOM.h | 21 +-
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp | 84 ++--
contrib/HighOrderMeshOptimizer/OptHomMesh.h | 2 +-
contrib/HighOrderMeshOptimizer/OptHomRun.cpp | 443 ++++++++++++++++--
contrib/HighOrderMeshOptimizer/OptHomRun.h | 12 +-
14 files changed, 751 insertions(+), 273 deletions(-)
diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index de8d0eac70..843508804f 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -53,13 +53,69 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
bool onlyVisible = (bool)o->butt[1]->value();
int nLayers = (int) o->value[2]->value();
- SetOrderN(GModel::current(), order, linear, incomplete, onlyVisible);
+ if (order == 1)
+ SetOrder1(GModel::current());
+ else
+ SetOrderN(GModel::current(), order, linear, incomplete, onlyVisible);
+
+ distanceFromMeshToGeometry_t dist;
+ computeDistanceFromMeshToGeometry (GModel::current(), dist);
+ for (std::map<GEntity*, double> ::iterator it = dist.d2.begin(); it !=dist.d2.end();++it){
+ printf ("GEntity %d of dim %d : dist %12.5E\n",it->first->tag(),it->first->dim(),it->second);
+ }
+
CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
drawContext::global()->draw();
}
+static void chooseopti_cb(Fl_Widget *w, void *data){
+ highOrderToolsWindow *o = FlGui::instance()->highordertools;
+ bool elastic = (bool)o->butt[3]->value();
+
+ if (elastic){
+ o->butt[4]->value(0);
+ o->choice[0]->hide();
+ for (int i=3;i<=6;i++)
+ o->value[i]->hide();
+ o->push[1]->hide();
+ }
+ else {
+ o->butt[3]->value(0);
+ o->push[0]->show();
+ o->choice[0]->show();
+ for (int i=3;i<=6;i++)
+ o->value[i]->show();
+ o->push[1]->show();
+ }
+
+}
+
+
+static void chooseopti2_cb(Fl_Widget *w, void *data){
+ highOrderToolsWindow *o = FlGui::instance()->highordertools;
+ bool elastic = !(bool)o->butt[4]->value();
+
+ if (elastic){
+ o->butt[4]->value(0);
+ o->choice[0]->hide();
+ for (int i=3;i<=6;i++)
+ o->value[i]->hide();
+ o->push[1]->hide();
+ }
+ else {
+ o->butt[3]->value(0);
+ o->push[0]->show();
+ o->choice[0]->show();
+ for (int i=3;i<=6;i++)
+ o->value[i]->show();
+ o->push[1]->show();
+ }
+
+}
+
+
static void highordertools_runelas_cb(Fl_Widget *w, void *data)
{
highOrderToolsWindow *o = FlGui::instance()->highordertools;
@@ -68,21 +124,29 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data)
bool linear = !(bool)o->butt[2]->value();
bool incomplete = (bool)o->butt[0]->value();
bool elastic = (bool)o->butt[3]->value();
- double threshold = o->value[1]->value();
+ double threshold_min = o->value[1]->value();
bool onlyVisible = (bool)o->butt[1]->value();
int nbLayers = (int) o->value[2]->value();
+ double threshold_max = o->value[8]->value();
- if(elastic)ElasticAnalogy(GModel::current(), threshold, onlyVisible);
+ if(elastic)ElasticAnalogy(GModel::current(), threshold_min, onlyVisible);
else {
OptHomParameters p;
p.nbLayers = nbLayers;
- p.BARRIER = threshold;
+ p.BARRIER_MIN = threshold_min;
+ p.BARRIER_MAX = threshold_max;
p.onlyVisible = onlyVisible;
p.dim = GModel::current()->getNumRegions() ? 3 : 2;
p.itMax = (int) o->value[3]->value();
+ p.weightFixed = o->value[5]->value();
+ p.weightFree = o->value[6]->value();
+ p.DistanceFactor = o->value[7]->value();
+ p.method = (int)o->choice[0]->value();
+ p.filter = (int)o->choice[1]->value();
HighOrderMeshOptimizer (GModel::current(),p);
}
+
CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
drawContext::global()->draw();
@@ -160,9 +224,8 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
butt[2]->type(FL_TOGGLE_BUTTON);
butt[2]->value(1);
- y += 1.2*BH;
push[0] = new Fl_Button
- (x, y, (int)(1.2*BB), BH, "Apply");
+ (x+1.9*IW, y, (int)(1.2*BB), BH, "Apply");
push[0]->callback(highordertools_runp_cb);
@@ -191,12 +254,20 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
value[1] = new Fl_Value_Input
- (x,y, IW, BH, "Distorsion threshold");
+ (x,y, IW/2.0, BH);
value[1]->minimum(0);
value[1]->maximum(1);
value[1]->step(.01);
value[1]->align(FL_ALIGN_RIGHT);
value[1]->value(0.1);
+
+ value[8] = new Fl_Value_Input
+ (x+IW/2.0,y, IW/2.0, BH, "Jacobian Range");
+ value[8]->minimum(1);
+ value[8]->maximum(10);
+ value[8]->step(.01);
+ value[8]->align(FL_ALIGN_RIGHT);
+ value[8]->value(2);
y += 1.2*BH;
value[2] = new Fl_Value_Input
@@ -205,7 +276,17 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
value[2]->maximum(20);
value[2]->step(1);
value[2]->align(FL_ALIGN_RIGHT);
- value[2]->value(3);
+ value[2]->value(6);
+
+ y += 1.2*BH;
+ value[7] = new Fl_Value_Input
+ (x,y, IW, BH, "Distance Factor");
+ value[7]->minimum(1);
+ value[7]->maximum(20000);
+ value[7]->step(1);
+ value[7]->align(FL_ALIGN_RIGHT);
+ value[7]->value(12);
+
}
{
y += BH;
@@ -213,6 +294,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
(x,y, IW, BH, "Relocate mesh vertices through Elastic Analogy");
butt[3]->type(FL_TOGGLE_BUTTON);
butt[3]->value(0);
+ butt[3]->callback(chooseopti_cb);
}
{
@@ -220,32 +302,47 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
butt[4] = new Fl_Check_Button
(x,y, IW, BH, "Use Optimization for relocating mesh vertices");
butt[4]->type(FL_TOGGLE_BUTTON);
- butt[4]->value(0);
+ butt[4]->value(1);
+ butt[4]->callback(chooseopti2_cb);
- static Fl_Menu_Item menu_opti[] = {
- {"IpOpt", 0, 0, 0},
- {"LBFGS", 0, 0, 0},
- {"Conjugate Gradients", 0, 0, 0},
+ static Fl_Menu_Item menu_objf[] = {
+ {"CAD + FIXBND", 0, 0, 0},
+ {"CAD + FREEBND", 0, 0, 0},
+ {"MESH ONLY", 0, 0, 0},
{0}
};
- static Fl_Menu_Item menu_objective[] = {
- {"Log Barriers", 0, 0, 0},
- {"Polynomial", 0, 0, 0},
+
+ static Fl_Menu_Item menu_filter[] = {
+ {"GLOBAL", 0, 0, 0},
+ {"BLOBS ", 0, 0, 0},
{0}
};
-
+
y += BH;
choice[0] = new Fl_Choice
- (x,y, IW, BH, "Optimization Package");
- choice[0]->menu(menu_opti);
+ (x,y, IW, BH, "Objective Function");
+ choice[0]->menu(menu_objf);
choice[0]->align(FL_ALIGN_RIGHT);
+
y += BH;
- choice[0] = new Fl_Choice
- (x,y, IW, BH, "Objective Function");
- choice[0]->menu(menu_objective);
- choice[0]->align(FL_ALIGN_RIGHT);
+ choice[1] = new Fl_Choice
+ (x,y, IW, BH, "Filter");
+ choice[1]->menu(menu_filter);
+ choice[1]->align(FL_ALIGN_RIGHT);
+
+
+ y += 1.2*BH;
+ value[5] = new Fl_Value_Input
+ (x,y, IW, BH, "W fixed");
+ value[5]->align(FL_ALIGN_RIGHT);
+ value[5]->value(1.e+5);
+
+ value[6] = new Fl_Value_Input
+ (x+IW*1.5,y, IW, BH, "W free");
+ value[6]->align(FL_ALIGN_RIGHT);
+ value[6]->value(1.e+2);
y += 1.2*BH;
value[3] = new Fl_Value_Input
@@ -256,6 +353,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
value[3]->align(FL_ALIGN_RIGHT);
value[3]->value(300);
+
y += 1.2*BH;
value[4] = new Fl_Value_Input
(x,y, IW, BH, "Tolerance");
@@ -265,11 +363,9 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
value[4]->align(FL_ALIGN_RIGHT);
value[4]->value(1.e-4);
- y += 1.2*BH;
- push[0] = new Fl_Button
- (x, y, (int)(1.2*BB), BH, "Apply");
- push[0]->callback(highordertools_runelas_cb);
-
+ push[1] = new Fl_Button
+ (x+1.9*IW, y, (int)(1.2*BB), BH, "Apply");
+ push[1]->callback(highordertools_runelas_cb);
}
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 7e4548598c..e2a4b40e6b 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -343,11 +343,11 @@ GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
// printf("parameter %g as an initial guess (dist = %g)\n",topt,DMIN);
if (topt == tMin)
- t = goldenSectionSearch (this, q, topt, topt + DT/2, topt + DT, 1.e-7);
+ t = goldenSectionSearch (this, q, topt, topt + DT/2, topt + DT, 1.e-9);
else if (topt == tMax)
- t = goldenSectionSearch (this, q, topt - DT, topt - DT/2 , topt, 1.e-7);
+ t = goldenSectionSearch (this, q, topt - DT, topt - DT/2 , topt, 1.e-9);
else
- t = goldenSectionSearch (this, q, topt - DT, topt, topt + DT, 1.e-7);
+ t = goldenSectionSearch (this, q, topt - DT, topt, topt + DT, 1.e-9);
const SVector3 dp = q - position(t);
// const double D = dp.norm();
@@ -356,6 +356,48 @@ GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
return point(t);
}
+bool GEdge::computeDistanceFromMeshToGeometry (double &d2, double &dmax) {
+ d2 = 0.0; dmax = 0.0;
+ if (geomType() == Line) return true;
+ if (!lines.size())return false;
+ IntPt *pts;
+ int npts;
+ lines[0]->getIntegrationPoints(2*lines[0]->getPolynomialOrder(), &npts, &pts);
+
+ for (int i=0;i<lines.size();i++){
+ MLine *l = lines[i];
+ double t[256];
+
+ for (int j=0; j< l->getNumVertices();j++){
+ MVertex *v = l->getVertex(j);
+ if (v->onWhat() == getBeginVertex()){
+ t[j] = getLowerBound();
+ }
+ else if (v->onWhat() == getEndVertex()){
+ t[j] = getUpperBound();
+ }
+ else {
+ v->getParameter(0,t[j]);
+ }
+ }
+ for (int j=0;j<npts;j++){
+ SPoint3 p;
+ l->pnt(pts[j].pt[0],0,0,p);
+ double tinit = l->interpolate(t,pts[j].pt[0],0,0);
+ GPoint pc = closestPoint(p, tinit);
+ if (!pc.succeeded())continue;
+ double dsq =
+ (pc.x()-p.x())*(pc.x()-p.x()) +
+ (pc.y()-p.y())*(pc.y()-p.y()) +
+ (pc.z()-p.z())*(pc.z()-p.z());
+ d2 += pts[i].weight * fabs(l->getJacobianDeterminant(pts[j].pt[0],0,0)) * dsq;
+ dmax = std::max(dmax,sqrt(dsq));
+ }
+ }
+ d2 = sqrt(d2);
+ return true;
+}
+
double GEdge::parFromPoint(const SPoint3 &P) const
{
double t;
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 55fe4a6dca..e86c44c509 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -199,6 +199,8 @@ class GEdge : public GEntity {
std::vector<MLine*> lines;
void addLine(MLine *line){ lines.push_back(line); }
+
+ bool computeDistanceFromMeshToGeometry (double &d2, double &dmax);
};
#endif
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 7dff4ec5f2..f0af7c17b8 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -149,6 +149,27 @@ void MElement::getThirdDerivativeShapeFunctions(double u, double v, double w, do
else Msg::Error("Function space not implemented for this type of element");
};
+SPoint3 MElement::barycenter_infty ()
+{
+ double xmin = getVertex(0)->x();
+ double xmax = xmin;
+ double ymin = getVertex(0)->y();
+ double ymax = ymin;
+ double zmin = getVertex(0)->z();
+ double zmax = zmin;
+ int n = getNumVertices();
+ for(int i = 0; i < n; i++) {
+ MVertex *v = getVertex(i);
+ xmin = std::min(xmin,v->x());
+ xmax = std::max(xmax,v->x());
+ ymin = std::min(ymin,v->y());
+ ymax = std::max(ymax,v->y());
+ zmin = std::min(zmin,v->z());
+ zmax = std::max(zmax,v->z());
+ }
+ return SPoint3(0.5*(xmin+xmax),0.5*(ymin+ymax),0.5*(zmin+zmax));
+}
+
SPoint3 MElement::barycenter()
{
SPoint3 p(0., 0., 0.);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index f4dc4dfde8..65c48e5a36 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -195,6 +195,8 @@ class MElement
// compute the barycenter
virtual SPoint3 barycenter();
virtual SPoint3 barycenterUVW();
+ // compute the barycenter in infinity norm
+ virtual SPoint3 barycenter_infty();
// revert the orientation of the element
virtual void revert(){}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 26fb4a802b..fc18842764 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1361,3 +1361,16 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
Msg::StatusBar(2, true, "Done meshing order %d (%g s)", order, t2 - t1);
}
+
+void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t &dist) {
+ for (GModel::eiter itEdge = m->firstEdge(); itEdge != m->lastEdge(); ++itEdge) {
+ double d2,dmax;
+ (*itEdge)->computeDistanceFromMeshToGeometry (d2,dmax);
+ dist.d2[*itEdge] = d2;
+ dist.d_max[*itEdge] = dmax;
+ }
+
+ for (GModel::fiter itFace = m->firstFace(); itFace != m->lastFace(); ++itFace) {
+
+ }
+}
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index b1d09169f7..f8f5090f0f 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -29,4 +29,10 @@ MTriangle* setHighOrder(MTriangle *t, GFace *gf,
void checkHighOrderTriangles(const char* cc, GModel *m,
std::vector<MElement*> &bad, double &minJGlob);
+struct distanceFromMeshToGeometry_t {
+ std::map<GEntity*, double> d_max, d2;
+};
+
+void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t &dist);
+
#endif
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 1133b968c2..42ce669dad 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -892,7 +892,7 @@ void bowyerWatson(GFace *gf, int MAXPNT)
{
FieldManager *fields = gf->model()->getFields();
BoundaryLayerField *blf = 0;
- if(fields->getBackgroundField() > 0){
+ if(fields->getBoundaryLayerField() > 0){
Field *bl_field = fields->get(fields->getBoundaryLayerField());
blf = dynamic_cast<BoundaryLayerField*> (bl_field);
if (blf && !blf->iRecombine) quadsToTriangles(gf,10000);
@@ -1208,7 +1208,7 @@ void bowyerWatsonFrontal(GFace *gf)
{
FieldManager *fields = gf->model()->getFields();
BoundaryLayerField *blf = 0;
- if(fields->getBackgroundField() > 0){
+ if(fields->getBoundaryLayerField() > 0){
Field *bl_field = fields->get(fields->getBoundaryLayerField());
blf = dynamic_cast<BoundaryLayerField*> (bl_field);
if (blf && !blf->iRecombine)quadsToTriangles(gf,10000);
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 807a5ed675..0197263f11 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -2,6 +2,7 @@
#include <algorithm>
#include "OptHomMesh.h"
#include "OptHOM.h"
+#include "GmshMessage.h"
#include "GmshConfig.h"
#ifdef HAVE_BFGS
@@ -15,18 +16,15 @@
// Constructor
-OptHOM::OptHOM(GEntity *ge,const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method) :
- mesh(ge, els, toFix, method)
+OptHOM::OptHOM(GEntity *ge, std::set<MVertex*> & toFix, int method) :
+ mesh(ge, toFix, method)
{
-}
+};
// Contribution of the element Jacobians to the objective function value and gradients (2D version)
bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
{
- minJac = 1.e300;
- maxJac = -1.e300;
-
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
mesh.scaledJac(iEl,sJ);
@@ -37,8 +35,6 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
const double f1 = compute_f1(sJ[l]);
for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++)
gradObj[mesh.indPCEl(iEl,iPC)] += f1*gSJ[mesh.indGSJ(iEl,l,iPC)];
- minJac = std::min(minJac,sJ[l]);
- maxJac = std::max(maxJac,sJ[l]);
}
}
@@ -52,22 +48,13 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj)
{
- maxDist = 0;
- avgDist = 0;
- int nbBnd = 0;
-
for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
const double Factor = invLengthScaleSq*(mesh.forced(iFV) ? Fact : Fact2);
- const double dSq = mesh.distSq(iFV), dist = sqrt(dSq);
- Obj += Factor * dSq;
+ Obj += Factor * mesh.distSq(iFV);
std::vector<double> gDSq(mesh.nPCFV(iFV));
mesh.gradDistSq(iFV,gDSq);
for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) gradObj[mesh.indPCFV(iFV,iPC)] += Factor*gDSq[iPC];
- maxDist = std::max(maxDist, dist);
- avgDist += dist;
- nbBnd++;
}
- if (nbBnd != 0) avgDist /= nbBnd;
return true;
@@ -83,11 +70,8 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
Obj = 0.;
for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
- if (minJac > barrier) printf("INFO: reached minimum Jacobian requirement, setting null gradient\n");
- else {
- addJacObjGrad(Obj, gradObj);
- addDistObjGrad(lambda, lambda2, Obj, gradObj);
- }
+ addJacObjGrad(Obj, gradObj);
+ addDistObjGrad(lambda, lambda2, Obj, gradObj);
}
@@ -95,30 +79,38 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj, void *HOInst)
{
- (static_cast<OptHOM*>(HOInst))->evalObjGrad(x, Obj, gradObj);
+ OptHOM &HO = *static_cast<OptHOM*> (HOInst);
+ HO.evalObjGrad(x, Obj, gradObj);
+ double distMaxBnd, distAvgBnd, minJac, maxJac;
+ HO.getDistances(distMaxBnd, distAvgBnd, minJac, maxJac);
+ if (minJac > HO.barrier_min && maxJac < HO.barrier_max) {
+ for (int i = 0; i < gradObj.length(); ++i) {
+ gradObj[i] = 0;
+ }
+ }
}
-void OptHOM::recalcJacDist()
+void OptHOM::getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac)
{
- maxDist = 0;
- avgDist = 0;
+ distMaxBND = 0;
+ distAvgBND = 0;
int nbBnd = 0;
for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
if (mesh.forced(iFV)) {
double dSq = mesh.distSq(iFV);
- maxDist = std::max(maxDist, sqrt(dSq));
- avgDist += sqrt(dSq);
+ distMaxBND = std::max(distMaxBND, sqrt(dSq));
+ distAvgBND += sqrt(dSq);
nbBnd++;
}
}
- if (nbBnd != 0) avgDist /= nbBnd;
+ if (nbBnd != 0) distAvgBND /= nbBnd;
- minJac = 1.e300;
- maxJac = -1.e300;
+ minJac = 1000;
+ maxJac = -1000;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
mesh.scaledJac(iEl,sJ);
@@ -136,105 +128,83 @@ void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj)
{
iter++;
- if (iter % progressInterv == 0) printf("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E maxJ = %12.5E Max D = %12.5E Avg D = %12.5E\n", iter, Obj, Obj/initObj, minJac, maxJac, maxDist, avgDist);
-
-}
+ if (iter % progressInterv == 0) {
+ double maxD, avgD, minJ, maxJ;
+ getDistances(maxD, avgD, minJ, maxJ);
+ printf("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E maxJ = %12.5E Max D = %12.5E Avg D = %12.5E\n", iter, Obj, Obj/initObj, minJ, maxJ, maxD, avgD);
+ Msg::Debug("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E maxJ = %12.5E Max D = %12.5E Avg D = %12.5E", iter, Obj, Obj/initObj, minJ, maxJ, maxD, avgD);
+ }
-void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst)
-{
- ((OptHOM*)HOInst)->printProgress(x,Obj);
}
-void OptHOM::calcScale(alglib::real_1d_array &scale)
+void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst)
{
-
- scale.setlength(mesh.nPC());
-
- // Calculate scale
- for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
- std::vector<double> scaleFV(mesh.nPCFV(iFV),1.);
- mesh.pcScale(iFV,scaleFV);
- for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) scale[mesh.indPCFV(iFV,iPC)] = scaleFV[iPC];
- }
-
- // Normalize scale vector (otherwise ALGLIB routines may fail)
- double scaleNormSq = 0.;
- for (int i = 0; i < mesh.nPC(); i++) scaleNormSq += scale[i]*scale[i];
- const double scaleNorm = sqrt(scaleNormSq);
- for (int i = 0; i < mesh.nPC(); i++) scale[i] /= scaleNorm;
-
+ ((OptHOM*)HOInst)->printProgress(x,Obj);
}
-
void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax)
{
static const double EPSG = 0.;
- static const double EPSF = 0.;
+ static const double EPSF = 1.e-8;
+// static const double EPSF = 0.;
static const double EPSX = 0.;
- static int OPTMETHOD = 1;
+// const double EPSX = x.length()*1.e-4/sqrt(invLengthScaleSq);
+// std::cout << "DEBUG: EPSX = " << EPSX << ", EPSX/x.length() = " << EPSX/x.length() << std::endl;
- std::cout << "--- Optimization pass with jacBar = " << jacBar << ", lambda = " << lambda << ", lambda2 = " << lambda2 << std::endl;
+// double iGONorm = 0;
+// for (int i=0; i<initGradObj.length(); i++) iGONorm += initGradObj[i]*initGradObj[i];
+// const double EPSG = 1.e-2*sqrt(iGONorm);
- iter = 0;
+ Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar);
- int iterationscount = 0, nfev = 0, terminationtype = -1;
- if (OPTMETHOD == 1) {
- alglib::mincgstate state;
- alglib::mincgreport rep;
- mincgcreate(x, state);
- alglib::real_1d_array scale;
- calcScale(scale);
- mincgsetscale(state,scale);
- mincgsetprecscale(state);
- mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
- mincgsetxrep(state, true);
- alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
- mincgresults(state, x, rep);
- iterationscount = rep.iterationscount;
- nfev = rep.nfev;
- terminationtype = rep.terminationtype;
- }
- else {
- alglib::minlbfgsstate state;
- alglib::minlbfgsreport rep;
- minlbfgscreate(3, x, state);
- alglib::real_1d_array scale;
- calcScale(scale);
- minlbfgssetscale(state,scale);
- minlbfgssetprecscale(state);
- minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
- minlbfgssetxrep(state, true);
- alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
- minlbfgsresults(state, x, rep);
- iterationscount = rep.iterationscount;
- nfev = rep.nfev;
- terminationtype = rep.terminationtype;
- }
+// alglib::minlbfgsstate state;
+// alglib::minlbfgsreport rep;
+ alglib::mincgstate state;
+ alglib::mincgreport rep;
- std::cout << "Optimization finalized after " << iterationscount << " iterations (" << nfev << " functions evaluations)";
- switch(int(terminationtype)) {
- case 1: std::cout << ", because relative function improvement is no more than EpsF"; break;
- case 2: std::cout << ", because relative step is no more than EpsX"; break;
- case 4: std::cout << ", because gradient norm is no more than EpsG"; break;
- case 5: std::cout << ", because the maximum number of steps was taken"; break;
- default: std::cout << " with code " << int(terminationtype); break;
- }
- std::cout << "." << std::endl;
+// minlbfgscreate(3, x, state);
+// minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
+// minlbfgssetxrep(state, true);
+ mincgcreate(x, state);
+ mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
+ mincgsetxrep(state, true);
+
+ iter = 0;
+
+// alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
+ alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
+
+// minlbfgsresults(state, x, rep);
+ mincgresults(state, x, rep);
+
+ // std::cout << "Optimization finalized after " << rep.iterationscount << " iterations (" << rep.nfev << " functions evaluations)";
+ // switch(int(rep.terminationtype)) {
+// case -2: std::cout << ", because rounding errors prevented further improvement"; break;
+// case -1: std::cout << ", because incorrect parameters were specified"; break;
+// case 1: std::cout << ", because relative function improvement is no more than EpsF"; break;
+// case 2: std::cout << ", because relative step is no more than EpsX"; break;
+// case 4: std::cout << ", because gradient norm is no more than EpsG"; break;
+// case 5: std::cout << ", because the maximum number of steps was taken"; break;
+// case 7: std::cout << ", because stopping conditions are too stringent, further improvement is impossible"; break;
+// default: std::cout << " with code " << int(rep.terminationtype); break;
+// }
+// std::cout << "." << std::endl;
}
-int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int pInt, int itMax)
+int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double b_max, int pInt, int itMax, double &minJ, double &maxJ)
{
- barrier = barrier_;
+ barrier_min = b_min;
+ barrier_max = b_max;
progressInterv = pInt;
// powM = 4;
// powP = 3;
@@ -253,13 +223,14 @@ int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int
mesh.getUvw(x.getcontent());
// Calculate initial performance
- recalcJacDist();
- initMaxDist = maxDist;
- initAvgDist = avgDist;
+ // double minJ, maxJ;
+ double initMaxD, initAvgD;
+ getDistances(initMaxD, initAvgD, minJ, maxJ);
- const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
+ const double jacBarStart = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
jacBar = jacBarStart;
setBarrierTerm(jacBarStart);
+ // std::cout << "DEBUG: jacBarStart = " << jacBarStart << std::endl;
// Calculate initial objective function value and gradient
initObj = 0.;
@@ -268,27 +239,30 @@ int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int
for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
evalObjGrad(x, initObj, gradObj);
- std::cout << "Initial mesh: Obj = " << initObj << ", minJ = " << minJac << ", maxJ = " << maxJac << ", maxD = " << initMaxDist << ", avgD = " << initAvgDist << std::endl;
- std::cout << "Start optimizing " << mesh.nEl() << " elements (" << mesh.nVert() << " vertices, "
- << mesh.nFV() << " free vertices, " << mesh.nPC() << " variables) with barrier = " << barrier << std::endl;
+ // std::cout << "Start optimizing with barrier = " << barrier << std::endl;
- while (minJac < barrier) {
+ int ITER = 0;
+ while (minJ < barrier_min) {
OptimPass(x, gradObj, itMax);
- jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
+ double dumMaxD, dumAvgD;
+ getDistances(dumMaxD, dumAvgD, minJ, maxJ);
+ jacBar = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
setBarrierTerm(jacBar);
+ if (ITER ++ > 5) break;
}
-// for (int i = 0; i<3; i++) {
-// lambda *= 10;
-// OptimPass(x, gradObj, itMax);
-// }
+ // for (int i = 0; i<3; i++) {
+ // lambda *= 100;
+ // OptimPass(x, gradObj, itMax);
+ // }
- double finalObj = 0.;
+ double finalObj = 0., finalMaxD, finalAvgD;
evalObjGrad(x, finalObj, gradObj);
- std::cout << "Final mesh: Obj = " << finalObj << ", maxD = " << maxDist << ", avgD = " << avgDist << ", minJ = " << minJac << ", maxJ = " << maxJac << std::endl;
+ getDistances(finalMaxD, finalAvgD, minJ, maxJ);
+ Msg::Info("Optimization finished : Avg distance to bnd = %12.5E MinJac %12.5E MaxJac %12.5E",finalAvgD,minJ,maxJ);
- return 0;
+ return (minJ > barrier_min && maxJ < barrier_max );
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index 39e49aabf3..b05b5d8fa6 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -21,15 +21,18 @@ public:
Mesh mesh;
- OptHOM(GEntity *gf, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method);
- void recalcJacDist();
- inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
- int optimize(double lambda, double lambda2, double barrier, int pInt, int itMax); // optimize one list of elements
+ OptHOM(GEntity *gf, std::set<MVertex*> & toFix, int method);
+ void getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac);
+ int optimize(double lambda, double lambda2, double barrier_min, double barrier_max,int pInt, int itMax, double &minJ, double &maxJ); // optimize one list of elements
+ // OptHOM(GEntity *gf, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method);
+ // void recalcJacDist();
+ // inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
+ // int optimize(double lambda, double lambda2, double barrier, int pInt, int itMax); // optimize one list of elements
void updateMesh(const alglib::real_1d_array &x);
void evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj);
void printProgress(const alglib::real_1d_array &x, double Obj);
- double barrier;
+ double barrier_min, barrier_max;
private:
@@ -81,10 +84,10 @@ inline double OptHOM::compute_f1(double v)
-void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
-{
- minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist;
-}
+//void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
+//{
+// minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist;
+//}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 0a598eeb11..7ea06cfa7e 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -1,17 +1,15 @@
+#include "GmshMessage.h"
#include "GRegion.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "ParamCoord.h"
#include "OptHomMesh.h"
-
-
+#include "JacobianBasis.h"
std::map<int, std::vector<double> > Mesh::_jacBez;
-
-
-std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
+static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
{
int nbNodes = lagrange->points.size1();
@@ -39,7 +37,7 @@ std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezie
for (int j = 0; j < nbNodes; j++) {
double Jij = dPsi(i, 0) * dPsi(j, 1) - dPsi(i, 1) * dPsi(j,0);
for (int l = 0; l < bezier->points.size1(); l++) {
- JB[indJB2DBase(nbNodes,l,i,j)] += bezier->matrixLag2Bez(l, k) * Jij;
+ JB[(l * nbNodes + i) * nbNodes + j] += bezier->matrixLag2Bez(l, k) * Jij;
}
}
}
@@ -53,7 +51,7 @@ std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezie
+ (dPsi(j, 2) * dPsi(m, 0) - dPsi(j, 0) * dPsi(m, 2)) * dPsi(i, 1)
+ (dPsi(j, 0) * dPsi(m, 1) - dPsi(j, 1) * dPsi(m, 0)) * dPsi(i, 2);
for (int l = 0; l < bezier->points.size1(); l++) {
- JB[indJB3DBase(nbNodes,l,i,j,m)] += bezier->matrixLag2Bez(l, k) * Jijm;
+ JB[(l * nbNodes + (i * nbNodes + j)) * nbNodes + m] += bezier->matrixLag2Bez(l, k) * Jijm;
}
}
}
@@ -63,9 +61,7 @@ std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezie
return JB;
}
-
-
-Mesh::Mesh(GEntity *ge,const std::set<MElement*> &elements, std::set<MVertex*> & toFix, int method):
+Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
_ge(ge)
{
_dim = _ge->dim();
@@ -73,41 +69,39 @@ Mesh::Mesh(GEntity *ge,const std::set<MElement*> &elements, std::set<MVertex*> &
if (method & METHOD_PHYSCOORD) {
if (_dim == 2) _pc = new ParamCoordPhys2D;
else _pc = new ParamCoordPhys3D;
- std::cout << "METHOD: Using physical coordinates" << std::endl;
+ Msg::Debug("METHOD: Using physical coordinates");
}
else if (method & METHOD_SURFCOORD) {
if (_dim == 2) {
_pc = new ParamCoordSurf(_ge);
- std::cout << "METHOD: Using surface parametric coordinates" << std::endl;
+ Msg::Debug("METHOD: Using surface parametric coordinates");
}
- else std::cout << "ERROR: Surface parametric coordinates only for 2D optimization" << std::endl;
+ Msg::Error("ERROR: Surface parametric coordinates only for 2D optimization");
}
else {
_pc = new ParamCoordParent;
- std::cout << "METHOD: Using parent parametric coordinates" << std::endl;
+ Msg::Debug("METHOD: Using parent parametric coordinates");
}
- if (method & METHOD_RELAXBND) std::cout << "METHOD: Relaxing boundary vertices" << std::endl;
- else if (method & METHOD_FIXBND) std::cout << "METHOD: Fixing all boundary vertices" << std::endl;
- else std::cout << "METHOD: Fixing vertices on geometric points and \"toFix\" boundary" << std::endl;
+ 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
_nPC = 0;
- int nElements = elements.size();
- _el.resize(nElements);
- _el2FV.resize(nElements);
- _el2V.resize(nElements);
- _nBezEl.resize(nElements);
- _nNodEl.resize(nElements);
- _indPCEl.resize(nElements);
- int iEl = 0;
- for (std::set<MElement *>::iterator it = elements.begin(); it != elements.end(); ++it, ++iEl) {
- MElement *el = *it;
+ _el.resize(_ge->getNumMeshElements());
+ _el2FV.resize(_ge->getNumMeshElements());
+ _el2V.resize(_ge->getNumMeshElements());
+ _nBezEl.resize(_ge->getNumMeshElements());
+ _nNodEl.resize(_ge->getNumMeshElements());
+ _indPCEl.resize(_ge->getNumMeshElements());
+ for (int iEl = 0; iEl < nEl(); iEl++) {
+ MElement *el = _ge->getMeshElement(iEl);
_el[iEl] = el;
const polynomialBasis *lagrange = el->getFunctionSpace();
const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
if (_jacBez.find(lagrange->type) == _jacBez.end()) {
- _jacBez[lagrange->type] = computeJB(lagrange, bezier);
+ _jacBez[lagrange->type] = _computeJB(lagrange, bezier);
}
_nBezEl[iEl] = bezier->points.size1();
_nNodEl[iEl] = lagrange->points.size1();
@@ -151,16 +145,14 @@ Mesh::Mesh(GEntity *ge,const std::set<MElement*> &elements, std::set<MVertex*> &
}
if ((_dim == 2) && (method & METHOD_PROJJAC)) {
projJac = true;
- std::cout << "METHOD: Using projected Jacobians" << std::endl;
+ Msg::Debug("METHOD: Using projected Jacobians");
}
else {
projJac = false;
- std::cout << "METHOD: Using usual Jacobians" << std::endl;
+ Msg::Debug("METHOD: Using usual Jacobians");
}
}
-
-
SVector3 Mesh::getNormalEl(int iEl)
{
@@ -264,7 +256,6 @@ void Mesh::updateMesh(const double *it)
}
-
void Mesh::distSqToStraight(std::vector<double> &dSq)
{
std::vector<SPoint3> sxyz(nVert());
@@ -333,7 +324,7 @@ void Mesh::scaledJac(int iEl, std::vector<double> &sJ)
}
sJ[l] *= n.z();
}
- }
+ }
else {
for (int l = 0; l < _nBezEl[iEl]; l++) {
sJ[l] = 0.;
@@ -430,8 +421,8 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
_pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < _nBezEl[iEl]; l++) {
gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
- if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
- if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
+ if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][1];
+ if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][2];
}
iPC += _nPCFV[iFVi];
}
@@ -440,27 +431,6 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
}
-
-
-void Mesh::pcScale(int iFV, std::vector<double> &scale)
-{
-
- // Calc. derivative of x, y & z w.r.t. parametric coordinates
- const SPoint3 dX(1.,0.,0.), dY(0.,1.,0.), dZ(0.,0.,1.);
- SPoint3 gX, gY, gZ;
- _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dX,gX);
- _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dY,gY);
- _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dZ,gZ);
-
- // Scale = inverse norm. of vector (dx/du, dy/du, dz/du)
- scale[0] = 1./sqrt(gX[0]*gX[0]+gY[0]*gY[0]+gZ[0]*gZ[0]);
- if (_nPCFV[iFV] >= 2) scale[1] = 1./sqrt(gX[1]*gX[1]+gY[1]*gY[1]+gZ[1]*gZ[1]);
- if (_nPCFV[iFV] == 3) scale[2] = 1./sqrt(gX[2]*gX[2]+gY[2]*gY[2]+gZ[2]*gZ[2]);
-
-}
-
-
-
void Mesh::writeMSH(const char *filename)
{
FILE *f = fopen(filename, "w");
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 5636311267..2e2741fbb4 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(GEntity *ge, std::set<MVertex*> & toFix, int method);
inline const int &nPC() { return _nPC; }
inline int nVert() { return _vert.size(); }
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 0c171c6308..88acf1dd74 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -8,41 +8,351 @@
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MLine.h"
#include "highOrderTools.h"
#include "OptHomMesh.h"
#include "OptHOM.h"
#include <stack>
-static std::set<MVertex*> filter(GEntity *gf, int nbLayers, double _qual, std::set<MElement*> &elements, bool onlytheworst = false)
+// get all elements that are neighbors
+
+double distMaxStraight (MElement *el){
+ const polynomialBasis *lagrange = el->getFunctionSpace();
+ const polynomialBasis *lagrange1 = 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];
+ lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1), lagrange->points(i, 2), 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");
+
+ int numVertices = gm->indexMeshVertices(true);
+ std::vector<GEntity*> entities;
+ gm->getEntities(entities);
+ fprintf(f,"%d %d\n",numVertices,dim);
+ for(unsigned int i = 0; i < entities.size(); i++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
+ MVertex *v = entities[i]->mesh_vertices[j];
+ if (dim == 2)fprintf(f,"%d %22.15E %22.15E\n",v->getIndex(),v->x(),v->y());
+ else if (dim == 3)fprintf(f,"%d %22.15E %22.15E %22.5E\n",v->getIndex(),v->x(),v->y(),v->z());
+ }
+
+ if (dim == 2){
+ int nt = 0;
+ int order = 0;
+ for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf){
+ std::vector<MTriangle*> &tris = (*itf)->triangles;
+ nt += tris.size();
+ if (tris.size())order = tris[0]->getPolynomialOrder();
+ }
+ fprintf(f,"%d %d\n",nt,(order+1)*(order+2)/2);
+ int count = 1;
+ for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf){
+ std::vector<MTriangle*> &tris = (*itf)->triangles;
+ for (int i=0;i<tris.size();i++){
+ MTriangle *t = tris[i];
+ fprintf(f,"%d ",count++);
+ for (int j=0;j<t->getNumVertices();j++){
+ fprintf(f,"%d ",t->getVertex(j)->getIndex());
+ }
+ fprintf(f,"\n");
+ }
+ }
+ int ne = 0;
+ for (GModel::eiter ite = gm->firstEdge(); ite != gm->lastEdge(); ++ite){
+ std::vector<MLine*> &l = (*ite)->lines;
+ ne += l.size();
+ }
+ fprintf(f,"%d %d\n",ne,(order+1));
+ count = 1;
+ for (GModel::eiter ite = gm->firstEdge(); ite != gm->lastEdge(); ++ite){
+ std::vector<MLine*> &l = (*ite)->lines;
+ for (int i=0;i<l.size();i++){
+ MLine *t = l[i];
+ fprintf(f,"%d ",count++);
+ for (int j=0;j<t->getNumVertices();j++){
+ fprintf(f,"%d ",t->getVertex(j)->getIndex());
+ }
+ fprintf(f,"%d \n",(*ite)->tag());
+ }
+ }
+ }
+ fclose(f);
+}
+
+
+
+static void getTopologicalNeighbors(int nbLayers,
+ const std::vector<MElement*> &all,
+ const std::vector<MElement*> &elements,
+ std::set<MElement*> &result){
+
+
+ std::set<MVertex*> touched;
+
+ for (int i = 0; i < elements.size(); i++)
+ for (int j=0;j<elements[i]->getNumVertices();j++)
+ touched.insert(elements[i]->getVertex(j));
+
+ for (int layer = 0; layer < nbLayers; layer++) {
+ for (int i = 0; i < all.size(); i++) {
+ MElement *t = all[i];
+ bool add_ = false;
+ for (int j = 0; j < t->getNumVertices(); j++) {
+ if (touched.find(t->getVertex(j)) != touched.end()) {
+ add_ = true;
+ }
+ }
+ if (add_) result.insert(t);
+ }
+ for (std::set<MElement*>::iterator it = result.begin(); it != result.end(); ++it) {
+ for (int j = 0; j < (*it)->getNumVertices(); j++) {
+ touched.insert((*it)->getVertex(j));
+ }
+ }
+ }
+ // printf("%d %d %d %d\n",all.size(),elements.size(), result.size(),nbLayers);
+ // exit(1);
+}
+
+static void getGeometricalNeighbors (MElement *e, const std::vector<MElement*> &all, double F, std::set<MElement*> &result) {
+
+ double R = distMaxStraight (e);
+
+ SPoint3 p = e->barycenter_infty();
+ for (int i = 0; i < all.size(); i++) {
+ MElement *e = all[i];
+ double dist = p.distance(e->barycenter_infty());
+ if (dist < R*F)result.insert(e);
+ }
+}
+
+static void intersection (const std::set<MElement*> &a, const std::set<MElement*> &b, std::set<MElement*> &result){
+ for (std::set<MElement*>::const_iterator it = a.begin() ; it != a.end() ; ++it ){
+ if (b.find(*it) != b.end()){
+ result.insert(*it);
+ }
+ }
+}
+
+static void computeVertices (const std::set<MElement*> &a, std::set<MVertex*> &result){
+ for (std::set<MElement*>::const_iterator it = a.begin() ; it != a.end() ; ++it ){
+ for (int i=0;i<(*it)->getNumVertices();++i){
+ result.insert((*it)->getVertex(i));
+ }
+ }
+}
+
+static MElement * compare_worst (MElement *a, MElement *b){
+ if (!a)return b;
+ if (!b)return a;
+ if (a->distoShapeMeasure() < b->distoShapeMeasure()) return a;
+ return b;
+}
+
+template <class T>
+MElement* getTheWorstElement (const std::vector<T*> &a) {
+ T *worst = 0;
+ double q = 1.e22;
+ for (int i=0;i<a.size();i++){
+ T *t = a[i];
+ double Q = t->distoShapeMeasure();
+ if (Q < q) {
+ worst = t;q = Q;
+ }
+ }
+ // printf("worst = %12.5E\n",q);
+ return worst;
+}
+
+std::set<MVertex*> filter3D(GRegion *gr, int nbLayers, double _qual)
{
- elements.clear();
std::set<MVertex*> touched;
std::set<MVertex*> not_touched;
- double minQ = 1e22;
- MElement *worst;
- for (unsigned int i = 0; i < gf->getNumMeshElements(); ++i) {
- MElement *el = gf->getMeshElement(i);
- double Q = el->distoShapeMeasure();
- if (Q < _qual || Q > 1.01) {
- if (Q < minQ) {
- worst = el;
- minQ = Q;
+ std::set<MTetrahedron *> ts;
+ for (int i=0;i<gr->tetrahedra.size();i++){
+ MTetrahedron *t = gr->tetrahedra[i];
+ if (t->distoShapeMeasure() < _qual){
+ ts.insert(t);
+ for(int j=0;j<t->getNumVertices();j++)touched.insert(t->getVertex(j));
+ }
+ if (ts.size() == 1)break;
+ }
+
+ printf("FILTER3D : %d bad tets found among %6d\n",ts.size(),gr->tetrahedra.size());
+
+ // add layers of elements around bad elements
+ for (int layer = 0;layer<nbLayers; layer++){
+ for (int i=0;i<gr->tetrahedra.size();i++){
+ MTetrahedron *t = gr->tetrahedra[i];
+ bool add_ = false;
+ for(int j=0;j<t->getNumVertices();j++){
+ if(touched.find(t->getVertex(j)) != touched.end()){
+ add_ = true;
+ }
}
- if (!onlytheworst) {
- elements.insert(el);
+ if (add_)ts.insert(t);
+ }
+
+ for (std::set<MTetrahedron*>::iterator it = ts.begin(); it != ts.end() ; ++it){
+ for(int j=0;j<(*it)->getNumVertices();j++){
+ touched.insert((*it)->getVertex(j));
}
}
}
- if (onlytheworst) {
- elements.insert(worst);
+
+ for (int i=0;i<gr->tetrahedra.size();i++){
+ if (ts.find(gr->tetrahedra[i]) == ts.end()){
+ for(int j=0;j<gr->tetrahedra[i]->getNumVertices();j++){
+ if (touched.find(gr->tetrahedra[i]->getVertex(j)) != touched.end())
+ not_touched.insert(gr->tetrahedra[i]->getVertex(j));
+ }
+ }
+ }
+
+ gr->tetrahedra.clear();
+ gr->tetrahedra.insert(gr->tetrahedra.begin(),ts.begin(),ts.end());
+
+ printf("FILTER3D : AFTER FILTERING %d tets remain %d nodes on the boundary\n",gr->tetrahedra.size(),not_touched.size());
+ return not_touched;
+}
+
+
+std::set<MVertex*> filter2D_boundaryLayer(GFace *gf,
+ int nbLayers,
+ double _qual_min,
+ double _qual_max,
+ double F ) {
+ MElement *worst = compare_worst (getTheWorstElement(gf->triangles),
+ getTheWorstElement(gf->quadrangles));
+
+ std::vector<MElement*> vworst; vworst.push_back(worst);
+ std::vector<MElement*> all;
+ all.insert(all.begin(),gf->triangles.begin(),gf->triangles.end());
+ all.insert(all.begin(),gf->quadrangles.begin(),gf->quadrangles.end());
+ std::set<MElement*> result1;
+ getTopologicalNeighbors(nbLayers, all, vworst,result1);
+ std::set<MElement*> result2;
+ getGeometricalNeighbors (worst, all, F, result2);
+ std::set<MElement*> result;
+ intersection (result1,result2,result);
+ // printf("intsersection(%d,%d) = %d\n",result1.size(),result2.size(),result.size());
+
+
+ gf->triangles.clear();
+ gf->quadrangles.clear();
+ for (std::set<MElement*>::iterator it = result.begin(); it != result.end(); ++it){
+ MElement *e = *it;
+ if (e->getType() == TYPE_QUA)gf->quadrangles.push_back((MQuadrangle*)e);
+ else if (e->getType() == TYPE_TRI)gf->triangles.push_back((MTriangle*)e);
+ }
+
+ std::set<MVertex*> vs;
+ for (int i = 0; i < all.size(); i++) {
+ if (result.find(all[i]) == result.end()) {
+ for (int j = 0; j < all[i]->getNumVertices(); j++) {
+ vs.insert(all[i]->getVertex(j));
+ }
+ }
+ }
+ return vs;
+}
+
+
+std::set<MVertex*> filter3D_boundaryLayer(GRegion *gr,
+ int nbLayers,
+ double _qual_min,
+ double _qual_max,
+ double F ) {
+ MElement *worst = compare_worst (getTheWorstElement(gr->tetrahedra),
+ getTheWorstElement(gr->prisms));
+ worst = compare_worst (worst,getTheWorstElement(gr->hexahedra));
+
+ std::vector<MElement*> vworst; vworst.push_back(worst);
+ std::vector<MElement*> all;
+ all.insert(all.begin(),gr->tetrahedra.begin(),gr->tetrahedra.end());
+ all.insert(all.begin(),gr->prisms.begin(),gr->prisms.end());
+ all.insert(all.begin(),gr->hexahedra.begin(),gr->hexahedra.end());
+ std::set<MElement*> result1;
+ getTopologicalNeighbors(nbLayers, all, vworst,result1);
+ std::set<MElement*> result2;
+ getGeometricalNeighbors (worst, all, F, result2);
+ std::set<MElement*> result;
+ intersection (result1,result2,result);
+ // printf("intsersection(%d,%d) = %d\n",result1.size(),result2.size(),result.size());
+
+
+ gr->tetrahedra.clear();
+ gr->hexahedra.clear();
+ gr->prisms.clear();
+ for (std::set<MElement*>::iterator it = result.begin(); it != result.end(); ++it){
+ MElement *e = *it;
+ if (e->getType() == TYPE_TET)gr->tetrahedra.push_back((MTetrahedron*)e);
+ else if (e->getType() == TYPE_HEX)gr->hexahedra.push_back((MHexahedron*)e);
+ else if (e->getType() == TYPE_PRI)gr->prisms.push_back((MPrism*)e);
}
- for (std::set<MElement*>::iterator it = elements.begin(); it != elements.end(); ++it) {
- MElement &el = **it;
- for (int j = 0; j < el.getNumVertices(); j++)
- touched.insert(el.getVertex(j));
+
+ std::set<MVertex*> vs;
+ for (int i = 0; i < all.size(); i++) {
+ if (result.find(all[i]) == result.end()) {
+ for (int j = 0; j < all[i]->getNumVertices(); j++) {
+ vs.insert(all[i]->getVertex(j));
+ }
+ }
}
- Msg::Info("FILTER2D : %d bad elements found among %6d\n", elements.size(), gf->getNumMeshElements());
+ return vs;
+}
+
+std::set<MVertex*> filter2D(GFace *gf,
+ int nbLayers,
+ double _qual_min,
+ double _qual_max,
+ double F = 1.e22)
+{
+ std::set<MVertex*> touched;
+ std::set<MVertex*> not_touched;
+ std::set<MElement*> elements;
+
+ for (int i = 0; i < gf->triangles.size(); i++) {
+ MTriangle *t = gf->triangles[i];
+ double Q = t->distoShapeMeasure();
+ if (Q < _qual_min || Q > _qual_max) {
+ elements.insert(t);
+ for (int j = 0; j < t->getNumVertices(); j++)
+ touched.insert(t->getVertex(j));
+ }
+ }
+ for (int i = 0; i < gf->quadrangles.size(); i++) {
+ MQuadrangle *q = gf->quadrangles[i];
+ double Q = q->distoShapeMeasure();
+ if (Q < _qual_min || Q > _qual_max) {
+ elements.insert(q);
+ for (int j = 0; j < q->getNumVertices(); j++)
+ touched.insert(q->getVertex(j));
+ }
+ }
+
// add layers of elements around bad elements
for (int layer = 0; layer < nbLayers; layer++) {
for (unsigned int i = 0; i < gf->getNumMeshElements(); ++i) {
@@ -61,16 +371,15 @@ static std::set<MVertex*> filter(GEntity *gf, int nbLayers, double _qual, std::s
}
}
}
- for (unsigned int i = 0; i < gf->getNumMeshElements(); ++i) {
- MElement &el = *gf->getMeshElement(i);
- for (int j = 0; j < el.getNumVertices(); ++j) {
- if (touched.find(el.getVertex(j)) == touched.end()) {
- not_touched.insert(el.getVertex(j));
+
+ for (int i = 0; i < gf->getNumMeshElements(); i++) {
+ if (elements.find(gf->getMeshElement(i)) == elements.end()) {
+ for (int j = 0; j < gf->getMeshElement(i)->getNumVertices(); j++) {
+ if (touched.find(gf->getMeshElement(i)->getVertex(j)) != touched.end()) not_touched.insert(gf->getMeshElement(i)->getVertex(j));
}
}
}
- // printf("FILTER2D : AFTER FILTERING %d triangles %d quads remain %d nodes on the boundary\n", gf->triangles.size(), gf->quadrangles.size(), not_touched.size());
return not_touched;
}
@@ -120,49 +429,83 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
// int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_FIXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_FIXBND | Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
-// int method = Mesh::METHOD_FIXBND;
+ int method;
+ 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;
+
+ printf("p.method = %d\n",p.method);
+
// int method = Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_SURFCOORD | Mesh::METHOD_PROJJAC;
- int method = Mesh::METHOD_PROJJAC;
+// int method = Mesh::METHOD_PROJJAC;
SVector3 BND(gm->bounds().max(), gm->bounds().min());
double SIZE = BND.norm();
+ Msg::Info("High order mesh optimizer starts");
+
+
if (p.dim == 2) {
for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
-
if (p.onlyVisible && !(*itf)->getVisibility())continue;
- std::set<MElement*> elements;
- std::set<MVertex*> toFix = filter(*itf, p.nbLayers, p.BARRIER, elements);
- std::vector<std::set<MElement*> > splitted = splitConnex(elements);
- for (size_t i = 0; i < splitted.size(); ++i) {
- OptHOM *temp = new OptHOM(*itf, splitted[i], toFix, method);
- double distMaxBND, distAvgBND, minJac, maxJac;
- temp->recalcJacDist();
- temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
- if (distMaxBND < 1.e-8 * SIZE && minJac > p.BARRIER) continue;
- p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax);
- temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
- temp->mesh.updateGEntityPositions();
- delete temp;
+ int ITER = 0;
+ Msg::Info("Model face %4d is considered",(*itf)->tag());
+ p.SUCCESS = true;
+ while (1){
+ std::vector<MTriangle*> tris = (*itf)->triangles;
+ std::vector<MQuadrangle*> quas = (*itf)->quadrangles;
+ std::set<MVertex*> toFix;
+
+ if (p.filter == 1) toFix = filter2D_boundaryLayer(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor);
+ else toFix = filter2D(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor);
+
+ OptHOM *temp = new OptHOM(*itf, toFix, method);
+ double distMaxBND, distAvgBND, minJac, maxJac;
+ temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac);
+ if (minJac < 1.e2)Msg::Info("Optimizing a blob of %4d elements minJ %12.5E -- maxJ %12.5E",(*itf)->getNumMeshElements(), minJac, maxJac);
+ (*itf)->triangles = tris;
+ (*itf)->quadrangles = quas;
+ std::ostringstream ossI;
+ ossI << "initial_" << (*itf)->tag() << "ITER_" << ITER << ".msh";
+ temp->mesh.writeMSH(ossI.str().c_str());
+ if (minJac > p.BARRIER_MIN && maxJac < p.BARRIER_MAX) break;
+
+ p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax, p.minJac, p.maxJac);
+ temp->mesh.updateGEntityPositions();
+ if (!p.SUCCESS) break;
+ ITER ++;
+
+ // std::ostringstream ossF;
+ // ossF << "final_" << (*itf)->tag() << ".msh";
+ // temp->mesh.writeMSH(ossF.str().c_str());
}
+ Msg::Info("Model face %4d is done (%d subdomains were computed) SUCCESS %1d",(*itf)->tag(),ITER,p.SUCCESS);
+ Msg::Info("----------------------------------------------------------------");
}
+ exportMeshToDassault (gm,gm->getName() + "_opti.das", 2);
}
else if (p.dim == 3) {
for (GModel::riter itr = gm->firstRegion(); itr != gm->lastRegion(); ++itr) {
if (p.onlyVisible && !(*itr)->getVisibility())continue;
std::vector<MTetrahedron*> tets = (*itr)->tetrahedra;
- std::set<MElement*> elements;
- std::set<MVertex*> toFix = filter(*itr, p.nbLayers, p.BARRIER, elements);
- OptHOM *temp = new OptHOM(*itr, elements, toFix, method);
+ std::set<MVertex*> toFix;
+
+ if (p.filter == 1) toFix = filter3D_boundaryLayer(*itr, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor);
+ else toFix = filter3D(*itr, p.nbLayers, p.BARRIER_MIN);
+
+ OptHOM *temp = new OptHOM(*itr, toFix, method);
double distMaxBND, distAvgBND, minJac, maxJac;
- temp->recalcJacDist();
- temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+ temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac);
+ // temp->recalcJacDist();
+ // temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
// temp->mesh.writeMSH("initial.msh");
// printf("--> Mesh Loaded for GRegion %d : minJ %12.5E -- maxJ %12.5E\n", (*itr)->tag(), minJac, maxJac);
- if (distMaxBND < 1.e-8 * SIZE && minJac > p.BARRIER) continue;
- p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax);
- temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
+ if (minJac > p.BARRIER_MIN && maxJac < p.BARRIER_MAX) continue;
+ p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax, p.minJac, p.maxJac);
temp->mesh.updateGEntityPositions();
(*itr)->tetrahedra = tets;
// temp->mesh.writeMSH("final.msh");
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index c77229145d..76970347a1 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -5,14 +5,20 @@ class GModel;
struct OptHomParameters {
// INPUT ------>
- double BARRIER ; // minimum scaled jcaobian
+ double BARRIER_MIN ; // minimum scaled jcaobian
+ double BARRIER_MAX ; // maximum scaled jcaobian
double weightFixed ; // weight of the energy for fixed nodes
double weightFree ; // weight of the energy for free nodes
+ double STOP;
int nbLayers ; // number of layers taken around a bad element
int dim ; // which dimension to optimize
int itMax ; // max number of iterations in the optimization process
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
+ // 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
+ int filter ; // 0--> standard 1--> boundary layers
// OUTPUT ------>
int SUCCESS ; // 0 --> success , 1 --> Not converged
double minJac, maxJac; // after optimization, range of jacobians
@@ -20,8 +26,8 @@ struct OptHomParameters {
OptHomParameters ()
// default values
- : BARRIER (0.1) , weightFixed (10.), weightFree (1.e-3),
- nbLayers (6) , dim(3) , itMax(1000), onlyVisible(true)
+ : STOP (0.01) , 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), filter(1)
{
}
};
--
GitLab