From 67e8a047902ce6c75239dcefa1b69155fefbb8a7 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 13 Jul 2013 18:13:24 +0000
Subject: [PATCH] more ho cleanup
---
Common/CommandLine.cpp | 19 +++---
Common/Context.cpp | 2 +-
Common/Context.h | 2 +-
Common/DefaultOptions.h | 8 +--
Common/Options.cpp | 8 +--
Common/Options.h | 2 +-
Fltk/optionWindow.cpp | 4 +-
Mesh/HighOrder.cpp | 38 ++----------
Mesh/highOrderTools.cpp | 60 +++++--------------
Solver/multiscaleLaplace.cpp | 3 -
.../HighOrderMeshOptimizer/OptHomElastic.cpp | 0
.../HighOrderMeshOptimizer/OptHomElastic.h | 0
doc/VERSIONS.txt | 2 +
13 files changed, 48 insertions(+), 100 deletions(-)
rename Mesh/highOrderSmoother.cpp => contrib/HighOrderMeshOptimizer/OptHomElastic.cpp (100%)
rename Mesh/highOrderSmoother.h => contrib/HighOrderMeshOptimizer/OptHomElastic.h (100%)
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index f251a1a924..506ddd80c7 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -76,11 +76,8 @@ std::vector<std::pair<std::string, std::string> > GetUsage()
"delquad, del3d, front3d, mmg3d)"));
s.push_back(mp("-smooth int", "Set number of mesh smoothing steps"));
s.push_back(mp("-order int", "Set mesh order (1, ..., 5)"));
- s.push_back(mp("-hoOptimize", "Optimize high order meshes"));
- s.push_back(mp("-hoMindisto float", "Min high-order element quality before optim (0.0->1.0)"));
- s.push_back(mp("-hoNLayers int", "Number of high order element layers to optimize"));
- s.push_back(mp("-hoElasticity float","Poisson ration for elasticity analogy (nu in [-1.0,0.5])"));
s.push_back(mp("-optimize[_netgen]", "Optimize quality of tetrahedral elements"));
+ s.push_back(mp("-optimize_ho", "Optimize high order meshes"));
s.push_back(mp("-optimize_lloyd", "Optimize 2D meshes using Lloyd algorithm"));
s.push_back(mp("-microstructure", "Generate polycrystal Voronoi geometry"));
s.push_back(mp("-clscale float", "Set global mesh element size scaling factor"));
@@ -459,25 +456,29 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.optimizeNetgen = 1;
i++;
}
- else if(!strcmp(argv[i] + 1, "hoOptimize")) {
+ else if(!strcmp(argv[i] + 1, "optimize_ho") ||
+ !strcmp(argv[i] + 1, "hoOptimize")) {
i++;
- opt_mesh_smooth_internal_edges(0, GMSH_SET, 1);
+ opt_mesh_ho_optimize(0, GMSH_SET, 1);
}
- else if(!strcmp(argv[i] + 1, "hoMindisto")) {
+ else if(!strcmp(argv[i] + 1, "ho_mindisto") ||
+ !strcmp(argv[i] + 1, "hoMindisto")) {
i++;
if(argv[i])
opt_mesh_ho_mindisto(0, GMSH_SET, atof(argv[i++]));
else
Msg::Fatal("Missing number");
}
- else if(!strcmp(argv[i] + 1, "hoElasticity")) {
+ else if(!strcmp(argv[i] + 1, "ho_poisson") ||
+ !strcmp(argv[i] + 1, "hoElasticity")) {
i++;
if(argv[i])
opt_mesh_ho_poisson(0, GMSH_SET, atof(argv[i++]));
else
Msg::Fatal("Missing number");
}
- else if(!strcmp(argv[i] + 1, "hoNlayers")) {
+ else if(!strcmp(argv[i] + 1, "ho_nlayers") ||
+ !strcmp(argv[i] + 1, "hoNlayers")) {
i++;
if(argv[i])
opt_mesh_ho_nlayers(0, GMSH_SET, atoi(argv[i++]));
diff --git a/Common/Context.cpp b/Common/Context.cpp
index e590c67d9e..9f2f33dc99 100644
--- a/Common/Context.cpp
+++ b/Common/Context.cpp
@@ -83,7 +83,7 @@ CTX::CTX() : gamepad(0)
mesh.prisms = mesh.pyramids = mesh.hexahedra = 0;
mesh.volumesEdges = mesh.volumesFaces = mesh.surfacesEdges = mesh.surfacesFaces = 0;
mesh.volumesFaces = mesh.surfacesEdges = mesh.surfacesFaces = 0;
- mesh.smoothInternalEdges = mesh.smoothNormals = mesh.reverseAllNormals = 0;
+ mesh.hoOptimize = mesh.smoothNormals = mesh.reverseAllNormals = 0;
mesh.explode = mesh.angleSmoothNormals = 0.;
mesh.numSubEdges = 0;
mesh.colorCarousel = 0;
diff --git a/Common/Context.h b/Common/Context.h
index b90539d770..c41cdc3907 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -41,7 +41,7 @@ struct contextMeshOptions {
int remeshParam, remeshAlgo;
int order, secondOrderLinear, secondOrderIncomplete;
int secondOrderExperimental, meshOnlyVisible;
- int smoothInternalEdges, minCircPoints, minCurvPoints;
+ int hoOptimize, minCircPoints, minCurvPoints;
int smoothNLayers;
double smoothDistoThreshold, smoothPoissonRatio;
int saveAll, saveTri, saveGroupsOfNodes, binary, bdfFieldFormat, saveParametric;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index ef6f041c0d..744de7d00c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -960,12 +960,12 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "Hexahedra" , opt_mesh_hexahedra , 1. ,
"Display mesh hexahedra?" },
{ F|0, "HighOrderNumLayers", opt_mesh_ho_nlayers, 3.,
- "Number of high order mesh elements to consider for optimization."},
- { F|O, "HighOrderOptimize" , opt_mesh_smooth_internal_edges , 0.,
- "Number of smoothing steps of internal edges for high order meshes" },
+ "Number of high order mesh elements to consider for optimization"},
+ { F|O, "HighOrderOptimize" , opt_mesh_ho_optimize , 0.,
+ "Optimize high order meshes?" },
{ F|0, "HighOrderPoissonRatio", opt_mesh_ho_poisson, 0.33,
"Poisson ratio of the material used in the elastic smoother for high order meshes."
- "Must be between -1.0 and 0.5, excluded."},
+ "Must be between -1.0 and 0.5, excluded"},
{ F|O, "HighOrderSmoothingThreshold", opt_mesh_ho_mindisto, 0.5,
"Distorition threshold when choosing which high order element to optimize."},
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 1a6c89827e..5e6e48f817 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5590,16 +5590,16 @@ double opt_mesh_order(OPT_ARGS_NUM)
return CTX::instance()->mesh.order;
}
-double opt_mesh_smooth_internal_edges(OPT_ARGS_NUM)
+double opt_mesh_ho_optimize(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
- CTX::instance()->mesh.smoothInternalEdges = (int)val;
+ CTX::instance()->mesh.hoOptimize = (int)val;
#if defined(HAVE_FLTK)
if(FlGui::available() && (action & GMSH_GUI))
FlGui::instance()->options->mesh.butt[3]->value
- (CTX::instance()->mesh.smoothInternalEdges);
+ (CTX::instance()->mesh.hoOptimize);
#endif
- return CTX::instance()->mesh.smoothInternalEdges;
+ return CTX::instance()->mesh.hoOptimize;
}
double opt_mesh_ho_nlayers(OPT_ARGS_NUM)
diff --git a/Common/Options.h b/Common/Options.h
index 46ca772bf6..3b156db98d 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -446,7 +446,7 @@ double opt_mesh_allow_swap_edge_angle(OPT_ARGS_NUM);
double opt_mesh_min_curv_points(OPT_ARGS_NUM);
double opt_mesh_multiple_passes(OPT_ARGS_NUM);
double opt_mesh_order(OPT_ARGS_NUM);
-double opt_mesh_smooth_internal_edges(OPT_ARGS_NUM);
+double opt_mesh_ho_optimize(OPT_ARGS_NUM);
double opt_mesh_ho_nlayers(OPT_ARGS_NUM);
double opt_mesh_ho_mindisto(OPT_ARGS_NUM);
double opt_mesh_ho_poisson(OPT_ARGS_NUM);
diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index 5c12f19805..8d3eeaf74b 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -492,7 +492,9 @@ static void mesh_options_ok_cb(Fl_Widget *w, void *data)
opt_mesh_optimize_netgen(0, GMSH_SET, o->mesh.butt[24]->value());
opt_mesh_remove_4_triangles(0,GMSH_SET, o->mesh.butt[25]->value());
opt_mesh_order(0, GMSH_SET, o->mesh.value[3]->value());
- opt_mesh_smooth_internal_edges(0, GMSH_SET, o->mesh.butt[3]->value());
+
+ opt_mesh_ho_optimize(0, GMSH_SET, o->mesh.butt[3]->value());
+
opt_mesh_second_order_incomplete(0, GMSH_SET, o->mesh.butt[4]->value());
opt_mesh_points(0, GMSH_SET, o->mesh.butt[6]->value());
opt_mesh_lines(0, GMSH_SET, o->mesh.butt[7]->value());
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 71dde47bea..ffcf0b3b23 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1252,7 +1252,7 @@ void getMeshInfoForHighOrder(GModel *gm, int &meshOrder, bool &complete,
}
}
-void ElasticAnalogy ( GModel *m, double threshold, bool onlyVisible)
+void ElasticAnalogy(GModel *m, double threshold, bool onlyVisible)
{
bool CAD, complete;
int meshOrder;
@@ -1266,15 +1266,14 @@ void ElasticAnalogy ( GModel *m, double threshold, bool onlyVisible)
checkHighOrderTriangles("Surface mesh", m, bad, worst);
{
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
- if (onlyVisible && !(*it)->getVisibility())continue;
+ if (onlyVisible && !(*it)->getVisibility()) continue;
std::vector<MElement*> v;
v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
- if (CAD)hot.applySmoothingTo(v, (*it));
+ if (CAD) hot.applySmoothingTo(v, (*it));
else hot.applySmoothingTo(v, 1.e32, false);
}
}
- // hot.ensureMinimumDistorsion(0.1);
checkHighOrderTriangles("Final surface mesh", m, bad, worst);
checkHighOrderTetrahedron("Volume Mesh", m, bad, worst);
@@ -1288,11 +1287,9 @@ void ElasticAnalogy ( GModel *m, double threshold, bool onlyVisible)
hot.applySmoothingTo(v,1.e32,false);
}
}
-
- // m->writeMSH("CORRECTED.msh");
+ checkHighOrderTetrahedron("File volume Mesh", m, bad, worst);
}
-
void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisible)
{
// replace all the elements in the mesh with second order elements
@@ -1346,16 +1343,6 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
}
- // highOrderTools hot(m);
-
- // now we smooth mesh the internal vertices of the faces
- // we do that model face by model face
- std::vector<MElement*> bad;
- double worst;
-
- // printJacobians(m, "smoothness_b.pos");
- // m->writeMSH("RAW.msh");
-
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
Msg::Info("Meshing volume %d order %d", (*it)->tag(), order);
Msg::ProgressMeter(++counter, nTot, false, msg);
@@ -1365,23 +1352,10 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
double t2 = Cpu();
- // printJacobians(m, "smoothness.pos");
-
+ std::vector<MElement*> bad;
+ double worst;
checkHighOrderTriangles("Surface mesh", m, bad, worst);
- if (!linear && CTX::instance()->mesh.smoothInternalEdges){
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
- if (onlyVisible && !(*it)->getVisibility())continue;
- std::vector<MElement*> v;
- v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
- v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
- //hot.applySmoothingTo(v, (*it));
- // hot.applySmoothingTo(v, .1,0);
- }
- // hot.ensureMinimumDistorsion(0.1);
- checkHighOrderTriangles("Final surface mesh", m, bad, worst);
- }
checkHighOrderTetrahedron("Volume Mesh", m, bad, worst);
- // m->writeMSH("CORRECTED.msh");
Msg::StatusBar(true, "Done meshing order %d (%g s)", order, t2 - t1);
}
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index 49200f113e..0ff338032e 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -174,19 +174,6 @@ void highOrderTools::ensureMinimumDistorsion (double threshold)
ensureMinimumDistorsion(v,threshold);
}
-/*
-void highOrderTools::applySmoothingTo(GRegion *gr)
-{
- std::vector<MElement*> v;
- v.insert(v.begin(), gr->tetrahedra.begin(),gr->tetrahedra.end());
- v.insert(v.begin(), gr->hexahedra.begin(),gr->hexahedra.end());
- v.insert(v.begin(), gr->prisms.begin(),gr->prisms.end());
- Msg::Info("Smoothing high order mesh : model region %d (%d elements)", gr->tag(),
- v.size());
- applySmoothingTo(v,gr);
-}
-*/
-
void highOrderTools::computeMetricInfo(GFace *gf,
MElement *e,
fullMatrix<double> &J,
@@ -260,7 +247,8 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
}
if (!v.size()) return;
- Msg::Info("Smoothing high order mesh : model face %d (%d elements considered in the elastic analogy, worst mapping %12.5E, %3d bad elements)", gf->tag(),
+ Msg::Info("Smoothing high order mesh : model face %d (%d elements considered in "
+ "the elastic analogy, worst mapping %12.5E, %3d bad elements)", gf->tag(),
v.size(),minD,numBad);
addOneLayer(all, v, layer);
@@ -323,7 +311,6 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
}
else{
SVector3 p = getTL(*it);
- // printf("%12.5E %12.5E %12.5E %d %d\n",p.x(),p.y(),p.z(),(*it)->onWhat()->dim(),(*it)->onWhat()->tag());
(*it)->x() = p.x();
(*it)->y() = p.y();
(*it)->z() = p.z();
@@ -416,7 +403,7 @@ highOrderTools::highOrderTools (GModel *gm, GModel *mesh, int order)
computeStraightSidedPositions();
}
-void highOrderTools::computeStraightSidedPositions ()
+void highOrderTools::computeStraightSidedPositions()
{
_clean();
// compute straigh sided positions that are actually X,Y,Z positions
@@ -552,12 +539,12 @@ void highOrderTools::computeStraightSidedPositions ()
// apply a displacement that does not create elements that are
// distorted over a value "thres"
-double highOrderTools::apply_incremental_displacement (double max_incr,
- std::vector<MElement*> & v,
- bool mixed,
- double thres,
- char *meshName,
- std::vector<MElement*> & disto)
+double highOrderTools::apply_incremental_displacement(double max_incr,
+ std::vector<MElement*> & v,
+ bool mixed,
+ double thres,
+ char *meshName,
+ std::vector<MElement*> & disto)
{
#ifdef HAVE_PETSC
// assume that the mesh is OK, yet already curved
@@ -592,16 +579,11 @@ double highOrderTools::apply_incremental_displacement (double max_incr,
}
}
- printf("coucou1\n");
-
for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
MVertex *vert = *it;
std::map<MVertex*,SVector3>::iterator itt = _targetLocation.find(vert);
// impose displacement @ boundary
- // if (!(vert->onWhat()->geomType() == GEntity::Line
- // && vert->onWhat()->tag() == 2)){
if (itt != _targetLocation.end() && vert->onWhat()->dim() < _dim){
- // printf("fixing motion of vertex %d to %g %g\n",vert->getNum(),itt->second.x()-vert->x(),itt->second.y()-vert->y());
myAssembler.fixVertex(vert, 0, _tag, itt->second.x()-vert->x());
myAssembler.fixVertex(vert, 1, _tag, itt->second.y()-vert->y());
myAssembler.fixVertex(vert, 2, _tag, itt->second.z()-vert->z());
@@ -624,28 +606,21 @@ double highOrderTools::apply_incremental_displacement (double max_incr,
}
}
- printf("coucou2\n");
-
- //+++++++++ Assembly & Solve ++++++++++++++++++++++++++++++++++++
if (myAssembler.sizeOfR()){
// assembly of the elasticity term on the
- double t1 = Cpu();
for (unsigned int i = 0; i < v.size(); i++){
SElement se(v[i]);
- if (mixed)El_mixed.addToMatrix(myAssembler, &se);
+ if (mixed) El_mixed.addToMatrix(myAssembler, &se);
else El.addToMatrix(myAssembler, &se);
}
- double t2 = Cpu();
// solve the system
- printf("coucou3 %12.5E\n", t2-t1);
lsys->systemSolve();
- double t3 = Cpu();
- printf("coucou4 %12.5E\n", t3-t2);
}
- //+++++++++ Move vertices @ maximum ++++++++++++++++++++++++++++++++++++
+ // Move vertices @ maximum
FILE *fd = Fopen ("d.msh","w");
- fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$NodeData\n1\n\"tr(sigma)\"\n1\n0.0\n3\n1\n3\n%d\n", (int) _vertices.size());
+ fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$NodeData\n1\n"
+ "\"tr(sigma)\"\n1\n0.0\n3\n1\n3\n%d\n", (int) _vertices.size());
for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
double ax, ay, az;
myAssembler.getDofValue(*it, 0, _tag, ax);
@@ -659,12 +634,10 @@ double highOrderTools::apply_incremental_displacement (double max_incr,
fprintf(fd,"$EndNodeData\n");
fclose(fd);
- //+------------------- Check now if elements are ok ---------------+++++++
+ // Check now if elements are ok
(*_vertices.begin())->onWhat()->model()->writeMSH(meshName);
- printf("coucou5\n");
-
double percentage = max_incr * 100.;
while(1){
std::vector<MElement*> disto;
@@ -684,7 +657,6 @@ double highOrderTools::apply_incremental_displacement (double max_incr,
}
else break;
}
- printf("coucou6\n");
delete lsys;
return percentage;
@@ -708,8 +680,8 @@ void highOrderTools::ensureMinimumDistorsion(std::vector<MElement*> &all,
}
}
-double highOrderTools::applySmoothingTo (std::vector<MElement*> &all,
- double threshold, bool mixed)
+double highOrderTools::applySmoothingTo(std::vector<MElement*> &all,
+ double threshold, bool mixed)
{
int ITER = 0;
double minD;
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index b7f0dcb4fd..6e0e82342c 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -756,9 +756,6 @@ static void printLevel_onlysmall(const char* fn,
multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
std::map<MVertex*, SPoint3> &allCoordinates)
{
- //To go through this execute gmsh with the option -optimize_hom
- //if (!CTX::instance()->mesh.smoothInternalEdges)return;
-
//Find the boundary loops
//The loop with the largest equivalent radius is the Dirichlet boundary
std::vector<std::pair<MVertex*,double> > boundaryNodes;
diff --git a/Mesh/highOrderSmoother.cpp b/contrib/HighOrderMeshOptimizer/OptHomElastic.cpp
similarity index 100%
rename from Mesh/highOrderSmoother.cpp
rename to contrib/HighOrderMeshOptimizer/OptHomElastic.cpp
diff --git a/Mesh/highOrderSmoother.h b/contrib/HighOrderMeshOptimizer/OptHomElastic.h
similarity index 100%
rename from Mesh/highOrderSmoother.h
rename to contrib/HighOrderMeshOptimizer/OptHomElastic.h
diff --git a/doc/VERSIONS.txt b/doc/VERSIONS.txt
index 3035f508ec..284cbf082d 100644
--- a/doc/VERSIONS.txt
+++ b/doc/VERSIONS.txt
@@ -1,3 +1,5 @@
+2.8.2 (July x, 2013): improved high order tools.
+
2.8.1 (July 11, 2013): improved compound surfaces and transfinite arrangements.
2.8.0 (July 8, 2013): improved Delaunay point insertion; fixed mesh orientation
--
GitLab