From f63c2cb590a27bf8b68976f19760521daa0852f9 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 22 Mar 2015 12:13:36 +0000
Subject: [PATCH] cleanup
---
Mesh/pointInsertion.cpp | 309 +++++++++-----------
Mesh/pointInsertion.h | 2 +-
Mesh/pointInsertionRTreeTools.cpp | 113 +++----
Mesh/pointInsertionRTreeTools.h | 469 +++++++++++-------------------
Mesh/surfaceFiller.h | 11 +-
Mesh/yamakawa.cpp | 11 +-
Mesh/yamakawa.h | 3 +-
7 files changed, 366 insertions(+), 552 deletions(-)
diff --git a/Mesh/pointInsertion.cpp b/Mesh/pointInsertion.cpp
index 9fbcab3656..b6945fb3bb 100644
--- a/Mesh/pointInsertion.cpp
+++ b/Mesh/pointInsertion.cpp
@@ -1,3 +1,10 @@
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
+
#include <iostream>
#include <fstream>
#include <sstream>
@@ -38,26 +45,31 @@ void print_nodal_info(string filename, map<MVertex*, T> &mapp)
out.close();
}
-
-bool shoot(const SPoint2 &start, const SPoint2 &dir, const double &h, SPoint2 &res){
-
+bool shoot(const SPoint2 &start, const SPoint2 &dir,
+ const double &h, SPoint2 &res)
+{
const int RK = 1;
if (RK==1){
res = start + (dir*h);
- // cout << "(" << start[0] << "," <<start[1] << ") -> (" << res[0] << "," <<res[1] << ") " << endl;
+ // cout << "(" << start[0] << "," <<start[1] << ") -> (" << res[0] << ","
+ // <<res[1] << ") " << endl;
return true;
}
return false;
}
-
-bool computeFourNeighbors (frameFieldBackgroundMesh2D *bgm, MVertex *v_center, // the vertex for which we want to generate 4 neighbors (real vertex (xyz), not parametric !!! )
- SPoint2 &midpoint,
- bool goNonLinear, // do we compute the position in the real surface which is nonlinear
- SPoint2 newP[4][NUMDIR], // look into other directions
- SMetric3 &metricField) // the mesh metric
+bool computeFourNeighbors (frameFieldBackgroundMesh2D *bgm,
+ MVertex *v_center, // the vertex for which we want to
+ // generate 4 neighbors (real
+ // vertex (xyz), not parametric!)
+ SPoint2 &midpoint,
+ bool goNonLinear, // do we compute the position in
+ // the real surface which is
+ // nonlinear
+ SPoint2 newP[4][NUMDIR], // look into other directions
+ SMetric3 &metricField) // the mesh metric
{
// we assume that v is on surface gf, and backgroundMesh2D has been created based on gf
@@ -69,42 +81,38 @@ bool computeFourNeighbors (frameFieldBackgroundMesh2D *bgm, MVertex *v_center, /
// get RK info on midpoint (infos in two directions...)
RK_form infos;
- bgm->compute_RK_infos(midpoint[0],midpoint[1],v_center->x(), v_center->y(),v_center->z(),infos);
+ bgm->compute_RK_infos(midpoint[0],midpoint[1],v_center->x(), v_center->y(),
+ v_center->z(), infos);
metricField = infos.metricField;
// shoot in four directions
SPoint2 param_vec;
double h;
- // cout << "shooting..." << endl;
for (int i=0;i<4;i++){// in four directions
switch (i){
- case 0:
- param_vec = infos.paramt1;
- h = infos.paramh.first;
- break;
- case 1:
- param_vec = infos.paramt2;
- h = infos.paramh.second;
- break;
- case 2:
- param_vec = infos.paramt1 * -1.;
- h = infos.paramh.first;
- break;
- case 3:
- param_vec = infos.paramt2 * -1.;
- h = infos.paramh.second;
- break;
+ case 0:
+ param_vec = infos.paramt1;
+ h = infos.paramh.first;
+ break;
+ case 1:
+ param_vec = infos.paramt2;
+ h = infos.paramh.second;
+ break;
+ case 2:
+ param_vec = infos.paramt1 * -1.;
+ h = infos.paramh.first;
+ break;
+ case 3:
+ param_vec = infos.paramt2 * -1.;
+ h = infos.paramh.second;
+ break;
}
shoot(midpoint,param_vec,h,newP[i][0]);
- // cout << "(" << midpoint[0] << "," <<midpoint[1] << ") -> (" << newP[i][0][0] << "," << newP[i][0][1] << ") " << endl;
+ // cout << "(" << midpoint[0] << "," <<midpoint[1] << ") -> (" <<
+ // newP[i][0][0] << "," << newP[i][0][1] << ") " << endl;
}
-
-
- // -------------------------------------------------
// the following comes from surfaceFiller.cpp...
- // -------------------------------------------------
-
const double EPS = 1.e-7;
for (int j=0;j<2;j++){
for (int i=0;i<4;i++){
@@ -112,8 +120,8 @@ bool computeFourNeighbors (frameFieldBackgroundMesh2D *bgm, MVertex *v_center, /
}
}
- // We could stop here. Yet, if the metric varies a lot, we can solve
- // a nonlinear problem in order to find a better approximation in the real
+ // We could stop here. Yet, if the metric varies a lot, we can solve a
+ // nonlinear problem in order to find a better approximation in the real
// surface
if (1 && goNonLinear){
double L = infos.localsize;
@@ -131,8 +139,8 @@ bool computeFourNeighbors (frameFieldBackgroundMesh2D *bgm, MVertex *v_center, /
// if (newPoint[i][1] > rangeV.high())newPoint[i][1] = rangeV.high();
GPoint pp = gf->point(newP[i][0]);
double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
- (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
- (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
+ (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
+ (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
ERR[i] = 100*fabs(D-L)/(D+L);
// printf("L = %12.5E D = %12.5E ERR = %12.5E\n",L,D,100*fabs(D-L)/(D+L));
}
@@ -142,46 +150,51 @@ bool computeFourNeighbors (frameFieldBackgroundMesh2D *bgm, MVertex *v_center, /
for (int i=0;i<4;i++){
if (ERR[i] > 12){
double uvt[3] = {newPoint[i][0],newPoint[i][1],0.0};
- // printf("Intersecting with circle N = %g %g %g dir = %g %g %g R = %g p = %g %g %g\n",n.x(),n.y(),n.z(),dirs[i].x(),dirs[i].y(),dirs[i].z(),L,v_center->x(),v_center->y(),v_center->z());
- curveFunctorCircle cf (dirs[i],infos.normal, SVector3(v_center->x(),v_center->y(),v_center->z()), L);
+ // printf("Intersecting with circle N = %g %g %g dir = %g %g %g R
+ // = %g p = %g %g
+ // %g\n",n.x(),n.y(),n.z(),dirs[i].x(),dirs[i].y(),dirs[i].z(),L,
+ // v_center->x(),v_center->y(),v_center->z());
+ curveFunctorCircle cf (dirs[i],infos.normal,
+ SVector3(v_center->x(),v_center->y(),v_center->z()), L);
if (intersectCurveSurface (cf,ss,uvt,infos.paramh.first*1.e-3)){
GPoint pp = gf->point(SPoint2(uvt[0],uvt[1]));
double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
- (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
- (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
+ (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
+ (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
double DP = sqrt ((newPoint[i][0]-uvt[0])*(newPoint[i][0]-uvt[0]) +
- (newPoint[i][1]-uvt[1])*(newPoint[i][1]-uvt[1]));
+ (newPoint[i][1]-uvt[1])*(newPoint[i][1]-uvt[1]));
double newErr = 100*fabs(D-L)/(D+L);
// if (v_center->onWhat() != gf && gf->tag() == 3){
// crossField2d::normalizeAngle (uvt[2]);
// printf("INTERSECT angle = %g DP %g\n",uvt[2],DP);
// }
if (newErr < 1 && DP < .1){
- // printf("%12.5E vs %12.5E : %12.5E %12.5E vs %12.5E %12.5E \n",ERR[i],newErr,newPoint[i][0],newPoint[i][1],uvt[0],uvt[1]);
- newPoint[i][0] = uvt[0]; //
- newPoint[i][1] = uvt[1]; //
- } //
+ // printf("%12.5E vs %12.5E : %12.5E %12.5E vs %12.5E %12.5E
+ // \n",ERR[i],newErr,newPoint[i][0],newPoint[i][1],uvt[0],uvt[1]);
+ newPoint[i][0] = uvt[0];
+ newPoint[i][1] = uvt[1];
+ }
// printf("OK\n");
}
else{
Msg::Debug("Cannot put a new point on Surface %d",gf->tag());
- // printf("NOT OK\n");
+ // printf("NOT OK\n");
}
}
- } //
+ }
// return the four new vertices
for (int i=0;i<4;i++){
newP[i][0] = SPoint2(newPoint[i][0],newPoint[i][1]);
}
- } /// end non linear -------------------------------------------------//
+ }
return true;
}
-
-void computeTwoNeighbors(frameFieldBackgroundMesh3D *bgm, MVertex *parent, vector<MVertex*> &spawns, SVector3 dir, double h){
-
+void computeTwoNeighbors(frameFieldBackgroundMesh3D *bgm, MVertex *parent,
+ vector<MVertex*> &spawns, SVector3 dir, double h)
+{
// using approximate size, RK1...
double x = parent->x();
double y = parent->y();
@@ -200,10 +213,9 @@ void computeTwoNeighbors(frameFieldBackgroundMesh3D *bgm, MVertex *parent, vecto
spawns[1] = new MVertex(newx,newy,newz,gr,0);
}
-
-
-void computeSixNeighbors(frameFieldBackgroundMesh3D *bgm, MVertex *parent, vector<MVertex*> &spawns, STensor3 dir, double h){
-
+void computeSixNeighbors(frameFieldBackgroundMesh3D *bgm, MVertex *parent,
+ vector<MVertex*> &spawns, STensor3 dir, double h)
+{
// using approximate size, RK1...
double x = parent->x();
double y = parent->y();
@@ -224,25 +236,26 @@ void computeSixNeighbors(frameFieldBackgroundMesh3D *bgm, MVertex *parent, vecto
}
}
-
double Filler2D::time_bgm_and_smoothing = 0.;
double Filler2D::time_insertion = 0.;
-
Filler2D::Filler2D(){}
-
-Filler2D::~Filler2D(){
+Filler2D::~Filler2D()
+{
cout << "FILLER2D timing:" << endl;
- cout << " ------- CUMULATIVE TIME2D bgm & smoothing : " << time_bgm_and_smoothing << " s." << endl;
- cout << " ------- CUMULATIVE TIME2D inserting points : " << time_insertion << " s." << endl;
- cout << " ------- TOTAL 2D TIME (new) : " << time_bgm_and_smoothing+time_insertion << " s." << endl;
+ cout << " ------- CUMULATIVE TIME2D bgm & smoothing : "
+ << time_bgm_and_smoothing << " s." << endl;
+ cout << " ------- CUMULATIVE TIME2D inserting points : "
+ << time_insertion << " s." << endl;
+ cout << " ------- TOTAL 2D TIME (new) : "
+ << time_bgm_and_smoothing+time_insertion << " s." << endl;
}
-
-void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMetric3> &metrics){
- // NB/ do not use the mesh in GFace, use the one in backgroundMesh2D !!!
-
+void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed,
+ vector<SMetric3> &metrics)
+{
+ // NB/ do not use the mesh in GFace, use the one in backgroundMesh2D!
// if(debug) cout << " ------------------ OLD -------------------" << endl;
// stringstream ssa;
@@ -252,21 +265,19 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
// backgroundMesh::current()->print(ssa.str(),gf,0);
//
//
- //
- //
+ //
+ //
// if(debug) cout << " ------------------ NEW -------------------" << endl;
// backgroundMesh2D *bgm2 = dynamic_cast<backgroundMesh2D*>(BGMManager::get(gf));
// stringstream ss2;
// ss2 << "basebg_sizefield_" << gf->tag() << ".pos";
// bgm2->exportSizeField(ss2.str());
- //
- //
- //
+ //
+ //
+ //
// return;
//
-
-
BGMManager::set_use_cross_field(true);
const bool goNonLinear = true;
@@ -274,13 +285,14 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
const bool export_stuff=true;
if (debug) cout << "ENTERING POINTINSERTION2D" << endl;
-
+
double a;
// acquire background mesh
if(debug) cout << "pointInsertion2D: recover BGM" << endl;
a=Cpu();
- frameFieldBackgroundMesh2D *bgm = dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(gf));
+ frameFieldBackgroundMesh2D *bgm =
+ dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(gf));
time_bgm_and_smoothing += (Cpu() - a);
if (!bgm){
@@ -319,11 +331,10 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
if(debug) cout << "pointInsertion2D : get bnd vertices " << endl;
set<MVertex*> bnd_vertices = bgm->get_vertices_of_maximum_dim(1);
- // put boundary vertices in a fifo queue
+ // put boundary vertices in a fifo queue
set<smoothness_point_pair, compareSurfacePointWithExclusionRegionPtr_Smoothness> fifo;
vector<surfacePointWithExclusionRegion*> vertices;
-
// initiate the rtree
if(debug) cout << "pointInsertion2D : initiate RTree " << endl;
RTree<surfacePointWithExclusionRegion*,double,2,double> rtree;
@@ -334,14 +345,8 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
for (; it != bnd_vertices.end() ; ++it){
SPoint2 midpoint;
computeFourNeighbors(bgm,*it, midpoint, goNonLinear, newp, metricField);
- surfacePointWithExclusionRegion *sp = new surfacePointWithExclusionRegion (*it, newp, midpoint,metricField);
-
- // if(debug){
- // cout << " treating parent uv=(" << midpoint[0] << "," << midpoint[1] << ")" << endl;
- // for (int ii=0;ii<4;ii++)
- // cout << " child is uv=(" << newp[ii][0][0] << "," << newp[ii][0][1] << ")" << endl;
- // }
-
+ surfacePointWithExclusionRegion *sp = new surfacePointWithExclusionRegion
+ (*it, newp, midpoint,metricField);
smoothness_point_pair mp;
mp.ptr = sp;
@@ -354,9 +359,6 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
rtree.Insert(_min,_max,sp);
}
-
-
-
// ---------- main loop -----------------
while(!fifo.empty()){
if(debug) cout << " -------- fifo.size() = " << fifo.size() << endl;
@@ -365,31 +367,20 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
surfacePointWithExclusionRegion * parent = (*fifo.begin()).ptr;
fifo.erase(fifo.begin());
- // if(debug) cout << " treating parent (" << parent->_v->x() << "," << parent->_v->y() << "," << parent->_v->z() << ") uv=(" << parent->_center[0] << "," << parent->_center[1] << ")" << endl;
-
for (int dir=0;dir<NUMDIR;dir++){
for (int i=0;i<4;i++){
- //if(debug) cout << "pointInsertion2D main loop: check child " << i << endl;
- // if(debug) cout << " treating child (" << gp.x() << "," << gp.y() << "," << gp.z() << ") uv=(" << gp.u() << "," << gp.v() << ")" << endl;
- //// if(debug) cout << " child is uv=(" << parent->_p[i][dir][0] << "," << parent->_p[i][dir][1] << ")" << endl;
-
if (!inExclusionZone (parent->_p[i][dir], rtree, vertices) ){
GPoint gp = gf->point(parent->_p[i][dir]);
-
- //if(debug) cout << "pointInsertion2D main loop: one child not in exclusion zone "<< endl;
MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
SPoint2 midpoint;
- //if(debug) cout << "pointInsertion2D main loop: computeFourNeighbors " << endl;
computeFourNeighbors(bgm,v, midpoint, goNonLinear, newp, metricField);
- //if(debug) cout << "pointInsertion2D main loop: create new surfacePointWithExclusionRegion " << endl;
- surfacePointWithExclusionRegion *sp = new surfacePointWithExclusionRegion (v, newp, midpoint, metricField, parent);
- //if(debug) cout << "pointInsertion2D main loop: get smoothness " << endl;
+ surfacePointWithExclusionRegion *sp = new surfacePointWithExclusionRegion
+ (v, newp, midpoint, metricField, parent);
smoothness_point_pair mp;mp.ptr = sp;mp.rank=(1.-bgm->get_smoothness(gp.u(),gp.v()));
if (debug) vert_priority[v] = priority_counter++;
- //if(debug) cout << "pointInsertion2D main loop: insert in RTree " << endl;
fifo.insert(mp);
vertices.push_back(sp);
double _min[2],_max[2];
@@ -397,11 +388,12 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
rtree.Insert(_min,_max,sp);
if (debug){
- cout << " adding node (" << sp->_v->x() << "," << sp->_v->y() << "," << sp->_v->z() << ")" << endl;
- cout << " ----------------------------- sub --- fifo.size() = " << fifo.size() << endl;
+ cout << " adding node (" << sp->_v->x() << "," << sp->_v->y()
+ << "," << sp->_v->z() << ")" << endl;
+ cout << " ----------------------------- sub --- fifo.size() = "
+ << fifo.size() << endl;
}
count_nbaddedpt++;
-
}
}
}
@@ -409,7 +401,6 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
}
time_insertion += (Cpu() - a);
-
if (debug){
stringstream ss;
ss << "priority_" << gf->tag() << ".pos";
@@ -417,7 +408,6 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
ss.clear();
}
-
// add the vertices as additional vertices in the
// surface mesh
char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
@@ -435,47 +425,24 @@ void Filler2D::pointInsertion2D(GFace* gf, vector<MVertex*> &packed, vector<SMe
}
fprintf(f,"};");
fclose(f);
-
-
-
-
-
}
-
-bool Filler3D::treat_region(GRegion *gr){
-
- //if (gr->tag()>=11) return false;
- // if (gr->tag()==11) return false;
- // if (gr->tag()==12) return false;
- // if (gr->tag()==13) return false;
- // if (gr->tag()==17) return false;
-
-
-
+bool Filler3D::treat_region(GRegion *gr)
+{
BGMManager::set_use_cross_field(true);
bool use_vectorial_smoothness;
bool use_fifo;
string algo;
-// readValue("param.dat","SMOOTHNESSALGO",algo);
+ // readValue("param.dat","SMOOTHNESSALGO",algo);
algo.assign("SCALAR");
if (!algo.compare("SCALAR")){
- compareSmoothnessVertexPairs::vectorial = false;
use_vectorial_smoothness = false;
use_fifo = false;
}
- else if (!algo.compare("VECTORIAL")){
- compareSmoothnessVertexPairs::vectorial = true;
- use_vectorial_smoothness = true;
- use_fifo = false;
- cout << "Vectorial smoothness algo: abandonné pour l'instant !" << endl;
- throw;
- }
else if (!algo.compare("FIFO")){
- compareSmoothnessVertexPairs::vectorial = false;
use_vectorial_smoothness = false;
use_fifo = true;
}
@@ -490,11 +457,11 @@ bool Filler3D::treat_region(GRegion *gr){
cout << "ENTERING POINTINSERTION3D" << endl;
-
// acquire background mesh
cout << "pointInsertion3D: recover BGM" << endl;
a = Cpu();
- frameFieldBackgroundMesh3D *bgm = dynamic_cast<frameFieldBackgroundMesh3D*>(BGMManager::get(gr));
+ frameFieldBackgroundMesh3D *bgm =
+ dynamic_cast<frameFieldBackgroundMesh3D*>(BGMManager::get(gr));
time_smoothing += (Cpu() - a);
if (!bgm){
@@ -502,7 +469,6 @@ bool Filler3D::treat_region(GRegion *gr){
throw;
}
-
// export BGM fields
if(export_stuff){
cout << "pointInsertion3D: export size field " << endl;
@@ -528,8 +494,6 @@ bool Filler3D::treat_region(GRegion *gr){
}
}
-
-
// ---------------- START FILLING NEW POINTS ----------------
cout << "pointInsertion3D : inserting points in region " << gr->tag() << endl;
@@ -553,16 +517,16 @@ bool Filler3D::treat_region(GRegion *gr){
map<MVertex*,double> smoothness_forplot;
MElement *element;
MVertex *vertex;
- list<GFace*> faces = gr->faces();
+ list<GFace*> faces = gr->faces();
for(list<GFace*>::iterator it=faces.begin();it!=faces.end();it++){// for all faces
GFace *gf = *it;
- // int limit = code_kesskessai(gf->tag());
+ // int limit = code_kesskessai(gf->tag());
for(unsigned int i=0;i<gf->getNumMeshElements();i++){
element = gf->getMeshElement(i);
for(int j=0;j<element->getNumVertices();j++){// for all vertices
vertex = element->getVertex(j);
temp.insert(vertex);
- // limits.insert(make_pair(vertex,limit));
+ // limits.insert(make_pair(vertex,limit));
}
}
}
@@ -580,7 +544,8 @@ bool Filler3D::treat_region(GRegion *gr){
y = boundary_vertices[i]->y();
z = boundary_vertices[i]->z();
- MVertex *closest = bgm->get_nearest_neighbor_on_boundary(boundary_vertices[i]);// "on boundary since working on boundary_vertices ...
+ // "on boundary since working on boundary_vertices ...
+ MVertex *closest = bgm->get_nearest_neighbor_on_boundary(boundary_vertices[i]);
h = bgm->size(closest);// get approximate size, closest vertex, faster ?!
fill_min_max(x,y,z,h,min,max);
@@ -592,7 +557,7 @@ bool Filler3D::treat_region(GRegion *gr){
svp->v = boundary_vertices[i];
svp->rank = bgm->get_smoothness(x,y,z);
svp->dir = 0;
- svp->layer = 0;
+ svp->layer = 0;
svp->size = h;
bgm->eval_approximate_crossfield(closest, svp->cf);
@@ -609,21 +574,19 @@ bool Filler3D::treat_region(GRegion *gr){
svp->v = boundary_vertices[i];
svp->rank = bgm->get_vectorial_smoothness(idir,x,y,z);
svp->dir = idir;
- svp->layer = 0;
+ svp->layer = 0;
svp->size = h;
svp->cf = temp;
for (int k=0;k<3;k++) svp->direction(k) = temp(k,idir);
- // cout << "fifo size=" << fifo->size() << " inserting " ;
+ // cout << "fifo size=" << fifo->size() << " inserting " ;
fifo->insert(svp);
- // cout << " -> fifo size=" << fifo->size() << endl;
+ // cout << " -> fifo size=" << fifo->size() << endl;
}
}
}
- // TODO: si fifo était list of *PTR -> pas de copies, gain temps ?
-
-
+ // TODO: si fifo était list of *PTR -> pas de copies, gain temps ?
Wrapper3D wrapper;
wrapper.set_bgm(bgm);
MVertex *parent,*individual;
@@ -646,15 +609,19 @@ bool Filler3D::treat_region(GRegion *gr){
vector<MVertex*> spawns;
if (!use_vectorial_smoothness){
spawns.resize(6);
- computeSixNeighbors(bgm,parent,spawns,fifo->get_first_crossfield(),fifo->get_first_size());
+ computeSixNeighbors(bgm,parent,spawns,fifo->get_first_crossfield(),
+ fifo->get_first_size());
}
else{
spawns.resize(2);
- computeTwoNeighbors(bgm,parent,spawns,fifo->get_first_direction(),fifo->get_first_size());
+ computeTwoNeighbors(bgm,parent,spawns,fifo->get_first_direction(),
+ fifo->get_first_size());
}
fifo->erase_first();
- // cout << "while, fifo->size()=" << fifo->size() << " parent=(" << parent->x() << "," << parent->y() << "," << parent->z() << ")" << endl;
+ // cout << "while, fifo->size()=" << fifo->size() << " parent=(" <<
+ // parent->x() << "," << parent->y() << "," << parent->z() << ")" <<
+ // endl;
for(unsigned int i=0;i<spawns.size();i++){
spawn_created = false;
@@ -662,7 +629,8 @@ bool Filler3D::treat_region(GRegion *gr){
x = individual->x();
y = individual->y();
z = individual->z();
- // cout << " working on candidate " << "(" << individual->x() << "," << individual->y() << "," << individual->z() << ")" << endl;
+ // cout << " working on candidate " << "(" << individual->x() << ","
+ // << individual->y() << "," << individual->z() << ")" << endl;
if(bgm->inDomain(x,y,z)){
// cout << " spawn " << i << " in domain" << endl;
@@ -691,7 +659,7 @@ bool Filler3D::treat_region(GRegion *gr){
svp->v = individual;
svp->rank=bgm->get_smoothness(individual->x(),individual->y(),individual->z());
svp->dir = 0;
- svp->layer = parent_layer+1;
+ svp->layer = parent_layer+1;
svp->size = h;
svp->cf = crossfield;
fifo->insert(svp);
@@ -708,7 +676,7 @@ bool Filler3D::treat_region(GRegion *gr){
svp->v = individual;
svp->rank = bgm->get_vectorial_smoothness(idir,x,y,z);
svp->dir = idir;
- svp->layer = parent_layer+1;
+ svp->layer = parent_layer+1;
svp->size = h;
for (int k=0;k<3;k++) svp->direction(k) = crossfield(k,idir);
svp->cf = crossfield;
@@ -733,7 +701,6 @@ bool Filler3D::treat_region(GRegion *gr){
time_insert_points += (Cpu() - a);
-
// --- output ---
if (debug){
stringstream ss;
@@ -755,7 +722,6 @@ bool Filler3D::treat_region(GRegion *gr){
int option = CTX::instance()->mesh.algo3d;
CTX::instance()->mesh.algo3d = ALGO_3D_DELAUNAY;
-
deMeshGRegion deleter;
deleter(gr);
std::vector<GRegion*> regions;
@@ -775,31 +741,32 @@ bool Filler3D::treat_region(GRegion *gr){
new_vertices.clear();
rtree.RemoveAll();
-
-
return true;
}
-
-int Filler3D::get_nbr_new_vertices(){
+int Filler3D::get_nbr_new_vertices()
+{
return new_vertices.size();
}
-
-MVertex* Filler3D::get_new_vertex(int i){
+MVertex* Filler3D::get_new_vertex(int i)
+{
return new_vertices[i];
}
-
Filler3D::Filler3D(){}
-
-Filler3D::~Filler3D(){
+Filler3D::~Filler3D()
+{
cout << "FILLER3D timing:" << endl;
- cout << " ------- CUMULATIVE TIME3D bgm & smoothing : " << time_smoothing << " s." << endl;
- cout << " ------- CUMULATIVE TIME3D inserting points : " << time_insert_points << " s." << endl;
- cout << " ------- CUMULATIVE TIME3D meshing region : " << time_meshing << " s." << endl;
- cout << " ------- CUMULATIVE TOTAL 3D TIME (new) : " << time_meshing+time_smoothing+time_insert_points << " s." << endl;
+ cout << " ------- CUMULATIVE TIME3D bgm & smoothing : "
+ << time_smoothing << " s." << endl;
+ cout << " ------- CUMULATIVE TIME3D inserting points : "
+ << time_insert_points << " s." << endl;
+ cout << " ------- CUMULATIVE TIME3D meshing region : "
+ << time_meshing << " s." << endl;
+ cout << " ------- CUMULATIVE TOTAL 3D TIME (new) : "
+ << time_meshing+time_smoothing+time_insert_points << " s." << endl;
}
std::vector<MVertex*> Filler3D::new_vertices;
diff --git a/Mesh/pointInsertion.h b/Mesh/pointInsertion.h
index ecda4439e7..156f0a8b25 100644
--- a/Mesh/pointInsertion.h
+++ b/Mesh/pointInsertion.h
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
diff --git a/Mesh/pointInsertionRTreeTools.cpp b/Mesh/pointInsertionRTreeTools.cpp
index d8845e1ad1..debf01c640 100644
--- a/Mesh/pointInsertionRTreeTools.cpp
+++ b/Mesh/pointInsertionRTreeTools.cpp
@@ -1,15 +1,20 @@
-#include "pointInsertionRTreeTools.h"
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Tristan Carrier and Paul-Emile Bernard
+#include "pointInsertionRTreeTools.h"
#include "BackgroundMesh.h"
#include "BackgroundMeshManager.h"
-
#include "pointInsertion.h"
-
#include "GEntity.h"
-
-
-surfacePointWithExclusionRegion::surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], SPoint2 &_mp, SMetric3 & meshMetric, surfacePointWithExclusionRegion *father){
+surfacePointWithExclusionRegion::surfacePointWithExclusionRegion
+ (MVertex *v, SPoint2 p[4][NUMDIR], SPoint2 &_mp, SMetric3 & meshMetric,
+ surfacePointWithExclusionRegion *father)
+{
_v = v;
_meshMetric = meshMetric;
_center = _mp;
@@ -28,8 +33,8 @@ surfacePointWithExclusionRegion::surfacePointWithExclusionRegion (MVertex *v, SP
}
}
-
-bool surfacePointWithExclusionRegion::inExclusionZone (const SPoint2 &p){
+bool surfacePointWithExclusionRegion::inExclusionZone (const SPoint2 &p)
+{
double mat[2][2];
double b[2] , uv[2];
mat[0][0]= _q[1].x()-_q[0].x();
@@ -53,16 +58,16 @@ bool surfacePointWithExclusionRegion::inExclusionZone (const SPoint2 &p){
return false;
}
-
-void surfacePointWithExclusionRegion::minmax (double _min[2], double _max[2]) const{
+void surfacePointWithExclusionRegion::minmax (double _min[2], double _max[2]) const
+{
_min[0] = std::min(std::min(std::min(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x());
_min[1] = std::min(std::min(std::min(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y());
_max[0] = std::max(std::max(std::max(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x());
_max[1] = std::max(std::max(std::max(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y());
}
-
-void surfacePointWithExclusionRegion::print (FILE *f, int i){
+void surfacePointWithExclusionRegion::print (FILE *f, int i)
+{
fprintf(f,"SP(%g,%g,%g){%d};\n",_center.x(),_center.y(),0.0,i);
fprintf(f,"SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
_q[0].x(),_q[0].y(),0.0,
@@ -72,15 +77,10 @@ void surfacePointWithExclusionRegion::print (FILE *f, int i){
}
-
-
my_wrapper::my_wrapper (SPoint2 sp) : _tooclose (false), _p(sp) {}
-
-
-
-
-bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){
+bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point)
+{
my_wrapper *w = static_cast<my_wrapper*>(point);
if (neighbour->inExclusionZone(w->_p)){
@@ -91,11 +91,11 @@ bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){
return true;
}
-
-
-bool inExclusionZone (SPoint2 &p,
- RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree,
- std::vector<surfacePointWithExclusionRegion*> & all ){
+bool inExclusionZone
+ (SPoint2 &p,
+ RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree,
+ std::vector<surfacePointWithExclusionRegion*> & all )
+{
// should assert that the point is inside the domain
// OLD BGM
if (old_algo_hexa()){
@@ -114,25 +114,18 @@ bool inExclusionZone (SPoint2 &p,
for (unsigned int i=0;i<all.size();++i){
if (all[i]->inExclusionZone(p)){
- // printf("%g %g is in exclusion zone of %g %g\n",p.x(),p.y(),all[i]._center.x(),all[i]._center.y());
+ // printf("%g %g is in exclusion zone of %g
+ // %g\n",p.x(),p.y(),all[i]._center.x(),all[i]._center.y());
return true;
}
}
return false;
}
-
-
-// ------------------------------------------------------------------------------------
-// --------------------------------- 3D --------------------------------------------
-// ------------------------------------------------------------------------------------
-
frameFieldBackgroundMesh3D* Wrapper3D::bgmesh = NULL;
-bool compareSmoothnessVertexPairs::vectorial = false;
-
-
-double infinity_distance_3D(const MVertex *v1,const MVertex *v2,STensor3 &cf){
+double infinity_distance_3D(const MVertex *v1,const MVertex *v2,STensor3 &cf)
+{
SPoint3 p1 = v1->point();
SPoint3 p2 = v2->point();
double x1=0.;
@@ -154,8 +147,8 @@ double infinity_distance_3D(const MVertex *v1,const MVertex *v2,STensor3 &cf){
return std::max(std::max(fabs(x2-x1),fabs(y2-y1)),fabs(z2-z1));// distance
};
-
-void fill_min_max(double x,double y,double z,double h,double *min,double *max){
+void fill_min_max(double x,double y,double z,double h,double *min,double *max)
+{
min[0] = x - sqrt3*h;
max[0] = x + sqrt3*h;
min[1] = y - sqrt3*h;
@@ -164,19 +157,18 @@ void fill_min_max(double x,double y,double z,double h,double *min,double *max){
max[2] = z + sqrt3*h;
};
-
-bool rtree_callback_3D(MVertex* neighbour,void* w){
+bool rtree_callback_3D(MVertex* neighbour,void* w)
+{
Wrapper3D* wrapper;
wrapper = static_cast<Wrapper3D*>(w);
const MVertex* individual = wrapper->get_individual();
const MVertex *parent = wrapper->get_parent();
if (parent==neighbour) return true;
-
-// frameFieldBackgroundMesh3D* bgm = wrapper->bgm();
-// const MVertex *closest = bgm->get_nearest_neighbor(individual);
-// const double h = bgm->size(closest);// get approximate size, closest vertex, faster ?!
-// STensor3 crossfield;
-// bgm->eval_approximate_crossfield(closest, crossfield);
+ // frameFieldBackgroundMesh3D* bgm = wrapper->bgm();
+ // const MVertex *closest = bgm->get_nearest_neighbor(individual);
+ // const double h = bgm->size(closest);// get approximate size, closest vertex, faster ?!
+ // STensor3 crossfield;
+ // bgm->eval_approximate_crossfield(closest, crossfield);
double *h = wrapper->get_size();
STensor3 *crossfield = wrapper->get_crossfield();
@@ -188,10 +180,10 @@ bool rtree_callback_3D(MVertex* neighbour,void* w){
return true;
};
-
-bool far_from_boundary_3D(frameFieldBackgroundMesh3D *bgm, MVertex* v, double h){
+bool far_from_boundary_3D(frameFieldBackgroundMesh3D *bgm, MVertex* v, double h)
+{
// check if the box (v->point +- k2*h) is in domain
-
+
const double x = v->x();
const double y = v->y();
const double z = v->z();
@@ -205,28 +197,3 @@ bool far_from_boundary_3D(frameFieldBackgroundMesh3D *bgm, MVertex* v, double h)
return true;
};
-
-
-// vient de l'ancien code hexa, mais... kesskessai ?
-//int code_kesskessai(int tag){
-// int limit;
-// std::string s;
-// std::stringstream s2;
-//
-// limit = -1;
-// s2 << tag;
-// s = s2.str();
-//
-// if(s.length()>=5){
-// if(s.at(0)=='1' && s.at(1)=='1' && s.at(2)=='1' && s.at(3)=='1' && s.at(4)=='1'){
-// limit = 0;
-// }
-// else if(s.at(0)=='2' && s.at(1)=='2' && s.at(2)=='2' && s.at(3)=='2' && s.at(4)=='2'){
-// limit = 1;
-// }
-// }
-//
-// return limit;
-//}
-
-
diff --git a/Mesh/pointInsertionRTreeTools.h b/Mesh/pointInsertionRTreeTools.h
index 22cf1ffe6a..8ac501c870 100644
--- a/Mesh/pointInsertionRTreeTools.h
+++ b/Mesh/pointInsertionRTreeTools.h
@@ -1,21 +1,21 @@
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Tristan Carrier and Paul-Emile Bernard
-#ifndef _POINTINSERTIONTOOLS_H_
-#define _POINTINSERTIONTOOLS_H_
-
+#ifndef _POINTINSERTIONRTREETOOLS_H_
+#define _POINTINSERTIONRTREETOOLS_H_
+#include <math.h>
+#include <queue>
#include "MVertex.h"
#include "STensor3.h"
-
#include "BackgroundMesh3D.h"
#include "GEntity.h"
-
#include "rtree.h"
-#include <math.h>
-#include <queue>
-
-//using namespace std;
-
static const double k1 = 0.61; //k1*h is the minimal distance between two nodes
static const double k2 = 0.5; //k2*h is the minimal distance to the boundary
static const double sqrt3 = 1.73205081;
@@ -27,18 +27,14 @@ static const double DIRS [NUMDIR] = {0.0};
//static const int NUMDIR = 3;
//static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.};
-
-
-
-
class surfacePointWithExclusionRegion {
- public:
- MVertex *_v;
- SPoint2 _center;
- SPoint2 _p[4][NUMDIR];
- SPoint2 _q[4];
- SMetric3 _meshMetric;
- double _distanceSummed;
+public:
+ MVertex *_v;
+ SPoint2 _center;
+ SPoint2 _p[4][NUMDIR];
+ SPoint2 _q[4];
+ SMetric3 _meshMetric;
+ double _distanceSummed;
/*
+ p3
p4 |
@@ -47,327 +43,212 @@ class surfacePointWithExclusionRegion {
+ p1
*/
+ surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR],
+ SPoint2 &_mp, SMetric3 & meshMetric,
+ surfacePointWithExclusionRegion *father = 0);
- surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], SPoint2 &_mp, SMetric3 & meshMetric, surfacePointWithExclusionRegion *father = 0);
-
- bool inExclusionZone (const SPoint2 &p);
- void minmax (double _min[2], double _max[2]) const;
- void print (FILE *f, int i);
+ bool inExclusionZone (const SPoint2 &p);
+ void minmax (double _min[2], double _max[2]) const;
+ void print (FILE *f, int i);
};
-
-
class my_wrapper {
- public:
- bool _tooclose;
- SPoint2 _p;
- my_wrapper (SPoint2 sp);
+public:
+ bool _tooclose;
+ SPoint2 _p;
+ my_wrapper (SPoint2 sp);
};
-
-
struct smoothness_point_pair{
double rank;
surfacePointWithExclusionRegion* ptr;
};
-
-
class compareSurfacePointWithExclusionRegionPtr_Smoothness
{
- public:
- inline bool operator () (const smoothness_point_pair &a, const smoothness_point_pair &b) const
- {
- if (a.rank == b.rank){
- if(a.ptr->_distanceSummed > b.ptr->_distanceSummed) return false;
- if(a.ptr->_distanceSummed < b.ptr->_distanceSummed) return true;
- return a.ptr<b.ptr;
- }
- // else
- return (a.rank < b.rank);
+public:
+ inline bool operator () (const smoothness_point_pair &a,
+ const smoothness_point_pair &b) const
+ {
+ if (a.rank == b.rank){
+ if(a.ptr->_distanceSummed > b.ptr->_distanceSummed) return false;
+ if(a.ptr->_distanceSummed < b.ptr->_distanceSummed) return true;
+ return a.ptr<b.ptr;
}
+ return (a.rank < b.rank);
+ }
};
-
-
-
class compareSurfacePointWithExclusionRegionPtr
{
- public:
- inline bool operator () (const surfacePointWithExclusionRegion *a, const surfacePointWithExclusionRegion *b) const
- {
- if(a->_distanceSummed > b->_distanceSummed) return false;
- if(a->_distanceSummed < b->_distanceSummed) return true;
- return a<b;
- }
+public:
+ inline bool operator () (const surfacePointWithExclusionRegion *a,
+ const surfacePointWithExclusionRegion *b) const
+ {
+ if(a->_distanceSummed > b->_distanceSummed) return false;
+ if(a->_distanceSummed < b->_distanceSummed) return true;
+ return a<b;
+ }
};
-
-
-
extern bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point);
-
-
-
-extern bool inExclusionZone (SPoint2 &p, RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree, std::vector<surfacePointWithExclusionRegion*> & all );
-
-
-
-
-
-
-
-// ------------------------------------------------------------------------------------
-// --------------------------------- 3D --------------------------------------------
-// ------------------------------------------------------------------------------------
-
-
+extern bool inExclusionZone (SPoint2 &p,
+ RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree,
+ std::vector<surfacePointWithExclusionRegion*> & all);
class Wrapper3D{
- private:
- static frameFieldBackgroundMesh3D* bgmesh;
- MVertex *individual,*parent;
- double *size;
- STensor3 *cf;
- SVector3 *vec;
- bool ok;
- // Node* individual;
- // Node* parent;
- public:
- Wrapper3D():ok(true){};
- Wrapper3D(MVertex* _i,MVertex* _p):individual(_i), parent(_p),ok(true){};
- ~Wrapper3D(){};
- void set_ok(bool b){ok=b;};
- void set_individual(MVertex *vertex){individual=vertex;};
- void set_parent(MVertex *vertex){parent=vertex;};
- void set_size(double *h){size=h;};
- void set_crossfield(STensor3 *_cf){cf=_cf;};
- void set_direction(SVector3 *_v){vec=_v;};
- bool get_ok(){return ok;};
- void set_bgm(frameFieldBackgroundMesh3D *bgm){bgmesh = bgm;};
- frameFieldBackgroundMesh3D * bgm(){return bgmesh;};
- MVertex* get_individual(){return individual;};
- MVertex* get_parent(){return parent;};
- STensor3* get_crossfield(){return cf;};
- SVector3* get_direction(){return vec;};
- double* get_size(){return size;};
+private:
+ static frameFieldBackgroundMesh3D* bgmesh;
+ MVertex *individual,*parent;
+ double *size;
+ STensor3 *cf;
+ SVector3 *vec;
+ bool ok;
+public:
+ Wrapper3D():ok(true){};
+ Wrapper3D(MVertex* _i,MVertex* _p):individual(_i), parent(_p),ok(true){};
+ ~Wrapper3D(){};
+ void set_ok(bool b){ok=b;};
+ void set_individual(MVertex *vertex){individual=vertex;};
+ void set_parent(MVertex *vertex){parent=vertex;};
+ void set_size(double *h){size=h;};
+ void set_crossfield(STensor3 *_cf){cf=_cf;};
+ void set_direction(SVector3 *_v){vec=_v;};
+ bool get_ok(){return ok;};
+ void set_bgm(frameFieldBackgroundMesh3D *bgm){bgmesh = bgm;};
+ frameFieldBackgroundMesh3D * bgm(){return bgmesh;};
+ MVertex* get_individual(){return individual;};
+ MVertex* get_parent(){return parent;};
+ STensor3* get_crossfield(){return cf;};
+ SVector3* get_direction(){return vec;};
+ double* get_size(){return size;};
};
-
-
-
extern double infinity_distance_3D(const MVertex *v1,const MVertex *v2,STensor3 &cf);
extern bool rtree_callback_3D(MVertex* neighbour,void* w);
extern bool far_from_boundary_3D(frameFieldBackgroundMesh3D *bgm, MVertex* v, double h);
-//extern int code_kesskessai(int tag);
extern void fill_min_max(double x,double y,double z,double h,double *min,double *max);
-
-
-
-// TODO: this is not a pair anymore, the name has to be changed
-// this will be in listOfPoints AND !!! in RTree: larger memory footprint but less CPU time...
+// TODO: this is not a pair anymore, the name has to be changed; this will be in
+// listOfPoints AND in RTree: larger memory footprint but less CPU time...
class smoothness_vertex_pair{
- public:
- smoothness_vertex_pair(){};
-// smoothness_vertex_pair(const smoothness_vertex_pair &other){
-// rank = other.rank;
-// v = other.v;
-// dir = other.dir;
-// layer = other.layer;
-// };
- ~smoothness_vertex_pair(){};
- STensor3 cf;
- SVector3 direction;
- double rank,size;
- MVertex *v;
- int dir,layer;
+public:
+ smoothness_vertex_pair(){};
+ ~smoothness_vertex_pair(){};
+ STensor3 cf;
+ SVector3 direction;
+ double rank, size;
+ MVertex *v;
+ int dir, layer;
};
-
-
class compareSmoothnessVertexPairs
{
- private:
- const double threshold;
- double roundit(const double &d)const{
- //return (round(d/threshold)*threshold);
- return ((int)(d/threshold + 0.5)*threshold);
- }
- public:
- static bool vectorial;
-
- compareSmoothnessVertexPairs():threshold(0.05){};
-
- inline bool operator () (const smoothness_vertex_pair *a, const smoothness_vertex_pair *b) const
- {
-
- // THIS IS USED :
- if (!vectorial){
- if (a->rank==b->rank) return (a->v<b->v);
- return (a->rank > b->rank);
- }
-
- // THE FOLLOWING IS USUALLY NOT USED... IT WAS AN ATTEMPT TO USE VECTORIAL SMOOTHNESS, SHOULD BE DELETED...
-
-
- // cout << " operator() " << endl;
- if (roundit(a->rank)==roundit(b->rank) && ((a->rank+b->rank)/2.>=0.9)){// if good smoothness and approximately equal
- // cout << " smoothness same,";
- if (a->layer!=b->layer){// layers are different
- // cout << "layers different, return " << (a->layer > b->layer) << endl;
- return (a->layer > b->layer);// priority to larger layer
- }
- //cout << "a b rankls: " << roundit(a->rank) << "," << roundit(b->rank) << " layers=" << a->layer << "," << b->layer << " -> return " << (a->layer > b->layer) << endl;
-
- if ((a->layer==0) && (b->layer==0)){// both on boundaries
- // cout << "both on boundaries,";
- if ((a->v->onWhat()->dim()==2) && (b->v->onWhat()->dim()==2)){// both on faces
- // cout << "both on faces,";
- if ((a->v->onWhat()->tag()) != (b->v->onWhat()->tag())){// with different tags
-
-
- // cout << "with different tags, return " << (a->v->onWhat()->tag() < b->v->onWhat()->tag()) << endl;
- return (a->v->onWhat()->tag() < b->v->onWhat()->tag());// priority to smaller tag
- }
- }
- }
- }
-
- if (a->rank==b->rank){// same smoothness
- // cout << " same smoothness";
- if (a->v==b->v){// same vertex
- // cout << " same vertex, return " << (a->dir < b->dir) << endl;
- return (a->dir < b->dir);
- }
- // cout << "return " << (a->v < b->v) << endl;
- return (a->v<b->v);
- }
- // cout << " comparing smoothness, return " << (a->rank > b->rank) << endl;
- // cout << "no roundit !!!" << endl;
-
- return (a->rank > b->rank);
- }
+private:
+ const double threshold;
+ double roundit(const double &d) const
+ {
+ //return (round(d/threshold)*threshold);
+ return ((int)(d/threshold + 0.5)*threshold);
+ }
+public:
+ compareSmoothnessVertexPairs():threshold(0.05){};
+
+ inline bool operator () (const smoothness_vertex_pair *a,
+ const smoothness_vertex_pair *b) const
+ {
+ if (a->rank==b->rank) return (a->v<b->v);
+ return (a->rank > b->rank);
+ }
};
-
-
class listOfPoints{
- public:
- listOfPoints(){};
- virtual ~listOfPoints(){};
- virtual void insert(smoothness_vertex_pair *svp)=0;
- virtual unsigned int size()=0;
- virtual MVertex* get_first_vertex()=0;
- virtual STensor3 get_first_crossfield()=0;
- virtual double get_first_size()=0;
- virtual int get_first_layer()=0;
- virtual SVector3 get_first_direction()=0;
- virtual void erase_first()=0;
- virtual bool empty()=0;
+public:
+ listOfPoints(){};
+ virtual ~listOfPoints(){};
+ virtual void insert(smoothness_vertex_pair *svp)=0;
+ virtual unsigned int size()=0;
+ virtual MVertex* get_first_vertex()=0;
+ virtual STensor3 get_first_crossfield()=0;
+ virtual double get_first_size()=0;
+ virtual int get_first_layer()=0;
+ virtual SVector3 get_first_direction()=0;
+ virtual void erase_first()=0;
+ virtual bool empty()=0;
};
-
class listOfPointsScalarSmoothness : public listOfPoints{
- public:
- listOfPointsScalarSmoothness(){
- cout << "USING SMOOTHNESS-BASED LIST" << endl;
- };
- virtual ~listOfPointsScalarSmoothness(){
- while (!empty())
- erase_first();
- };
- virtual void insert(smoothness_vertex_pair *svp){points.insert(svp);};
- virtual unsigned int size(){return points.size();};
- virtual MVertex* get_first_vertex(){return (*points.begin())->v;};
- virtual STensor3 get_first_crossfield(){return (*points.begin())->cf;};
- virtual double get_first_size(){return (*points.begin())->size;};
- virtual int get_first_layer(){return (*points.begin())->layer;};
- virtual SVector3 get_first_direction(){
- cout << "listOfPointsScalarSmoothness:: get_first_direction NOT applicable ! " << endl;
- throw;
- return SVector3(0.);
- };
- virtual void erase_first(){
- smoothness_vertex_pair *ptr = *(points.begin());
- points.erase(points.begin());
- delete ptr;
- };
- virtual bool empty(){return points.empty();};
-
- protected:
- set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
+public:
+ listOfPointsScalarSmoothness(){ };
+ virtual ~listOfPointsScalarSmoothness()
+ {
+ while (!empty())
+ erase_first();
+ };
+ virtual void insert(smoothness_vertex_pair *svp){ points.insert(svp); };
+ virtual unsigned int size(){ return points.size(); };
+ virtual MVertex* get_first_vertex(){ return (*points.begin())->v; };
+ virtual STensor3 get_first_crossfield(){ return (*points.begin())->cf; };
+ virtual double get_first_size(){ return (*points.begin())->size; };
+ virtual int get_first_layer(){ return (*points.begin())->layer; };
+ virtual SVector3 get_first_direction()
+ {
+ Msg::Error("listOfPointsScalarSmoothness::get_first_direction NOT applicable");
+ return SVector3(0.);
+ };
+ virtual void erase_first()
+ {
+ smoothness_vertex_pair *ptr = *(points.begin());
+ points.erase(points.begin());
+ delete ptr;
+ };
+ virtual bool empty(){ return points.empty(); };
+
+protected:
+ set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
};
-
class listOfPointsVectorialSmoothness : public listOfPointsScalarSmoothness{
- public:
- listOfPointsVectorialSmoothness(){
- cout << "USING VECTORIAL SMOOTHNESS-BASED LIST" << endl;
- };
- virtual ~listOfPointsVectorialSmoothness(){
- while (!empty())
- erase_first();
- };
- virtual SVector3 get_first_direction(){return (*points.begin())->direction;};
-
- protected:
- set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
+public:
+ listOfPointsVectorialSmoothness(){};
+ virtual ~listOfPointsVectorialSmoothness(){
+ while (!empty())
+ erase_first();
+ };
+ virtual SVector3 get_first_direction(){ return (*points.begin())->direction; };
+protected:
+ set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
};
-
class listOfPointsFifo : public listOfPoints{
- public:
- listOfPointsFifo(){
- cout << "USING FIFO LIST" << endl;
- };
- virtual ~listOfPointsFifo(){
- while (!empty())
- erase_first();
- };
- virtual void insert(smoothness_vertex_pair *svp){points.push(svp);};
- virtual unsigned int size(){return points.size();};
- virtual MVertex* get_first_vertex(){return (points.front())->v;};
- virtual STensor3 get_first_crossfield(){return (points.front())->cf;};
- virtual double get_first_size(){return (points.front())->size;};
- virtual int get_first_layer(){return (points.front())->layer;};
- virtual SVector3 get_first_direction(){
- cout << "listOfPointsFifo:: get_first_direction NOT applicable ! " << endl;
- throw;
- return SVector3(0.);
- };
- virtual void erase_first(){
- smoothness_vertex_pair *ptr = points.front();
- points.pop();
- delete ptr;
- };
- virtual bool empty(){return points.empty();};
-
- protected:
- queue<smoothness_vertex_pair*> points;
+public:
+ listOfPointsFifo(){};
+ virtual ~listOfPointsFifo(){
+ while (!empty())
+ erase_first();
+ };
+ virtual void insert(smoothness_vertex_pair *svp){ points.push(svp); };
+ virtual unsigned int size(){ return points.size(); };
+ virtual MVertex* get_first_vertex(){ return (points.front())->v; };
+ virtual STensor3 get_first_crossfield(){ return (points.front())->cf; };
+ virtual double get_first_size(){ return (points.front())->size; };
+ virtual int get_first_layer(){ return (points.front())->layer; };
+ virtual SVector3 get_first_direction()
+ {
+ Msg::Error("listOfPointsFifo::get_first_direction NOT applicable");
+ return SVector3(0.);
+ };
+ virtual void erase_first()
+ {
+ smoothness_vertex_pair *ptr = points.front();
+ points.pop();
+ delete ptr;
+ };
+ virtual bool empty(){ return points.empty(); };
+
+protected:
+ queue<smoothness_vertex_pair*> points;
};
-
-
-
-
-
-
-////class 3DNode{
-//// private:
-//// SPoint3 point;
-//// public:
-//// double min[3];
-//// double max[3];
-//// Node(SPoint3);
-//// Node(double x, double y, double z);
-//// ~Node();
-////};
-////
-
-
-
-
-
#endif
diff --git a/Mesh/surfaceFiller.h b/Mesh/surfaceFiller.h
index e7fba84897..dabb87040b 100644
--- a/Mesh/surfaceFiller.h
+++ b/Mesh/surfaceFiller.h
@@ -7,12 +7,15 @@
#ifndef _SURFACEFILLER_H_
#define _SURFACEFILLER_H_
-
-#include "STensor3.h"
#include <vector>
+#include "STensor3.h"
+
class GFace;
class MVertex;
-void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );
-void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );
+
+void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed,
+ std::vector<SMetric3> &metrics);
+void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed,
+ std::vector<SMetric3> &metrics);
#endif
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 24382028e6..ec356ca73d 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -3,11 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
//
-// Contributor(s):
-// Tristan Carrier
-
-// FIXME: This code should be cleaned up and made to conform with Gmsh's coding
-// style
+// Contributed by Tristan Carrier and Paul-Emile Bernard
#include <numeric>
#include <iterator>
@@ -29,7 +25,7 @@
void export_gregion_mesh(GRegion *gr, string filename)
{
- // FIXME: why not use MElement::writeMSH !?
+ // FIXME: use MElement::writeMSH!
// create set of all tets
map<MVertex*,int> vertices;
@@ -255,7 +251,8 @@ void clique_stop_criteria<T>::export_corresponding_mesh
}
// create MHexahedron, remove included tets from set "tets"
- for (typename graph_data_no_hash::const_iterator it = clique.begin();it!=clique.end();it++){
+ for (typename graph_data_no_hash::const_iterator it = clique.begin();
+ it!=clique.end(); it++){
typename map<T, std::set<MElement*> >::const_iterator itfind = hex_to_tet.find(*it);
if (itfind==hex_to_tet.end()){
cout << "clique_stop_criteria::void export_corresponding_mesh : not found !!!" << endl;
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index 2fc5405046..19f0d357b8 100644
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -3,8 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
//
-// Contributor(s):
-// Tristan Carrier
+// Contributed by Tristan Carrier
#ifndef _YAMAKAWA_H_
#define _YAMAKAWA_H_
--
GitLab