diff --git a/Fltk/FlGui.cpp b/Fltk/FlGui.cpp
index e1e00bdbfb5d4f4f52898e280232aa74f2f99ecd..cbbe33eb1859f0fe5f9fbf5951b6209810c55eb9 100644
--- a/Fltk/FlGui.cpp
+++ b/Fltk/FlGui.cpp
@@ -269,7 +269,8 @@ FlGui::FlGui(int argc, char **argv)
Fl_Tooltip::color(FL_LIGHT2);
#endif
- //Fl::option(Fl::OPTION_HIGH_RES_GL, 1);
+ // FIXME RETINA
+ //Fl::use_high_res_GL(1);
// register image formats not in core fltk library (jpeg/png)
fl_register_images();
diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index ec91132d2433318273392166744c7bb5ff4b9fec..e97194a79edf85e271ee0b953601cee0ef15eb59 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -3103,8 +3103,11 @@ void graphicWindow::detachMenu()
_tile->remove(_onelab);
_browser->resize(0, _browser->y(), _browser->w() + w, _browser->h());
for(unsigned int i = 0; i < gl.size(); i++){
- if(gl[i]->x() == w)
+ if(gl[i]->x() == w){
gl[i]->resize(0, gl[i]->y(), gl[i]->w() + w, gl[i]->h());
+ // FIXME RETINA
+ //gl[i]->redraw();
+ }
}
_tile->redraw();
@@ -3142,8 +3145,11 @@ void graphicWindow::attachMenu()
if(_browser->w() - w < 0) w = _browser->w() / 2;
_browser->resize(w, _browser->y(), _browser->w() - w, _browser->h());
for(unsigned int i = 0; i < gl.size(); i++){
- if(gl[i]->x() == 0)
+ if(gl[i]->x() == 0){
gl[i]->resize(w, gl[i]->y(), gl[i]->w() - w, gl[i]->h());
+ // FIXME RETINA
+ //gl[i]->redraw();
+ }
}
_onelab->box(GMSH_SIMPLE_RIGHT_BOX);
_tile->add(_onelab);
@@ -3331,8 +3337,11 @@ void graphicWindow::setMenuWidth(int w)
double dw = w - _onelab->w();
if(!dw) return;
for(unsigned int i = 0; i < gl.size(); i++){
- if(gl[i]->x() == _onelab->x() + _onelab->w())
+ if(gl[i]->x() == _onelab->x() + _onelab->w()){
gl[i]->resize(gl[i]->x() + dw, gl[i]->y(), gl[i]->w() - dw, gl[i]->h());
+ // FIXME RETINA
+ //gl[i]->redraw();
+ }
}
_browser->resize(_browser->x() + dw, _browser->y(), _browser->w() - dw, _browser->h());
_onelab->resize(_onelab->x(), _onelab->y(), _onelab->w() + dw, _onelab->h());
@@ -3373,8 +3382,11 @@ void graphicWindow::setMessageHeight(int h)
int dh = h - _browser->h();
if(!dh) return;
for(unsigned int i = 0; i < gl.size(); i++){
- if(gl[i]->y() + gl[i]->h() == _browser->y())
+ if(gl[i]->y() + gl[i]->h() == _browser->y()){
gl[i]->resize(gl[i]->x(), gl[i]->y(), gl[i]->w(), gl[i]->h() - dh);
+ // FIXME RETINA
+ //gl[i]->redraw();
+ }
}
_browser->resize(_browser->x(), _browser->y() - dh,
_browser->w(), _browser->h() + dh);
diff --git a/Fltk/openglWindow.cpp b/Fltk/openglWindow.cpp
index 775f20556e3eb38d4a55439e047f4488c1ea472a..aac3a42f8250a70aa062b790e270b1d390d1f184 100644
--- a/Fltk/openglWindow.cpp
+++ b/Fltk/openglWindow.cpp
@@ -159,7 +159,10 @@ void openglWindow::draw()
_ctx->viewport[3] = h();
glViewport(_ctx->viewport[0], _ctx->viewport[1],
_ctx->viewport[2], _ctx->viewport[3]);
- //2*_ctx->viewport[2], 2*_ctx->viewport[3]);
+
+ // FIXME RETINA
+ //glViewport(_ctx->viewport[0], _ctx->viewport[1],
+ // 2*_ctx->viewport[2], 2*_ctx->viewport[3]);
if(lassoMode) {
// draw the zoom or selection lasso on top of the current scene (without
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 53ecf79fdd89e12e917ea8ccb2160800799527d3..3f3c9ef0601ab685492aac9b2d4a648e932047e1 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -33,7 +33,7 @@
#include "Context.h"
#include "OS.h"
#include "GEdgeLoop.h"
-#include "MVertexPositionSet.h"
+#include "MVertexRTree.h"
#include "OpenFile.h"
#include "CreateFile.h"
#include "Options.h"
@@ -1485,16 +1485,8 @@ void GModel::checkMeshCoherence(double tolerance)
for(unsigned int i = 0; i < entities.size(); i++)
vertices.insert(vertices.end(), entities[i]->mesh_vertices.begin(),
entities[i]->mesh_vertices.end());
- MVertexPositionSet pos(vertices);
- for(unsigned int i = 0; i < vertices.size(); i++)
- pos.find(vertices[i]->x(), vertices[i]->y(), vertices[i]->z(), eps);
- int num = 0;
- for(unsigned int i = 0; i < vertices.size(); i++)
- if(!vertices[i]->getIndex()){
- Msg::Info("Duplicate vertex %d at (%.16g,%.16g,%.16g)", vertices[i]->getNum(),
- vertices[i]->x(), vertices[i]->y(), vertices[i]->z());
- num++;
- }
+ MVertexRTree pos(eps);
+ int num = pos.insert(vertices, true);
if(num) Msg::Error("%d duplicate vert%s", num, num > 1 ? "ices" : "ex");
}
@@ -1502,7 +1494,7 @@ void GModel::checkMeshCoherence(double tolerance)
{
Msg::Info("Checking for duplicate elements...");
std::vector<MVertex*> vertices;
- for(unsigned int i = 0; i < entities.size(); i++)
+ for(unsigned int i = 0; i < entities.size(); i++){
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(j);
double vol = e->getVolume();
@@ -1513,18 +1505,11 @@ void GModel::checkMeshCoherence(double tolerance)
SPoint3 p = e->barycenter();
vertices.push_back(new MVertex(p.x(), p.y(), p.z()));
}
- MVertexPositionSet pos(vertices);
+ }
+ MVertexRTree pos(eps);
+ int num = pos.insert(vertices, true);
for(unsigned int i = 0; i < vertices.size(); i++)
- pos.find(vertices[i]->x(), vertices[i]->y(), vertices[i]->z(), eps);
- int num = 0;
- for(unsigned int i = 0; i < vertices.size(); i++){
- if(!vertices[i]->getIndex()){
- Msg::Info("Duplicate element with barycenter (%.16g,%.16g,%.16g)",
- vertices[i]->x(), vertices[i]->y(), vertices[i]->z());
- num++;
- }
delete vertices[i];
- }
if(num) Msg::Error("%d duplicate element%s", num, num > 1 ? "s" : "");
}
@@ -1543,17 +1528,14 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
getEntities(entities);
std::vector<MVertex*> vertices;
- for(unsigned int i = 0; i < entities.size(); i++)
+ 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];
vertices.push_back(new MVertex(v->x(), v->y(), v->z()));
}
- MVertexPositionSet pos(vertices);
- for(unsigned int i = 0; i < vertices.size(); i++)
- pos.find(vertices[i]->x(), vertices[i]->y(), vertices[i]->z(), eps);
- int num = 0;
- for(unsigned int i = 0; i < vertices.size(); i++)
- if(!vertices[i]->getIndex()) num++;
+ }
+ MVertexRTree pos(eps);
+ int num = pos.insert(vertices);
Msg::Info("Found %d duplicate vertices ", num);
@@ -1571,7 +1553,7 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
std::vector<MVertex*> verts;
for(int k = 0; k < e->getNumVertices(); k++){
MVertex *v = e->getVertex(k);
- MVertex *v2 = pos.find(v->x(), v->y(), v->z(), eps);
+ MVertex *v2 = pos.find(v->x(), v->y(), v->z());
if(v2) verts.push_back(v2);
}
if((int)verts.size() == e->getNumVertices()){
@@ -2898,8 +2880,7 @@ static void computeDuplicates(GModel *model,
std::map<GVertex*, GVertex*> &Duplicates2Unique,
const double &eps)
{
- // FIXME: currently we use a greedy algorithm in n^2 (using a kd-tree:
- // cf. MVertexPositionSet)
+ // FIXME: currently we use a greedy algorithm in n^2 (use e.g. MVertexRTree)
// FIXME: add option to remove orphaned entities after duplicate check
std::list<GVertex*> v;
diff --git a/Geo/GModelIO_STL.cpp b/Geo/GModelIO_STL.cpp
index 229b249ffdec4f504f33b58c244c98f5b0140603..476477ff8f27156fb5cdbbb698f069bf571fc077 100644
--- a/Geo/GModelIO_STL.cpp
+++ b/Geo/GModelIO_STL.cpp
@@ -9,7 +9,7 @@
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
-#include "MVertexPositionSet.h"
+#include "MVertexRTree.h"
#include "discreteFace.h"
#include "StringUtils.h"
@@ -143,7 +143,8 @@ int GModel::readSTL(const std::string &name, double tolerance)
for(unsigned int j = 0; j < points[i].size(); j++)
vertices.push_back(new MVertex(points[i][j].x(), points[i][j].y(),
points[i][j].z()));
- MVertexPositionSet pos(vertices);
+ MVertexRTree pos(eps);
+ pos.insert(vertices);
std::set<MFace,Less_Face> unique;
int nbDuplic = 0;
@@ -154,7 +155,7 @@ int GModel::readSTL(const std::string &name, double tolerance)
double x = points[i][j + k].x();
double y = points[i][j + k].y();
double z = points[i][j + k].z();
- v[k] = pos.find(x, y, z, eps);
+ v[k] = pos.find(x, y, z);
}
MFace mf (v[0], v[1], v[2]);
if (unique.find(mf) == unique.end()){
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index b3a0d8bfb6b7f72b82faeea400ca067ccdf4c0bc..b2fc854b903e390be35ae7b05f888d3e20fcc447 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -17,7 +17,7 @@
#include "SBoundingBox3d.h"
#include "OS.h"
#include "MVertex.h"
-#include "MVertexPositionSet.h"
+#include "MVertexRTree.h"
#if defined(HAVE_SOLVER)
#include "linearSystem.h"
@@ -107,11 +107,12 @@ GRbf::GRbf(double sizeBox, int variableEps, int rbfFun,
}
// then create Mvertex position
- std::vector<MVertex*> vertices( allNodes.begin(), allNodes.end() );
- MVertexPositionSet pos(vertices);
+ std::vector<MVertex*> vertices(allNodes.begin(), allNodes.end());
+ MVertexRTree pos(tol);
+ pos.insert(vertices);
for(unsigned int i = 0; i < vertices.size(); i++){
MVertex *v = vertices[i];
- pos.find(v->x(), v->y(), v->z(), tol);
+ pos.find(v->x(), v->y(), v->z());
allCenters(i,0) = v->x()/sBox;
allCenters(i,1) = v->y()/sBox;
allCenters(i,2) = v->z()/sBox;
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 9e9552cd2b5e59f0f7cf231ffe4ef86c59bf9ce0..eef456456261c8186eb9de243ac767ea16953580 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -11,7 +11,7 @@
#include "GModel.h"
#include "GeoInterpolation.h"
#include "Context.h"
-#include "MVertexPositionSet.h"
+#include "MVertexRTree.h"
#if defined(HAVE_MESH)
#include "Field.h"
@@ -2379,14 +2379,15 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
v2V[v] = V;
}
List_Delete(tmp);
- MVertexPositionSet pos(all);
+ MVertexRTree pos(tol);
+ pos.insert(all);
// touch all points
tmp = Tree2List(GModel::current()->getGEOInternals()->Points);
for(int i = 0; i < List_Nbr(tmp); i++) {
Vertex *V;
List_Read(tmp, i, &V);
- pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol);
+ pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z);
}
List_Delete(tmp);
@@ -2397,21 +2398,21 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
Curve *c;
List_Read(tmp, i, &c);
// replace begin/end points
- c->beg = v2V[pos.find(c->beg->Pos.X, c->beg->Pos.Y, c->beg->Pos.Z, tol)];
- c->end = v2V[pos.find(c->end->Pos.X, c->end->Pos.Y, c->end->Pos.Z, tol)];
+ c->beg = v2V[pos.find(c->beg->Pos.X, c->beg->Pos.Y, c->beg->Pos.Z)];
+ c->end = v2V[pos.find(c->end->Pos.X, c->end->Pos.Y, c->end->Pos.Z)];
// replace control points
for(int j = 0; j < List_Nbr(c->Control_Points); j++) {
Vertex *V;
List_Read(c->Control_Points, j, &V);
List_Write(c->Control_Points, j,
- &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]);
+ &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z)]);
}
// replace extrusion sources
if(c->Extrude && c->Extrude->geo.Mode == EXTRUDED_ENTITY){
Vertex *V = FindPoint(std::abs(c->Extrude->geo.Source));
if(V) c->Extrude->geo.Source =
- v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]->Num;
+ v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z)]->Num;
}
}
List_Delete(tmp);
@@ -2426,7 +2427,7 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
Vertex *V;
List_Read(s->TrsfPoints, j, &V);
List_Write(s->TrsfPoints, j,
- &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]);
+ &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z)]);
}
}
List_Delete(tmp);
@@ -2441,7 +2442,7 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
Vertex *V;
List_Read(vol->TrsfPoints, j, &V);
List_Write(vol->TrsfPoints, j,
- &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]);
+ &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z)]);
}
}
List_Delete(tmp);
@@ -2456,7 +2457,7 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
List_Read(p->Entities, j, &num);
Vertex *V = FindPoint(std::abs(num));
if(V) List_Write(p->Entities, j,
- &(v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]->Num));
+ &(v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z)]->Num));
}
}
}
@@ -4287,7 +4288,7 @@ void setSurfaceEmbeddedCurves(Surface *s, List_T *curves)
List_Read(cToAddInSurf->Control_Points, l, &w1);
List_Read(cToAddInSurf->Control_Points, l+1, &w2);
- if (w1 == v1 || w1 == v2 || w2 == v1 || w2 == v2)continue;
+ if (w1 == v1 || w1 == v2 || w2 == v1 || w2 == v2)continue;
SPoint3 q1 = SPoint3(w1->Pos.X, w1->Pos.Y, w1->Pos.Z);
SPoint3 q2 = SPoint3(w2->Pos.X, w2->Pos.Y, w2->Pos.Z);
diff --git a/Geo/MVertexPositionSet.h b/Geo/MVertexPositionSet.h
deleted file mode 100644
index d8f3fc2c426416f6fab87582246ab03b17014ab2..0000000000000000000000000000000000000000
--- a/Geo/MVertexPositionSet.h
+++ /dev/null
@@ -1,91 +0,0 @@
-// Gmsh - Copyright (C) 1997-2014 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>.
-
-#ifndef _MVERTEX_POSITION_SET_
-#define _MVERTEX_POSITION_SET_
-
-#include <vector>
-#include "GmshMessage.h"
-#include "MVertex.h"
-
-#if defined(HAVE_ANN)
-#include "ANN/ANN.h"
-
-// Stores MVertices in a kd-tree so we can query unique vertices (up
-// to a prescribed tolerance). The constructor tags all the vertices
-// with 0; find() tags the returned vertex with -1; if no
-// negatively-tagged vertex exists, find() returns the closest vertex
-// up to the prescribed tolerance.
-class MVertexPositionSet{
- private:
- ANNkd_tree *_kdtree;
- ANNpointArray _zeronodes;
- ANNidxArray _index;
- ANNdistArray _dist;
- int _maxDuplicates;
- std::vector<MVertex*> &_vertices;
- public:
- MVertexPositionSet(std::vector<MVertex*> &vertices, int maxDuplicates=10)
- : _kdtree(0), _maxDuplicates(maxDuplicates), _vertices(vertices)
- {
- int totpoints = vertices.size();
- if(!totpoints) return;
- if(_maxDuplicates > totpoints) _maxDuplicates = totpoints;
- _zeronodes = annAllocPts(totpoints, 3);
- for(int i = 0; i < totpoints; i++){
- vertices[i]->setIndex(0);
- _zeronodes[i][0] = vertices[i]->x();
- _zeronodes[i][1] = vertices[i]->y();
- _zeronodes[i][2] = vertices[i]->z();
- }
- _kdtree = new ANNkd_tree(_zeronodes, totpoints, 3);
- _index = new ANNidx[_maxDuplicates];
- _dist = new ANNdist[_maxDuplicates];
- }
- ~MVertexPositionSet()
- {
- if(!_kdtree) return;
- delete _kdtree;
- annDeallocPts(_zeronodes);
- delete [] _index;
- delete [] _dist;
- }
- MVertex *find(double x, double y, double z, double tolerance)
- {
- if(!_kdtree) return 0;
- double xyz[3] = {x, y, z};
- _kdtree->annkSearch(xyz, _maxDuplicates, _index, _dist);
- for(int i = 0; i < _maxDuplicates; i++){
- if(_index[i] >= 0 && sqrt(_dist[i]) < tolerance &&
- _vertices[_index[i]]->getIndex() < 0)
- return _vertices[_index[i]];
- }
- if(_index[0] >= 0 && sqrt(_dist[0]) < tolerance){
- // printf("tol %g dist %g x %g %g %g\n",tolerance,_dist[0],x,y,z);
- _vertices[_index[0]]->setIndex(-1);
- return _vertices[_index[0]];
- }
- Msg::Error("Could not find vertex (%g,%g,%g) (tol %g)",
- x, y, z, tolerance);
- return 0;
- }
-};
-
-#else
-
-class MVertexPositionSet{
- public:
- MVertexPositionSet(std::vector<MVertex*> &vertices, int maxDuplicates=10){}
- ~MVertexPositionSet(){}
- MVertex *find(double x, double y, double z, double tolerance)
- {
- Msg::Error("Gmsh must be compiled with ANN to use MVertexPositionSet");
- return 0;
- }
-};
-
-#endif
-
-#endif
diff --git a/Geo/MVertexRTree.h b/Geo/MVertexRTree.h
new file mode 100644
index 0000000000000000000000000000000000000000..697c535b9f9cb911fceaafbf10005c15a93f7504
--- /dev/null
+++ b/Geo/MVertexRTree.h
@@ -0,0 +1,74 @@
+// Gmsh - Copyright (C) 1997-2014 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>.
+
+#ifndef _MVERTEX_RTREE_
+#define _MVERTEX_RTREE_
+
+#include <vector>
+#include "GmshMessage.h"
+#include "MVertex.h"
+#include "rtree.h"
+
+// Stores MVertex pointers in an R-Tree so we can query unique vertices by their
+// coordinates, up to a prescribed tolerance.
+class MVertexRTree{
+ private:
+ RTree<MVertex*, double, 3, double> *_rtree;
+ double _tol;
+ static bool rtree_callback(MVertex *v, void *ctx)
+ {
+ MVertex **out = static_cast<MVertex**>(ctx);
+ *out = v;
+ return false;
+ }
+ public:
+ MVertexRTree(double tolerance = 1.e-8)
+ {
+ _rtree = new RTree<MVertex*, double, 3, double>();
+ _tol = tolerance;
+ }
+ ~MVertexRTree()
+ {
+ _rtree->RemoveAll();
+ delete _rtree;
+ }
+ int insert(MVertex *v, bool warnIfExists=false)
+ {
+ MVertex *out;
+ double _min[3] = {v->x() - _tol, v->y() - _tol, v->z() - _tol};
+ double _max[3] = {v->x() + _tol, v->y() + _tol, v->z() + _tol};
+ if(!_rtree->Search(_min, _max, rtree_callback, &out)){
+ _rtree->Insert(_min, _max, v);
+ return 0;
+ }
+ else if(warnIfExists){
+ Msg::Warning("Vertex %d (%.16g, %.16g, %.16g) already exists in the "
+ "mesh with tolerance %g", _tol);
+ }
+ return 1;
+ }
+ int insert(std::vector<MVertex*> &v, bool warnIfExists=false)
+ {
+ int num = 0;
+ for(unsigned int i = 0; i < v.size(); i++)
+ num += insert(v[i], warnIfExists);
+ return num;
+ }
+ MVertex *find(double x, double y, double z)
+ {
+ MVertex *out;
+ double _min[3] = {x - _tol, y - _tol, z - _tol};
+ double _max[3] = {x + _tol, y + _tol, z + _tol};
+ if(_rtree->Search(_min, _max, rtree_callback, &out))
+ return out;
+ return 0;
+ }
+ unsigned int size()
+ {
+ return _rtree->Count();
+ }
+};
+
+#endif
diff --git a/Mesh/QuadTriExtruded2D.cpp b/Mesh/QuadTriExtruded2D.cpp
index 0af475873ddf6860bea7d2429210fa2a19aa1cd3..75a704842f5c7b8ebf86251e68c5f1f5176dd74b 100644
--- a/Mesh/QuadTriExtruded2D.cpp
+++ b/Mesh/QuadTriExtruded2D.cpp
@@ -1,36 +1,37 @@
-/**************************************************************************************************
+/*******************************************************************************
QuadTriExtruded2D.cpp
The code in this file was written by Dr. Trevor S. Strickler.
email: <trevor.strickler@gmail.com>
-This file is part of the QuadTri contribution to Gmsh. QuadTri allows the conformal interface
-of quadrangle faces to triangle faces using pyramids and other mesh elements.
+This file is part of the QuadTri contribution to Gmsh. QuadTri allows the
+conformal interface of quadrangle faces to triangle faces using pyramids and
+other mesh elements.
-See READMEQUADTRI.txt for more information. The license information is in LICENSE.txt.
+See READMEQUADTRI.txt for more information. The license information is in
+LICENSE.txt.
-Trevor S. Strickler hereby transfers copyright of QuadTri files to
-Christophe Geuzaine and J.-F. Remacle with the understanding that
-his contribution shall be cited appropriately.
+Trevor S. Strickler hereby transfers copyright of QuadTri files to Christophe
+Geuzaine and J.-F. Remacle with the understanding that his contribution shall be
+cited appropriately.
-All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
-Gmsh is available at: www.geuz.org/gmsh
+All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine,
+J.-F. Remacle Gmsh is available at: www.geuz.org/gmsh
-For Gmsh license information, see the LICENSE.txt file for license information. Please report all
-Gmsh bugs and problems to the public mailing list <gmsh@geuz.org>.
+For Gmsh license information, see the LICENSE.txt file for license
+information. Please report all Gmsh bugs and problems to the public mailing list
+<gmsh@geuz.org>.
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License, Version 2,
-as published by the Free Software Foundation, or (at your option)
-any later version, with or without the exception given in the
-LICENSE.txt file supplied with this code and with Gmsh.
+This program is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License, Version 2, as published by the Free
+Software Foundation, or (at your option) any later version, with or without the
+exception given in the LICENSE.txt file supplied with this code and with Gmsh.
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
+This program is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE. See the GNU General Public License for more details.
-****************************************************************************************************/
+********************************************************************************/
#include "QuadTriExtruded2D.h"
@@ -42,44 +43,41 @@ static void addTriangle(MVertex* v1, MVertex* v2, MVertex* v3,
to->triangles.push_back(newTri);
}
-
-
-// The function that tests whether a 2D surface is a lateral of a valid QuadToTri
-// region and whether there are conflicts. If surface is not part of valid QuadToTri region
-// or if there are QuadToTri conflicts, return 0. Note that RemoveDuplicateSurfaces()
-// makes this DIFFICULT. Also, the tri_quad_flag determins whether the surface
-// should be meshed with triangles or quadrangles:
+// The function that tests whether a 2D surface is a lateral of a valid
+// QuadToTri region and whether there are conflicts. If surface is not part of
+// valid QuadToTri region or if there are QuadToTri conflicts, return 0. Note
+// that RemoveDuplicateSurfaces() makes this DIFFICULT. Also, the tri_quad_flag
+// determins whether the surface should be meshed with triangles or quadrangles:
// tri_quad_values: 0 = no override, 1 = mesh as quads, 2 = mesh as triangles.
// Added 2010-12-09.
-int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToTriLateral )
+int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToTriLateral)
{
(*tri_quad_flag) = 0;
(*detectQuadToTriLateral) = false;
GModel *model = face->model();
-
ExtrudeParams *ep = face->meshAttributes.extrude;
- if( !ep || !ep->mesh.ExtrudeMesh || !(ep->geo.Mode == EXTRUDED_ENTITY) ){
+ if(!ep || !ep->mesh.ExtrudeMesh || !(ep->geo.Mode == EXTRUDED_ENTITY)){
Msg::Error("In IsValidQuadToTriLateral(), face %d is not a structured extrusion.",
face->tag() );
return 0;
}
- GEdge *face_source = model->getEdgeByTag( std::abs( ep->geo.Source ) );
- if( !face_source ){
+ GEdge *face_source = model->getEdgeByTag(std::abs(ep->geo.Source));
+ if(!face_source){
Msg::Error("In IsValidQuadToTriLateral(), face %d has no source edge.",
- face->tag() );
+ face->tag());
}
- // It seems the member pointers to neighboring regions for extruded lateral faces are not set.
- // For now, have to loop through
- // ALL the regions to see if the presently considered face belongs to the region and if
- // any neighboring region is QUADTRI.
- // The following loop will find all the regions that the face bounds, and determine
- // whether the face is a lateral of the region (including whether the region is even extruded).
- // After that information is determined, function can test for QuadToTri neighbor conflicts.
+ // It seems the member pointers to neighboring regions for extruded lateral
+ // faces are not set. For now, have to loop through ALL the regions to see if
+ // the presently considered face belongs to the region and if any neighboring
+ // region is QUADTRI. The following loop will find all the regions that the
+ // face bounds, and determine whether the face is a lateral of the region
+ // (including whether the region is even extruded). After that information is
+ // determined, function can test for QuadToTri neighbor conflicts.
std::vector<GRegion *> lateral_regions;
std::vector<GRegion *> adjacent_regions;
@@ -87,26 +85,26 @@ int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToT
int numLateralRegions = 0;
numRegions = GetNeighborRegionsOfFace(face, adjacent_regions);
- for( int i_reg = 0; i_reg < numRegions; i_reg++ ){
+ for(int i_reg = 0; i_reg < numRegions; i_reg++){
GRegion *region = adjacent_regions[i_reg];
// is region in the current model's region's or is it deleted?
- if( !FindVolume( ( region->tag() ) ) )
+ if(!FindVolume((region->tag())))
continue;
// is the region mesh extruded?
- if( !region->meshAttributes.extrude ||
- ( region->meshAttributes.extrude &&
- !region->meshAttributes.extrude->mesh.ExtrudeMesh ) )
+ if(!region->meshAttributes.extrude ||
+ (region->meshAttributes.extrude &&
+ !region->meshAttributes.extrude->mesh.ExtrudeMesh))
continue;
- if( region->meshAttributes.extrude->geo.Mode != EXTRUDED_ENTITY )
+ if(region->meshAttributes.extrude->geo.Mode != EXTRUDED_ENTITY)
continue;
// Test whether the face is a lateral
- if( IsSurfaceALateralForRegion(region, face) ){
+ if(IsSurfaceALateralForRegion(region, face)){
lateral_regions.push_back(region);
numLateralRegions++;
- if( region->meshAttributes.extrude->mesh.QuadToTri )
+ if(region->meshAttributes.extrude->mesh.QuadToTri)
(*detectQuadToTriLateral) = true;
}
@@ -114,7 +112,7 @@ int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToT
// MAIN test of whether this is even a quadToTri extrusion lateral
// the only return 0 path that is NOT an error
- if( !(*detectQuadToTriLateral) )
+ if(!(*detectQuadToTriLateral))
return 0;
// now will start conflict checks
@@ -127,18 +125,17 @@ int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToT
bool detect_conflict = false;
-
// Set the tri_quad_flag that lets ExtrudeMesh override ep->Recombine;
// tri_quad_values: 0 = no override, 1 = mesh as quads, 2 = mesh as triangles.
// if this face is a free surface:
if( adjacent_regions.size() == 1 ){
- if( lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1_RECOMB ||
- lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1_RECOMB ){
+ if(lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1_RECOMB ||
+ lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1_RECOMB ){
(*tri_quad_flag) = 1;
}
- if( lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1 ||
- lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1 ){
+ if(lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1 ||
+ lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1 ){
(*tri_quad_flag) = 2;
}
else
@@ -165,25 +162,25 @@ int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToT
// neighboring extrusion.
else if( adj_ep && adj_ep->mesh.ExtrudeMesh &&
model->getFaceByTag( std::abs( adj_ep->geo.Source ) ) == face ){
- if( lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1_RECOMB ||
- lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1_RECOMB )
+ if(lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1_RECOMB ||
+ lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1_RECOMB )
(*tri_quad_flag) = 1;
- else if( lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1 ||
- lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1 )
+ else if(lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_NOVERTS_1 ||
+ lateral_regions[0]->meshAttributes.extrude->mesh.QuadToTri == QUADTRI_ADDVERTS_1 )
(*tri_quad_flag) = 2;
else
(*tri_quad_flag) = 0;
}
// if both neighbors are structured but none of the previous apply:
- else if( adj_ep && adj_ep->mesh.ExtrudeMesh ){
- if( (adj_ep && !adj_ep->mesh.QuadToTri && adj_ep->mesh.Recombine) ||
- (ep && !ep->mesh.QuadToTri && ep->mesh.Recombine) )
+ else if(adj_ep && adj_ep->mesh.ExtrudeMesh){
+ if((adj_ep && !adj_ep->mesh.QuadToTri && adj_ep->mesh.Recombine) ||
+ (ep && !ep->mesh.QuadToTri && ep->mesh.Recombine))
(*tri_quad_flag) = 1;
- else if( (adj_ep && !adj_ep->mesh.QuadToTri && !adj_ep->mesh.Recombine) ||
- (ep && !ep->mesh.QuadToTri && !ep->mesh.Recombine) )
+ else if((adj_ep && !adj_ep->mesh.QuadToTri && !adj_ep->mesh.Recombine) ||
+ (ep && !ep->mesh.QuadToTri && !ep->mesh.Recombine))
(*tri_quad_flag) = 2;
// if both are quadToTri ALWAYS try to recombine
- else if( ep->mesh.QuadToTri && adj_ep && adj_ep->mesh.QuadToTri )
+ else if(ep->mesh.QuadToTri && adj_ep && adj_ep->mesh.QuadToTri)
(*tri_quad_flag) = 1;
else
(*tri_quad_flag) = 0;
@@ -199,24 +196,21 @@ int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToT
(*tri_quad_flag) = 0;
}
- if( detect_conflict )
+ if(detect_conflict)
return 0;
else
return 1;
}
-
-
-// The function that tests whether a surface is a QuadToTri top surface and whether
-// there are conflicts. If surface is not a top for a valid QuadToTri region or if
-// there are QuadToTri conflicts, return 0.
-// if the surface turns out to be the source of a toroidal loop extrusion (which will then
-// NOT have geo.Mode == COPIED_ENTITY), return 2 (this will require special meshing considerations).
-// Note that RemoveDuplicateSurfaces() makes this DIFFICULT.
-// Also, the type of QuadToTri interface is placed into the
-// pointer argument quadToTri. .
-// Added 2010-12-09.
+// The function that tests whether a surface is a QuadToTri top surface and
+// whether there are conflicts. If surface is not a top for a valid QuadToTri
+// region or if there are QuadToTri conflicts, return 0. if the surface turns
+// out to be the source of a toroidal loop extrusion (which will then NOT have
+// geo.Mode == COPIED_ENTITY), return 2 (this will require special meshing
+// considerations). Note that RemoveDuplicateSurfaces() makes this DIFFICULT.
+// Also, the type of QuadToTri interface is placed into the pointer argument
+// quadToTri. . Added 2010-12-09.
int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop)
{
(*quadToTri) = NO_QUADTRI;
@@ -226,18 +220,16 @@ int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop)
GModel *model = face->model();
- // First thing is first: determine if this is a toroidal quadtri extrusion. if so, can skip the rest
-
-
- // It seems the member pointers to neighboring regions for extruded top faces are not set.
- // For now, have to loop through
- // ALL the regions to see if the presently considered face belongs to the region.
- // The following loop will find all the regions that the face bounds, and determine
- // whether the face is a top face of the region (including whether the region is even extruded).
- // After that information is determined, function can test for QuadToTri neighbor conflicts.
-
-
+ // First thing is first: determine if this is a toroidal quadtri extrusion.
+ // if so, can skip the rest
+ // It seems the member pointers to neighboring regions for extruded top faces
+ // are not set. For now, have to loop through ALL the regions to see if the
+ // presently considered face belongs to the region. The following loop will
+ // find all the regions that the face bounds, and determine whether the face
+ // is a top face of the region (including whether the region is even
+ // extruded). After that information is determined, function can test for
+ // QuadToTri neighbor conflicts.
// first determine if this is toroidal quadtotri
is_toroidal_quadtri = IsInToroidalQuadToTri(face);
@@ -252,7 +244,7 @@ int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop)
int numTopRegions = 0;
std::set<GRegion *, GEntityLessThan>::iterator itreg;
- for( itreg = model->firstRegion(); itreg != model->lastRegion(); itreg++ )
+ for(itreg = model->firstRegion(); itreg != model->lastRegion(); itreg++)
all_regions.push_back( (*itreg) );
for(unsigned int i_reg = 0; i_reg < all_regions.size(); i_reg++ ){
@@ -278,7 +270,8 @@ int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop)
continue;
// is region a structured extruded?
- if( !(region->meshAttributes.extrude && region->meshAttributes.extrude->mesh.ExtrudeMesh &&
+ if( !(region->meshAttributes.extrude &&
+ region->meshAttributes.extrude->mesh.ExtrudeMesh &&
region->meshAttributes.extrude->geo.Mode == EXTRUDED_ENTITY) )
continue;
@@ -312,8 +305,9 @@ int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop)
}
if( ep->mesh.QuadToTri == 0){
- Msg::Error("In IsValidQuadToTriTop(), surface %d was determined to be the top surface "
- "for a QuadToTri extrusion, but does not have QuadToTri parameters set within itself.",
+ Msg::Error("In IsValidQuadToTriTop(), surface %d was determined to be the "
+ "top surface for a QuadToTri extrusion, but does not have "
+ "QuadToTri parameters set within itself.",
face->tag() );
return 0;
}
@@ -331,23 +325,22 @@ int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop)
return 0;
}
-
- if( top_regions.size() ){
+ if(top_regions.size()){
(*quadToTri) = top_regions[0]->meshAttributes.extrude->mesh.QuadToTri;
}
// Make sure that face is the top for only one region. if not, then there will likely
// be conflicts (two regions extruded into each other).
- if( top_regions.size() > 1 ){
+ if(top_regions.size() > 1){
Msg::Error("In IsValidQuadToTriTop(), QuadToTri top surface %d identified as top "
- "surface for more than one region. Likely conflict.", face->tag() );
+ "surface for more than one region. Likely conflict.", face->tag() );
return 0;
}
} // end of else that executes if NOT toroidal extrusion
-
- // this is technically redundant...but if changes are made, it's good to keep this here at the end for safety
+ // this is technically redundant...but if changes are made, it's good to
+ // keep this here at the end for safety
if( !(*detectQuadToTriTop) )
return 0;
@@ -360,69 +353,58 @@ int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop)
}
-
// this function specifically meshes a quadToTri top in an unstructured way
// return 1 if success, return 0 if failed.
// Added 2010-12-20
-static int MeshQuadToTriTopUnstructured(GFace *from, GFace *to,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+static int MeshQuadToTriTopUnstructured(GFace *from, GFace *to, MVertexRTree &pos)
{
-
// if the source is all triangles, then just return 1.
- if( from->triangles.size() && !from->quadrangles.size() )
+ if(from->triangles.size() && !from->quadrangles.size())
return 1;
- if( !to->meshAttributes.extrude || !to->meshAttributes.extrude->mesh.QuadToTri )
+ if(!to->meshAttributes.extrude || !to->meshAttributes.extrude->mesh.QuadToTri)
return 0;
-
// in weird case of NO quads and NO tri
- if( !from->triangles.size() && !from->quadrangles.size() )
+ if(!from->triangles.size() && !from->quadrangles.size())
return 0;
-
// make set of source edge vertices
- std::set<MVertex*, MVertexLessThanLexicographic> pos_src_edge;
+ MVertexRTree pos_src_edge(1.e-12 * CTX::instance()->lc);
QuadToTriInsertFaceEdgeVertices(from, pos_src_edge);
// Loop through all the quads and make the triangles with diagonals running
// in a selected direction.
to->triangles.reserve(to->triangles.size()+from->quadrangles.size()*2);
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
for(unsigned int i = 0; i < from->quadrangles.size(); i++){
std::vector<MVertex*> verts;
for(int j = 0; j < from->quadrangles[i]->getNumVertices(); j++){
MVertex *v = from->quadrangles[i]->getVertex(j);
- MVertex tmp(v->x(), v->y(), v->z(), 0, -1);
+ double x = v->x(), y = v->y(), z = v->z();
ExtrudeParams *ep = to->meshAttributes.extrude;
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
- tmp.x(), tmp.y(), tmp.z());
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()) {
+ x, y, z);
+ MVertex *tmp = pos.find(x, y, z);
+ if(!tmp) {
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
- to->triangles.reserve(to->triangles.size()+1);
+ x, y, z, to->tag());
+ to->triangles.reserve(to->triangles.size() + 1);
return 0;
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
-
- if( verts.size() != 4 ){
+ if(verts.size() != 4){
Msg::Error("During mesh of QuadToTri surface %d, %d vertices found "
"in quad of source surface %d.", to->tag(), verts.size(),
from->tag() );
return 0;
}
-
- // draw other diagonals to minimize difference in average edge length with diagonal length, in quadrature
+ // draw other diagonals to minimize difference in average edge length with
+ // diagonal length, in quadrature
double mag_sq_ave = 0.0;
for( int p = 0; p < 4; p++ ){
@@ -447,79 +429,73 @@ static int MeshQuadToTriTopUnstructured(GFace *from, GFace *to,
return 1;
}
-
-// This function meshes the top surface of a QuadToTri extrusion. It returns 0 if it is given a
-// non-quadToTri extrusion or if it fails.
+// This function meshes the top surface of a QuadToTri extrusion. It returns 0
+// if it is given a non-quadToTri extrusion or if it fails.
// Args:
-// 'GFace *to' is the top surface to mesh, 'from' is the source surface, 'pos' is a std::set
-// of vertex positions for the top surface.
-int MeshQuadToTriTopSurface( GFace *from, GFace *to, std::set<MVertex*,
- MVertexLessThanLexicographic> &pos )
+// 'GFace *to' is the top surface to mesh,
+// 'from' is the source surface
+// 'pos' is a tree of vertex positions for the top surface.
+int MeshQuadToTriTopSurface(GFace *from, GFace *to, MVertexRTree &pos)
{
- if( !to->meshAttributes.extrude || !to->meshAttributes.extrude->mesh.QuadToTri )
+ if(!to->meshAttributes.extrude || !to->meshAttributes.extrude->mesh.QuadToTri)
return 0;
// if the source is all triangles, then just let this function is not needed. Return 1.
- if( from->triangles.size() && !from->quadrangles.size() )
+ if(from->triangles.size() && !from->quadrangles.size())
return 1;
// in weird case of NO quads and NO tri
- if( !from->triangles.size() && !from->quadrangles.size() )
+ if(!from->triangles.size() && !from->quadrangles.size())
return 0;
-
ExtrudeParams *ep = to->meshAttributes.extrude;
- if( !ep || !ep->mesh.ExtrudeMesh || !(ep->geo.Mode == COPIED_ENTITY) ){
+ if(!ep || !ep->mesh.ExtrudeMesh || !(ep->geo.Mode == COPIED_ENTITY)){
Msg::Error("In MeshQuadToTriTopSurface(), incomplete or no "
- "extrude information for top face %d.", to->tag() );
+ "extrude information for top face %d.", to->tag() );
return 0;
}
// is this a quadtri extrusion with added vertices?
bool is_addverts = false;
- if( ep && (ep->mesh.QuadToTri == QUADTRI_ADDVERTS_1 || ep->mesh.QuadToTri == QUADTRI_ADDVERTS_1_RECOMB) )
+ if(ep && (ep->mesh.QuadToTri == QUADTRI_ADDVERTS_1 ||
+ ep->mesh.QuadToTri == QUADTRI_ADDVERTS_1_RECOMB))
is_addverts = true;
- // execute this section if
- // IF this is a 'no new vertices' quadToTri, mesh the surfaces according to this modified
- // least point value method: if a 3 boundary point quad, draw diagonals from middle corner toward
- // interior. If a a 2- or 1- point boundary quad, draw toward lowest pointer number NOT on boundary.
- // All interior quad, draw diagonal to vertex with lowest pointer number.
+ // execute this section if IF this is a 'no new vertices' quadToTri, mesh the
+ // surfaces according to this modified least point value method: if a 3
+ // boundary point quad, draw diagonals from middle corner toward interior. If
+ // a a 2- or 1- point boundary quad, draw toward lowest pointer number NOT on
+ // boundary. All interior quad, draw diagonal to vertex with lowest pointer
+ // number.
- if( !is_addverts ){
- std::set<MVertex*, MVertexLessThanLexicographic> pos_src_edge;
+ if(!is_addverts){
+ MVertexRTree pos_src_edge(1.e-12 * CTX::instance()->lc);
QuadToTriInsertFaceEdgeVertices(from, pos_src_edge);
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
// loop through each element source quadrangle and extrude
for(unsigned int i = 0; i < from->quadrangles.size(); i++){
std::vector<MVertex*> verts;
for(int j = 0; j < from->quadrangles[i]->getNumVertices(); j++){
MVertex *v = from->quadrangles[i]->getVertex(j);
- MVertex tmp(v->x(), v->y(), v->z(), 0, -1);
+ double x = v->x(), y = v->y(), z = v->z();
ExtrudeParams *ep = to->meshAttributes.extrude;
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
- tmp.x(), tmp.y(), tmp.z());
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()) {
+ x, y, z);
+ MVertex *tmp = pos.find(x, y, z);
+ if(!tmp) {
Msg::Error("In MeshQuadToTriTopSurface(), Could not find "
"extruded vertex (%.16g, %.16g, %.16g) in surface %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
- to->triangles.reserve(to->triangles.size()+1);
+ x, y, z, to->tag());
+ to->triangles.reserve(to->triangles.size() + 1);
return 0;
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
-
- if( verts.size() != 4 ){
+ if(verts.size() != 4){
Msg::Error("During mesh of QuadToTri surface %d, %d vertices found "
"in quad of source surface %d.", to->tag(), verts.size(),
- from->tag() );
+ from->tag());
return 0;
}
@@ -530,8 +506,10 @@ int MeshQuadToTriTopSurface( GFace *from, GFace *to, std::set<MVertex*,
int edge_verts_count = 0;
//int skip_index = 0;
int bnd_indices[4];
- for( int p = 0; p < element->getNumVertices(); p++ ){
- if( pos_src_edge.find( element->getVertex(p) ) != pos_src_edge.end() ){
+ for(int p = 0; p < element->getNumVertices(); p++){
+ if(pos_src_edge.find(element->getVertex(p)->x(),
+ element->getVertex(p)->y(),
+ element->getVertex(p)->z())){
edge_verts_count++;
bnd_indices[p] = 1;
}
@@ -543,21 +521,22 @@ int MeshQuadToTriTopSurface( GFace *from, GFace *to, std::set<MVertex*,
// Apply modified lowest vertex pointer diagonalization
int low_index = -1;
- if( edge_verts_count == 3 || edge_verts_count == 2 || edge_verts_count == 1 ){
- for( int p = 0; p < 4; p++ ){
- if( !bnd_indices[p] && verts[p] != element->getVertex(p) ){
- if( low_index < 0 )
+ if(edge_verts_count == 3 || edge_verts_count == 2 || edge_verts_count == 1){
+ for(int p = 0; p < 4; p++ ){
+ if(!bnd_indices[p] && verts[p] != element->getVertex(p)){
+ if(low_index < 0 )
low_index = p;
- else if( verts[p] < verts[low_index] )
+ else if(verts[p] < verts[low_index])
low_index = p;
}
}
- if( low_index < 0 ) // what if they are all degenerate? Avoid the out-of-bounds error.
+ if(low_index < 0) // what if they are all degenerate? Avoid the
+ // out-of-bounds error.
low_index = 0;
}
// lowest possible vertex pointer, regardless of if on edge or not
- else if( edge_verts_count == 4 || edge_verts_count == 0 )
+ else if(edge_verts_count == 4 || edge_verts_count == 0)
low_index = getIndexForLowestVertexPointer(verts);
addTriangle( verts[low_index],verts[(low_index+1)%verts.size()],
@@ -569,49 +548,44 @@ int MeshQuadToTriTopSurface( GFace *from, GFace *to, std::set<MVertex*,
}
- // AFTER THIS POINT IN FUNCTION, CODE IS ALL FOR 'ADD INTERNAL VERTEX' EXTRUSIONS (Less restrictive).
+ // AFTER THIS POINT IN FUNCTION, CODE IS ALL FOR 'ADD INTERNAL VERTEX'
+ // EXTRUSIONS (Less restrictive).
// if source face is unstructured, can try to make the top mesh a little neater
GFace *root_source = findRootSourceFaceForFace( from );
ExtrudeParams *ep_src = root_source->meshAttributes.extrude;
bool struct_root = false;
- if( root_source &&
- ( (ep_src && ep_src->mesh.ExtrudeMesh && ep_src->geo.Mode == EXTRUDED_ENTITY) ||
- root_source->meshAttributes.method == MESH_TRANSFINITE ) )
+ if(root_source &&
+ ((ep_src && ep_src->mesh.ExtrudeMesh && ep_src->geo.Mode == EXTRUDED_ENTITY) ||
+ root_source->meshAttributes.method == MESH_TRANSFINITE))
struct_root = true;
- if( !struct_root && MeshQuadToTriTopUnstructured(from, to, pos) )
+ if(!struct_root && MeshQuadToTriTopUnstructured(from, to, pos)){
return 1;
-
- // And top surface for the 'added internal vertex' method can be meshed quite easily
+ }
else{
- std::set<MVertex *, MVertexLessThanLexicographic >::iterator itp;
- // loop through each element source quadrangle and extrude
+ // And top surface for the 'added internal vertex' method can be meshed
+ // quite easily. Loop through each element source quadrangle and extrude
for(unsigned int i = 0; i < from->quadrangles.size(); i++){
std::vector<MVertex*> verts;
for(int j = 0; j < from->quadrangles[i]->getNumVertices(); j++){
MVertex *v = from->quadrangles[i]->getVertex(j);
- MVertex tmp(v->x(), v->y(), v->z(), 0, -1);
+ double x = v->x(), y = v->y(), z = v->z();
ExtrudeParams *ep = to->meshAttributes.extrude;
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
- tmp.x(), tmp.y(), tmp.z());
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()) {
+ x, y, z);
+ MVertex *tmp = pos.find(x, y, z);
+ if(!tmp) {
Msg::Error("In MeshQuadToTriTopSurface(), Could not find "
"extruded vertex (%.16g, %.16g, %.16g) in surface %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
- to->triangles.reserve(to->triangles.size()+1);
+ x, y, z, to->tag());
+ to->triangles.reserve(to->triangles.size() + 1);
return 0;
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
-
- if( verts.size() != 4 ){
+ if(verts.size() != 4){
Msg::Error("During mesh of QuadToTri surface %d, %d vertices found "
"in quad of source surface %d.", to->tag(), verts.size(),
from->tag() );
@@ -626,5 +600,4 @@ int MeshQuadToTriTopSurface( GFace *from, GFace *to, std::set<MVertex*,
}
return 0;
-
}
diff --git a/Mesh/QuadTriExtruded2D.h b/Mesh/QuadTriExtruded2D.h
index 3b1a05d8dc3ef02de035892f3365e576d0673dea..9a6b26fe3eb68054fb81a884243bb773c5976383 100644
--- a/Mesh/QuadTriExtruded2D.h
+++ b/Mesh/QuadTriExtruded2D.h
@@ -1,36 +1,37 @@
-/**************************************************************************************************
+/*******************************************************************************
QuadTriExtruded2D.h
The code in this file was written by Dr. Trevor S. Strickler.
-email: <trevor.strickler@gmail.com>
+email: <trevor.strickler@gmail.com>
-This file is part of the QuadTri contribution to Gmsh. QuadTri allows the conformal interface
-of quadrangle faces to triangle faces using pyramids and other mesh elements.
+This file is part of the QuadTri contribution to Gmsh. QuadTri allows the
+conformal interface of quadrangle faces to triangle faces using pyramids and
+other mesh elements.
-See READMEQUADTRI.txt for more information. The license information is in LICENSE.txt.
+See READMEQUADTRI.txt for more information. The license information is in
+LICENSE.txt.
-Trevor S. Strickler hereby transfers copyright of QuadTri files to
-Christophe Geuzaine and J.-F. Remacle with the understanding that
-his contribution shall be cited appropriately.
+Trevor S. Strickler hereby transfers copyright of QuadTri files to Christophe
+Geuzaine and J.-F. Remacle with the understanding that his contribution shall be
+cited appropriately.
-All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
-Gmsh is available at: www.geuz.org/gmsh
+All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine,
+J.-F. Remacle Gmsh is available at: www.geuz.org/gmsh
-For Gmsh license information, see the LICENSE.txt file for license information. Please report all
-Gmsh bugs and problems to the public mailing list <gmsh@geuz.org>.
+For Gmsh license information, see the LICENSE.txt file for license
+information. Please report all Gmsh bugs and problems to the public mailing list
+<gmsh@geuz.org>.
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License, Version 2,
-as published by the Free Software Foundation, or (at your option)
-any later version, with or without the exception given in the
-LICENSE.txt file supplied with this code and with Gmsh.
+This program is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License, Version 2, as published by the Free
+Software Foundation, or (at your option) any later version, with or without the
+exception given in the LICENSE.txt file supplied with this code and with Gmsh.
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
+This program is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE. See the GNU General Public License for more details.
+********************************************************************************/
-****************************************************************************************************/
#if !defined(_QTEXTR2D_H_)
#define _QTEXTR2D_H_
@@ -53,33 +54,27 @@ GNU General Public License for more details.
#include <math.h>
#include "QuadTriUtils.h"
-
-
-// The function that tests whether a 2D surface is a lateral of a valid QuadToTri
-// region and whether there are conflicts. If surface is not part of valid QuadToTri region
-// or if there are QuadToTri conflicts, return 0. Note that RemoveDuplicateSurfaces()
-// makes this DIFFICULT. Also, the tri_quad_flag determins whether the surface
-// should be meshed with triangles or quadrangles:
-// tri_quad_values: 0 = no override, 1 = mesh as quads, 2 = mesh as triangles.
+// The function that tests whether a 2D surface is a lateral of a valid
+// QuadToTri region and whether there are conflicts. If surface is not part of
+// valid QuadToTri region or if there are QuadToTri conflicts, return 0. Note
+// that RemoveDuplicateSurfaces() makes this DIFFICULT. Also, the tri_quad_flag
+// determins whether the surface should be meshed with triangles or quadrangles:
+// tri_quad_values: 0 = no override, 1 = mesh as quads, 2 = mesh as triangles.
// Added 2010-12-09.
-int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToTriLateral );
+int IsValidQuadToTriLateral(GFace *face, int *tri_quad_flag, bool *detectQuadToTriLateral);
-
-// The function that tests whether a surface is a QuadToTri top surface and whether
-// there are conflicts. If surface is not a top for a valid QuadToTri region or if
-// there are QuadToTri conflicts, return 0. Note that RemoveDuplicateSurfaces()
-// makes this DIFFICULT. Also, the type of QuadToTri interface is placed into the
-// pointer argument quadToTri. .
-// Added 2010-12-09.
+// The function that tests whether a surface is a QuadToTri top surface and
+// whether there are conflicts. If surface is not a top for a valid QuadToTri
+// region or if there are QuadToTri conflicts, return 0. Note that
+// RemoveDuplicateSurfaces() makes this DIFFICULT. Also, the type of QuadToTri
+// interface is placed into the pointer argument quadToTri. . Added 2010-12-09.
int IsValidQuadToTriTop(GFace *face, int *quadToTri, bool *detectQuadToTriTop);
-// This function meshes the top surface of a QuadToTri extrusion. It returns 0 if it is given a
-// non-quadToTri extrusion or if it fails.
-// Args:
-// 'GFace *to' is the top surface to mesh, 'from' is the source surface, 'pos' is a std::set
-// of vertex positions for the top surface.
-int MeshQuadToTriTopSurface( GFace *from, GFace *to, std::set<MVertex*,
- MVertexLessThanLexicographic> &pos );
+// This function meshes the top surface of a QuadToTri extrusion. It returns 0
+// if it is given a non-quadToTri extrusion or if it fails. Args: 'GFace *to'
+// is the top surface to mesh, 'from' is the source surface, 'pos' is a tree
+// of vertex positions for the top surface.
+int MeshQuadToTriTopSurface(GFace *from, GFace *to, MVertexRTree &pos);
#endif
diff --git a/Mesh/QuadTriExtruded3D.cpp b/Mesh/QuadTriExtruded3D.cpp
index cab245c3022a0b455742b793ed89a328e294f9a8..39f4fd3c5c392c732fafd0d82c95844fc728f3c9 100644
--- a/Mesh/QuadTriExtruded3D.cpp
+++ b/Mesh/QuadTriExtruded3D.cpp
@@ -1,36 +1,37 @@
-/**************************************************************************************************
+/*******************************************************************************
QuadTriExtruded3D.cpp
The code in this file was written by Dr. Trevor S. Strickler.
email: <trevor.strickler@gmail.com>
-This file is part of the QuadTri contribution to Gmsh. QuadTri allows the conformal interface
-of quadrangle faces to triangle faces using pyramids and other mesh elements.
+This file is part of the QuadTri contribution to Gmsh. QuadTri allows the
+conformal interface of quadrangle faces to triangle faces using pyramids and
+other mesh elements.
-See READMEQUADTRI.txt for more information. The license information is in LICENSE.txt.
+See READMEQUADTRI.txt for more information. The license information is in
+LICENSE.txt.
-Trevor S. Strickler hereby transfers copyright of QuadTri files to
-Christophe Geuzaine and J.-F. Remacle with the understanding that
-his contribution shall be cited appropriately.
+Trevor S. Strickler hereby transfers copyright of QuadTri files to Christophe
+Geuzaine and J.-F. Remacle with the understanding that his contribution shall be
+cited appropriately.
-All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
-Gmsh is available at: www.geuz.org/gmsh
+All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine,
+J.-F. Remacle Gmsh is available at: www.geuz.org/gmsh
-For Gmsh license information, see the LICENSE.txt file for license information. Please report all
-Gmsh bugs and problems to the public mailing list <gmsh@geuz.org>.
+For Gmsh license information, see the LICENSE.txt file for license
+information. Please report all Gmsh bugs and problems to the public mailing list
+<gmsh@geuz.org>.
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License, Version 2,
-as published by the Free Software Foundation, or (at your option)
-any later version, with or without the exception given in the
-LICENSE.txt file supplied with this code and with Gmsh.
+This program is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License, Version 2, as published by the Free
+Software Foundation, or (at your option) any later version, with or without the
+exception given in the LICENSE.txt file supplied with this code and with Gmsh.
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
+This program is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE. See the GNU General Public License for more details.
-****************************************************************************************************/
+********************************************************************************/
#include "QuadTriExtruded3D.h"
@@ -67,7 +68,6 @@ static void addHexahedron(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
to->hexahedra.push_back(newElem);
}
-
// Does the pair of MVertex pointers v1 and v2 exist in the set 'edges'?
static int edgeExists(MVertex *v1, MVertex *v2,
std::set<std::pair<MVertex*, MVertex*> > &edges)
@@ -76,7 +76,6 @@ static int edgeExists(MVertex *v1, MVertex *v2,
return edges.count(p);
}
-
// Create the pair of MVertex pointers v1 and v2 exist in the set 'edges.'
static void createEdge(MVertex *v1, MVertex *v2,
std::set<std::pair<MVertex*, MVertex*> > &edges)
@@ -85,11 +84,11 @@ static void createEdge(MVertex *v1, MVertex *v2,
edges.insert(p);
}
-
-// Create the entry for a forbidden edge in forbidden_edges (note that all four verts are
-// forbidden, but only store two, using lowest vertex pointer diagonal).
+// Create the entry for a forbidden edge in forbidden_edges (note that all four
+// verts are forbidden, but only store two, using lowest vertex pointer
+// diagonal).
static void createForbidden(std::vector<MVertex*> v,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges)
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges)
{
if( v.size() != 4 ){
Msg::Error("In createForbidden(), number of vertices not equal 4.");
@@ -105,9 +104,10 @@ static void createForbidden(std::vector<MVertex*> v,
}
-// Is the given vector of quad vertices forbidden to diagonalize (it is in forbidden_edges)?
+// Is the given vector of quad vertices forbidden to diagonalize (it is in
+// forbidden_edges)?
static int forbiddenExists(std::vector<MVertex*> v,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges)
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges)
{
if( v.size() != 4 ){
@@ -124,7 +124,6 @@ static int forbiddenExists(std::vector<MVertex*> v,
return forbidden_edges.count(pair);
}
-
// delete a pair of vertex pointers v1 and v2 from 'edges.'
static void deleteEdge(MVertex *v1, MVertex *v2,
std::set<std::pair<MVertex*, MVertex*> > &edges)
@@ -133,16 +132,15 @@ static void deleteEdge(MVertex *v1, MVertex *v2,
edges.erase(p);
}
-
-// Get the two mesh vertices extruded from vertices v0 and v1 on a lateral face at layer j, element k.
-// Added 2010-01-26
-static std::vector<MVertex*> getExtrudedLateralVertices( MVertex *v0, MVertex *v1, GEntity *entity,
- unsigned int j, unsigned int k, ExtrudeParams *loop_ep,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+// Get the two mesh vertices extruded from vertices v0 and v1 on a lateral face
+// at layer j, element k. Added 2010-01-26
+static std::vector<MVertex*> getExtrudedLateralVertices(MVertex *v0, MVertex *v1,
+ GEntity *entity,
+ unsigned int j, unsigned int k,
+ ExtrudeParams *loop_ep,
+ MVertexRTree &pos)
{
-
std::vector<MVertex*> verts;
- std::set<MVertex *, MVertexLessThanLexicographic>::iterator itp;
double x[4] = {v0->x(), v1->x(), v0->x(), v1->x()};
double y[4] = {v0->y(), v1->y(), v0->y(), v1->y()};
double z[4] = {v0->z(), v1->z(), v0->z(), v1->z()};
@@ -151,37 +149,28 @@ static std::vector<MVertex*> getExtrudedLateralVertices( MVertex *v0, MVertex *v
loop_ep->Extrude(j, k + 1, x[p + 2], y[p + 2], z[p + 2]);
}
for(int p = 0; p < 4; p++){
- MVertex tmp(x[p], y[p], z[p], 0, -1);
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()){
+ MVertex *tmp = pos.find(x[p], y[p], z[p]);
+ if(!tmp){
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g) in geometrical entity %d",
- tmp.x(), tmp.y(), tmp.z(), entity->tag());
-
+ x[p], y[p], z[p], entity->tag());
verts.clear();
return verts;
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
return verts;
-
}
-
-// Get the extruded vertices from MElement *elem at layer j, element k.
-// Added 2010-01-26
-static int get2DExtrudedVertices( MElement *elem, ExtrudeParams *ep, unsigned int j, unsigned int k,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
- std::vector<MVertex *> &verts )
+// Get the extruded vertices from MElement *elem at layer j, element k. Added
+// 2010-01-26
+static int get2DExtrudedVertices(MElement *elem, ExtrudeParams *ep,
+ unsigned int j, unsigned int k,
+ MVertexRTree &pos, std::vector<MVertex *> &verts)
{
std::vector<MVertex*> source_verts;
- elem->getVertices( source_verts );
+ elem->getVertices(source_verts);
- std::set<MVertex *, MVertexLessThanLexicographic>::iterator itp;
int sz = source_verts.size();
std::vector<double> x(sz), y(sz), z(sz);
for( int p = 0; p < sz; p++ ){
@@ -191,35 +180,24 @@ static int get2DExtrudedVertices( MElement *elem, ExtrudeParams *ep, unsigned in
}
for(int p = 0; p < sz; p++){
ep->Extrude(j, k, x[p], y[p], z[p]);
- MVertex tmp(x[p], y[p], z[p], 0, -1);
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()){
+ MVertex *tmp = pos.find(x[p], y[p], z[p]);
+ if(!tmp){
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g).",
- tmp.x(), tmp.y(), tmp.z() );
-
+ x[p], y[p], z[p]);
verts.clear();
return verts.size();
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
return verts.size();
-
}
-
-// Copied from meshGRegionExtruded.cpp, By Geuzaine, Remacle...
-// Extrudes a set of source vertices in 3D
-// added 2010-01-18
+// Copied from meshGRegionExtruded.cpp, By Geuzaine, Remacle... Extrudes a set
+// of source vertices in 3D added 2010-01-18
static int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
- std::vector<MVertex*> &verts)
+ MVertexRTree &pos, std::vector<MVertex*> &verts)
{
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
double x[8], y[8], z[8];
int n = ele->getNumVertices();
for(int p = 0; p < n; p++){
@@ -233,26 +211,19 @@ static int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
ep->Extrude(j, k + 1, x[p + n], y[p + n], z[p + n]);
}
for(int p = 0; p < 2 * n; p++){
- MVertex tmp(x[p], y[p], z[p], 0, -1);
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end())
+ MVertex *tmp = pos.find(x[p], y[p], z[p]);
+ if(!tmp)
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g)",
- tmp.x(), tmp.y(), tmp.z());
+ x[p], y[p], z[p]);
else
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
return verts.size();
}
-
// Determines whether the region is a valid QuadToTri region. Performs some
-// basic checks, including whether there is a valid top, valid source,
-// and that the surfaces serving as laterals are structured
-// Added 2010-12-30
+// basic checks, including whether there is a valid top, valid source, and that
+// the surfaces serving as laterals are structured Added 2010-12-30
bool IsValidQuadToTriRegion(GRegion *region, bool *allNonGlobalSharedLaterals)
{
ExtrudeParams *ep = region->meshAttributes.extrude;
@@ -263,22 +234,22 @@ bool IsValidQuadToTriRegion(GRegion *region, bool *allNonGlobalSharedLaterals)
GModel *model = region->model();
// find source face
- GFace *reg_source = model->getFaceByTag( std::abs( ep->geo.Source ) );
+ GFace *reg_source = model->getFaceByTag(std::abs(ep->geo.Source));
if( !reg_source ){
Msg::Error("In IsValidQuadToTriRegion(), could not find source face "
- "%d for region %d.", std::abs( ep->geo.Source ),
- region->tag() );
+ "%d for region %d.", std::abs(ep->geo.Source),
+ region->tag());
return false;
}
-
+
bool is_toroidal = IsInToroidalQuadToTri(reg_source);
GFace *root = findRootSourceFaceForFace(reg_source);
-
+
// Find a source surface. Then find a COPIED_ENTITY that is the top surface.
- // Then determine if all the laterals are either all quad or all triangle.
- // If shared laterals are all static (quad or non subdivide triangles),
- // set the allNonGlobalSharedLaterals argument to true.
- // If any lateral is unstructured, error.
+ // Then determine if all the laterals are either all quad or all triangle. If
+ // shared laterals are all static (quad or non subdivide triangles), set the
+ // allNonGlobalSharedLaterals argument to true. If any lateral is
+ // unstructured, error.
bool foundTop = false, foundSource = false,
foundNoStruct = false, foundRoot = false;
@@ -288,7 +259,7 @@ bool IsValidQuadToTriRegion(GRegion *region, bool *allNonGlobalSharedLaterals)
(*allNonGlobalSharedLaterals) = true;
-
+
for( it = faces.begin(); it != faces.end(); it++ ){
ExtrudeParams *face_tmp_ep = (*it)->meshAttributes.extrude;
if( (*it) == root )
@@ -309,7 +280,7 @@ bool IsValidQuadToTriRegion(GRegion *region, bool *allNonGlobalSharedLaterals)
else if( top_source_tmp == reg_source &&
!IsSurfaceALateralForRegion(region, *it) )
foundTop = true;
-
+
}
// This is a check to see if there are lateral surface triangles that need to be edged globally in subdivide operation
else if( IsSurfaceALateralForRegion(region, *it) ){
@@ -332,7 +303,7 @@ bool IsValidQuadToTriRegion(GRegion *region, bool *allNonGlobalSharedLaterals)
// if didn't find the copied entity, maybe this is toroidal and the top has been replaced
if( is_toroidal && !foundTop && foundRoot && root != reg_source )
foundTop = true;
-
+
// test for errors
bool detectConflict = false;
if( !foundTop ){
@@ -367,15 +338,15 @@ bool IsValidQuadToTriRegion(GRegion *region, bool *allNonGlobalSharedLaterals)
// Ordering of faces starts with the lateral face containing verts[0] and verts[1], then goes in order
// of increasing vertex index around element. Finally, the bottom, then the top.
// Added 2010-01-21
-static std::map<std::string, std::vector<int> > getFaceTypes(GRegion *gr, MElement *elem, int j, int k,
- std::vector<MVertex *> &verts,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos_src_edge,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
- std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
- std::vector<bool> &vert_bnd, std::vector<int> &nfix1,
- std::vector<int> &nfix2, std::vector<int> &nadj1,
- std::vector<int> &nadj2, std::vector<int> &free_flag )
+static std::map<std::string, std::vector<int> >
+getFaceTypes(GRegion *gr, MElement *elem, int j, int k,
+ std::vector<MVertex *> &verts,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
+ std::vector<bool> &vert_bnd, std::vector<int> &nfix1,
+ std::vector<int> &nfix2, std::vector<int> &nadj1,
+ std::vector<int> &nadj2, std::vector<int> &free_flag )
{
std::map<std::string, std::vector<int> > face_types;
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -878,12 +849,11 @@ static void bruteForceEdgeQuadToTriPrism( GRegion *gr, MElement *elem,
return;
-} // end of bruteForceEdgeQuadToTriPrism()
-
+}
// Divide hexahedron degenerated at two points (degenerate face is a line) by brute force
-static void addEdgesForQuadToTriTwoPtDegenHexa( GRegion *gr, MElement *elem, ExtrudeParams *ep,
+static void addEdgesForQuadToTriTwoPtDegenHexa(GRegion *gr, MElement *elem, ExtrudeParams *ep,
int j, int k, std::vector<MVertex *> verts,
std::map<std::string, std::vector<int> > &face_types,
std::set<std::pair<MVertex*, MVertex*> > &edges_new,
@@ -893,7 +863,6 @@ static void addEdgesForQuadToTriTwoPtDegenHexa( GRegion *gr, MElement *elem, Ext
std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
std::vector<int> nfix1, std::vector<int> nfix2,
std::vector<int> nadj1, std::vector<int> nadj2,
std::vector<int> free_flag )
@@ -1204,7 +1173,7 @@ static void addEdgesForQuadToTriTwoPtDegenHexa( GRegion *gr, MElement *elem, Ext
return;
-} // end of addEdgesForQuadToTriTwoPtDegenHexa()
+}
// Divide a hexahedron degenerate at one point (one degenerate corner) by brute force.
@@ -1218,7 +1187,6 @@ static void addEdgesForQuadToTriOnePtDegenHexa( GRegion *gr, MElement *elem, Ext
std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
std::vector<int> nfix1, std::vector<int> nfix2,
std::vector<int> nadj1, std::vector<int> nadj2,
std::vector<int> free_flag )
@@ -1538,12 +1506,11 @@ static void addEdgesForQuadToTriOnePtDegenHexa( GRegion *gr, MElement *elem, Ext
problems_new[elem].insert(jkpair);
}
-} // end of addEdgesForQuadToTriOnePtDegenHexa()
-
+}
// Divide a fully non-degenerate hexahedron by brute force.
-static void addEdgesForQuadToTriFullHexa( GRegion *gr, MElement *elem, ExtrudeParams *ep,
+static void addEdgesForQuadToTriFullHexa(GRegion *gr, MElement *elem, ExtrudeParams *ep,
int j, int k, std::vector<MVertex *> verts,
std::map<std::string, std::vector<int> > &face_types,
std::set<std::pair<MVertex*, MVertex*> > &edges_new,
@@ -1553,10 +1520,9 @@ static void addEdgesForQuadToTriFullHexa( GRegion *gr, MElement *elem, ExtrudePa
std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
std::vector<int> nfix1, std::vector<int> nfix2,
std::vector<int> nadj1, std::vector<int> nadj2,
- std::vector<int> free_flag )
+ std::vector<int> free_flag)
{
// There are 4 main possibilities for minimum diags to guarantee slice:
@@ -2064,26 +2030,25 @@ static void addEdgesForQuadToTriFullHexa( GRegion *gr, MElement *elem, ExtrudePa
problems_new[elem].insert(jkpair);
}
-} // end of addEdgesForQuadToTriFullHexa();
+}
// Generate face diagonals to subdivide hexahedra by BRUTE FORCE. Not recommended for general use, but it
// is required for some elements which have all vertices on an external region boundary.
// Added 2010-01-29
-static void bruteForceEdgeQuadToTriHexa( GRegion *gr, MElement *elem,
- int j, int k, std::vector<MVertex *> verts,
- std::map<std::string, std::vector<int> > &face_types,
- std::set<std::pair<MVertex*, MVertex*> > &edges_new,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_new,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
- std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
- std::vector<int> nfix1, std::vector<int> nfix2,
- std::vector<int> nadj1, std::vector<int> nadj2,
- std::vector<int> free_flag )
+static void bruteForceEdgeQuadToTriHexa(GRegion *gr, MElement *elem,
+ int j, int k, std::vector<MVertex *> verts,
+ std::map<std::string, std::vector<int> > &face_types,
+ std::set<std::pair<MVertex*, MVertex*> > &edges_new,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_new,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ std::vector<int> nfix1, std::vector<int> nfix2,
+ std::vector<int> nadj1, std::vector<int> nadj2,
+ std::vector<int> free_flag)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -2220,10 +2185,10 @@ static void bruteForceEdgeQuadToTriHexa( GRegion *gr, MElement *elem,
// First shape: A PRISM
// Only two possibilities. Either the bottom can be divided along same diagonal as the top,
// or the one quad side can be divided joining the top diagonal. Otherwise, problem.
- if( face_types["degen"].size() ){
- addEdgesForQuadToTriTwoPtDegenHexa( gr, elem, ep, j, k, verts, face_types, edges_new, forbidden_new,
- quadToTri_edges, forbidden_edges, lat_tri_diags, problems_new,
- problems, pos, nfix1, nfix2, nadj1, nadj2, free_flag );
+ if(face_types["degen"].size()){
+ addEdgesForQuadToTriTwoPtDegenHexa(gr, elem, ep, j, k, verts, face_types, edges_new, forbidden_new,
+ quadToTri_edges, forbidden_edges, lat_tri_diags, problems_new,
+ problems, nfix1, nfix2, nadj1, nadj2, free_flag);
return;
}
@@ -2231,9 +2196,9 @@ static void bruteForceEdgeQuadToTriHexa( GRegion *gr, MElement *elem,
// Since this shape can be divided to create legitimate tets and/or pyramids,
// let us try it.
if( !face_types["degen"].size() && face_types["single_tri"].size() == 2 ){
- addEdgesForQuadToTriOnePtDegenHexa( gr, elem, ep, j, k, verts, face_types, edges_new, forbidden_new,
- quadToTri_edges, forbidden_edges, lat_tri_diags, problems_new,
- problems, pos, nfix1, nfix2, nadj1, nadj2, free_flag );
+ addEdgesForQuadToTriOnePtDegenHexa(gr, elem, ep, j, k, verts, face_types, edges_new, forbidden_new,
+ quadToTri_edges, forbidden_edges, lat_tri_diags, problems_new,
+ problems, nfix1, nfix2, nadj1, nadj2, free_flag);
return;
}
@@ -2243,30 +2208,27 @@ static void bruteForceEdgeQuadToTriHexa( GRegion *gr, MElement *elem,
// 2. One PAIR of opposite diagonals, ABSOLUTELY ALL OTHER FACES FORBIDDEN
// 2. Two PAIRS of opposite diags
// 3. One PAIR of opposite diags PLUS a vertex with two ADDITIONAL diags meeting on it.
- if( !face_types["single_tri"].size() && !face_types["degen"].size() ){
-
- addEdgesForQuadToTriFullHexa( gr, elem, ep, j, k, verts, face_types, edges_new, forbidden_new,
- quadToTri_edges, forbidden_edges, lat_tri_diags, problems_new,
- problems, pos, nfix1, nfix2, nadj1, nadj2, free_flag );
+ if(!face_types["single_tri"].size() && !face_types["degen"].size()){
+ addEdgesForQuadToTriFullHexa(gr, elem, ep, j, k, verts, face_types, edges_new, forbidden_new,
+ quadToTri_edges, forbidden_edges, lat_tri_diags, problems_new,
+ problems, nfix1, nfix2, nadj1, nadj2, free_flag);
return;
}
-
-} // end of bruteForceEdgeQuadToTriHexa()
-
+}
// This is a shortcut function to simply copy face diagonals all the way up a vertical column of extruded elements.
// This function may not save very many operations, but it is going to stay...
-static int ExtrudeDiags( GRegion *gr, std::vector<MVertex*> v, unsigned int j_start,
- unsigned int k_start, unsigned int j_top, unsigned int k_top,
- MElement *elem, ExtrudeParams *loop_ep,
- std::set<std::pair<MVertex*, MVertex*> > &edges_new,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_new,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static int ExtrudeDiags(GRegion *gr, std::vector<MVertex*> v, unsigned int j_start,
+ unsigned int k_start, unsigned int j_top, unsigned int k_top,
+ MElement *elem, ExtrudeParams *loop_ep,
+ std::set<std::pair<MVertex*, MVertex*> > &edges_new,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_new,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ MVertexRTree &pos)
{
if( !loop_ep || !loop_ep->mesh.QuadToTri || !loop_ep->mesh.ExtrudeMesh ){
@@ -2366,10 +2328,10 @@ static int ExtrudeDiags( GRegion *gr, std::vector<MVertex*> v, unsigned int j_st
// Fixed edges include top surface diagonals and any lateral surface diagonals that cannot be swapped.
// Added 2010-01-24
static bool QuadToTriGetRegionDiags(GRegion *gr,
- std::set<std::pair<MVertex*,MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*,MVertex*> > &forbidden_edges,
- std::set<std::pair<MVertex*,MVertex*> > &lat_tri_diags,
- std::set<MVertex *, MVertexLessThanLexicographic> &pos)
+ std::set<std::pair<MVertex*,MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*,MVertex*> > &forbidden_edges,
+ std::set<std::pair<MVertex*,MVertex*> > &lat_tri_diags,
+ MVertexRTree &pos)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -2400,7 +2362,7 @@ static bool QuadToTriGetRegionDiags(GRegion *gr,
bool is_toroidal = IsInToroidalQuadToTri(reg_source);
if( is_toroidal )
root_face = findRootSourceFaceForFace(reg_source);
-
+
for( it = faces.begin(); it != faces.end(); it++ ){
ExtrudeParams *face_tmp_ep = (*it)->meshAttributes.extrude;
if( (*it) == root_face )
@@ -2451,26 +2413,27 @@ static bool QuadToTriGetRegionDiags(GRegion *gr,
// take care of forbidden edges
// test whether this surface is a lateral bounded by two quadtri regions.
- // if so, and the other is not already meshed,
- // then don't make these forbidden. This is worked out in a
+ // if so, and the other is not already meshed,
+ // then don't make these forbidden. This is worked out in a
// lateral remesh later
std::vector<GRegion *> adj_regions;
int numNeighbors = 0;
numNeighbors = GetNeighborRegionsOfFace((*it), adj_regions);
int ind_notcurrent = adj_regions[0] == gr ? 1 : 0;
- if( !( numNeighbors == 2 && adj_regions[0]->meshAttributes.extrude && adj_regions[1]->meshAttributes.extrude &&
- adj_regions[0]->meshAttributes.extrude->mesh.ExtrudeMesh && adj_regions[1]->meshAttributes.extrude->mesh.ExtrudeMesh &&
- adj_regions[0]->meshAttributes.extrude->geo.Mode == EXTRUDED_ENTITY && adj_regions[1]->meshAttributes.extrude->geo.Mode == EXTRUDED_ENTITY &&
- adj_regions[0]->meshAttributes.extrude->mesh.QuadToTri && adj_regions[1]->meshAttributes.extrude->mesh.QuadToTri &&
- IsSurfaceALateralForRegion(adj_regions[ind_notcurrent], *it) &&
- !adj_regions[ind_notcurrent]->getNumMeshElements() ) ){
- for( unsigned int i = 0; i < (*it)->quadrangles.size(); i++){
+ if(!(numNeighbors == 2 && adj_regions[0]->meshAttributes.extrude && adj_regions[1]->meshAttributes.extrude &&
+ adj_regions[0]->meshAttributes.extrude->mesh.ExtrudeMesh && adj_regions[1]->meshAttributes.extrude->mesh.ExtrudeMesh &&
+ adj_regions[0]->meshAttributes.extrude->geo.Mode == EXTRUDED_ENTITY &&
+ adj_regions[1]->meshAttributes.extrude->geo.Mode == EXTRUDED_ENTITY &&
+ adj_regions[0]->meshAttributes.extrude->mesh.QuadToTri && adj_regions[1]->meshAttributes.extrude->mesh.QuadToTri &&
+ IsSurfaceALateralForRegion(adj_regions[ind_notcurrent], *it) &&
+ !adj_regions[ind_notcurrent]->getNumMeshElements())){
+ for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++){
std::vector<MVertex*> v;
(*it)->quadrangles[i]->getVertices(v);
createForbidden(v, forbidden_edges );
}
}
-
+
// at this point, if there are no triangles, continue
if( !(*it)->triangles.size() )
continue;
@@ -2592,7 +2555,7 @@ static bool QuadToTriGetRegionDiags(GRegion *gr,
}
else if( !(*it)->quadrangles.size() )
Msg::Error("In QuadToTriGetRegionDiags(), failed to find a diagonal in lateral surface %d.", (*it)->tag() );
-
+
index_guess += 2;
/*
@@ -2620,14 +2583,14 @@ static bool QuadToTriGetRegionDiags(GRegion *gr,
else
createEdge( elemEdge_tmp.first, elemEdge_tmp.second, lat_tri_diags );
*/
-
+
}
}
}
}
}
-
+
// Insert diagonals of the top surface into quadToTri_edges;
unsigned int index_guess = reg_source->triangles.size();
if( reg_top->quadrangles.size() && !is_toroidal ){
@@ -2636,27 +2599,25 @@ static bool QuadToTriGetRegionDiags(GRegion *gr,
return false;
}
- for( unsigned int i = 0; i < reg_source->quadrangles.size(); i++){
+ for(unsigned int i = 0; i < reg_source->quadrangles.size(); i++){
int j_top = ep->mesh.NbLayer-1;
int k_top = ep->mesh.NbElmLayer[ep->mesh.NbLayer-1];
-
MElement *elem = reg_source->quadrangles[i];
std::vector<MVertex *> verts;
- get2DExtrudedVertices( elem, ep, j_top, k_top, pos, verts );
- if( verts.size() != 4 ) break;
- if( !is_toroidal ){
+ get2DExtrudedVertices(elem, ep, j_top, k_top, pos, verts);
+ if(verts.size() != 4) break;
+ if(!is_toroidal){
// Find diagonal:
std::pair<int,int> diag(0,0);
- diag = FindDiagonalEdgeIndices( verts, reg_top, false, index_guess );
- if( diag.first || diag.second )
- createEdge( verts[diag.first], verts[diag.second], quadToTri_edges );
+ diag = FindDiagonalEdgeIndices(verts, reg_top, false, index_guess);
+ if(diag.first || diag.second)
+ createEdge( verts[diag.first], verts[diag.second], quadToTri_edges);
else
- Msg::Error("In QuadToTriGetRegionDiags(), failed to find a diagonal on top surface %d, but should have.", reg_top->tag() );
+ Msg::Error("In QuadToTriGetRegionDiags(), failed to find a diagonal on top surface %d, but should have.", reg_top->tag() );
index_guess += 2;
}
else
createForbidden(verts, forbidden_edges);
-
}
return true;
}
@@ -2665,16 +2626,14 @@ static bool QuadToTriGetRegionDiags(GRegion *gr,
// For use in QuadToTriEdgeGenerator: Controls BRUTE FORCE edging of elements with ALL vertices
// on a lateral boundary surface.
// Added 04/08/2011
-static int makeEdgesForElemsWithAllVertsOnBnd( GRegion *gr, bool is_addverts,
- CategorizedSourceElements &cat_src_elems,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
- std::map<MElement*, std::set<std::pair<unsigned int,
- unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static int makeEdgesForElemsWithAllVertsOnBnd(GRegion *gr, bool is_addverts,
+ CategorizedSourceElements &cat_src_elems,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ MVertexRTree &pos)
{
-
ExtrudeParams *ep = gr->meshAttributes.extrude;
if( !ep || !ep->mesh.QuadToTri || !ep->mesh.ExtrudeMesh ){
@@ -2685,29 +2644,29 @@ static int makeEdgesForElemsWithAllVertsOnBnd( GRegion *gr, bool is_addverts,
}
GModel *model = gr->model();
- if( !model ){
+ if(!model){
Msg::Error("In makeEdgesForElemsWithAllVertsOnBnd(), invalid model for region "
"%d.", gr->tag() );
return 0;
}
- GFace *reg_source = model->getFaceByTag( std::abs( ep->geo.Source ) );
- if( !reg_source ){
+ GFace *reg_source = model->getFaceByTag(std::abs(ep->geo.Source));
+ if(!reg_source){
Msg::Error("In makeEdgesForElemsWithAllVertsOnBnd(), invalid source face for region "
- "%d.", gr->tag() );
+ "%d.", gr->tag());
return 0;
}
- if( gr != cat_src_elems.region || reg_source != cat_src_elems.source_face ){
- Msg::Error("In makeEdgesForElemsWithAllVertsOnBnd(), too many elements in the CategorizedSourceElements "
- " structure for given source face %d.", reg_source->tag() );
+ if(gr != cat_src_elems.region || reg_source != cat_src_elems.source_face){
+ Msg::Error("In makeEdgesForElemsWithAllVertsOnBnd(), too many elements in the "
+ "CategorizedSourceElements structure for given source face %d.",
+ reg_source->tag());
return 0;
}
// find edge verts of source face
- std::set<MVertex *, MVertexLessThanLexicographic> pos_src_edge;
- QuadToTriInsertFaceEdgeVertices(reg_source, pos_src_edge );
-
+ MVertexRTree pos_src_edge(1.e-12 * CTX::instance()->lc);
+ QuadToTriInsertFaceEdgeVertices(reg_source, pos_src_edge);
// while Loop to diagonalize 3-boundary point triangles and 4-boundary point quadrangles:
bool finish_all_tri = false;
@@ -2808,7 +2767,7 @@ static int makeEdgesForElemsWithAllVertsOnBnd( GRegion *gr, bool is_addverts,
}
if( !s ){
- face_types = getFaceTypes( gr, elem, j, k, verts, pos_src_edge, quadToTri_edges,
+ face_types = getFaceTypes( gr, elem, j, k, verts, quadToTri_edges,
forbidden_edges, lat_tri_diags, vert_bnd, nfix1, nfix2,
nadj1, nadj2, free_flag);
@@ -2854,14 +2813,14 @@ static int makeEdgesForElemsWithAllVertsOnBnd( GRegion *gr, bool is_addverts,
problems.erase( elem );
break;
}
- face_types = getFaceTypes( gr, elem, j, k, verts, pos_src_edge, quadToTri_edges,
- forbidden_edges, lat_tri_diags, vert_bnd, nfix1, nfix2,
- nadj1, nadj2, free_flag);
-
- bruteForceEdgeQuadToTriHexa( gr, elem, j, k, verts, face_types, edges_new,
- forbidden_new, quadToTri_edges, forbidden_edges,
- lat_tri_diags, problems_new, problems, pos, nfix1, nfix2,
- nadj1, nadj2, free_flag);
+ face_types = getFaceTypes(gr, elem, j, k, verts, quadToTri_edges,
+ forbidden_edges, lat_tri_diags, vert_bnd, nfix1, nfix2,
+ nadj1, nadj2, free_flag);
+
+ bruteForceEdgeQuadToTriHexa(gr, elem, j, k, verts, face_types, edges_new,
+ forbidden_new, quadToTri_edges, forbidden_edges,
+ lat_tri_diags, problems_new, problems, nfix1, nfix2,
+ nadj1, nadj2, free_flag);
}
//IMPORTANT
num_edged++;
@@ -2903,21 +2862,18 @@ static int makeEdgesForElemsWithAllVertsOnBnd( GRegion *gr, bool is_addverts,
return 1;
-} // end of makeEdgesForElemsWithAllVertsOnBnd
-
-
+}
// For use in QuadToTriEdgeGenerator: Does the edging of prisms with some but not all vertices
// on a lateral boundary surface.
// Added 04/08/2011
-static int makeEdgesForOtherBndPrisms( GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
- std::map<MElement*, std::set<std::pair<unsigned int,
- unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static int makeEdgesForOtherBndPrisms(GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ MVertexRTree &pos)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -3031,20 +2987,18 @@ static int makeEdgesForOtherBndPrisms( GRegion *gr, bool is_addverts, Categorize
return 1;
-} // end of makeEdgesForOtherBndPrisms()
-
+}
// For use in QuadToTriEdgeGenerator: Does the edging of hexahedra with some but not all vertices
// on a lateral boundary surface.
// Added 04/08/2011
-static int makeEdgesForOtherBndHexa( GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
- std::map<MElement*, std::set<std::pair<unsigned int,
- unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static int makeEdgesForOtherBndHexa(GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ MVertexRTree &pos)
{
// Edge creation for extruded quadrangles with some but not all vertices on a boundary.
@@ -3276,9 +3230,9 @@ static int makeEdgesForOtherBndHexa( GRegion *gr, bool is_addverts, CategorizedS
createEdge( verts[skip+elem_size], verts[(skip+2)%elem_size+elem_size], quadToTri_edges );
createEdge( verts[skip+elem_size], verts[(skip+2)%elem_size+elem_size], edges_new );
// copies the diags in edges_new on up the extrusion but NOT on the top surface
- ExtrudeDiags( gr, verts, j_start, k_start, j_top_start, k_top_start,
- elem, ep, edges_new, forbidden_new,
- quadToTri_edges, forbidden_edges, problems_new, problems, pos);
+ ExtrudeDiags(gr, verts, j_start, k_start, j_top_start, k_top_start,
+ elem, ep, edges_new, forbidden_new,
+ quadToTri_edges, forbidden_edges, problems_new, problems, pos);
}
// finally, top diagonal and other two laterals for 1 bnd point quad
@@ -3355,18 +3309,18 @@ static int makeEdgesForOtherBndHexa( GRegion *gr, bool is_addverts, CategorizedS
return 1;
-} // end of makeEdgesForOtherBndHexa()
+}
// For use in QuadToTriEdgeGenerator: Does the lateral edging of internal elements that touch
// a pivot vertex of a hexahedral element that has ONE SOURCE vertex on a lateral boundary surface.
// See inside function for a definition of the "pivot vertex."
// Added 04/08/2011
- static int makeEdgesForElemsTouchPivotVert( GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
- std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static int makeEdgesForElemsTouchPivotVert(GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
+ std::set<std::pair<MVertex*, MVertex*> > &forbidden_edges,
+ MVertexRTree &pos)
{
// Draw diagonals toward the "pivot vertex" of a hexahedron whose source quad has only one
@@ -3582,15 +3536,15 @@ static int makeEdgesForOtherBndHexa( GRegion *gr, bool is_addverts, CategorizedS
return 1;
-} // end of makeEdgesForElemsTouchPivotVert()
+}
// For use in QuadToTriEdgeGenerator: Does the lateral edging of internal elements in the top extrusion
// layer by lowest vertex pointer value in TOP FACE.
// Added 04/08/2011
-static int makeEdgesInternalTopLayer( GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static int makeEdgesInternalTopLayer(GRegion *gr, bool is_addverts, CategorizedSourceElements &cat_src_elems,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ MVertexRTree &pos)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -3673,7 +3627,7 @@ static int makeEdgesInternalTopLayer( GRegion *gr, bool is_addverts, Categorized
return 1;
-} // End of makeEdgesInternalTopLayer()
+}
// Generate the set of QuadToTri diagonal edges to subdivide elements,
@@ -3683,7 +3637,7 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
std::set<std::pair<MVertex*, MVertex*> > &lat_tri_diags,
std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+ MVertexRTree &pos)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -3702,7 +3656,7 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
return 0;
}
-
+
// number of extrusion layers
int num_layers = 0;
for( int p = 0; p < ep->mesh.NbLayer; p++ )
@@ -3725,8 +3679,8 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
// need for toroidal loop extrusions...top layer treated specially
bool is_toroidal = IsInToroidalQuadToTri(reg_source);
-
-
+
+
std::list<GFace *> reg_faces = gr->faces();
std::list<GFace *>::iterator itf = reg_faces.begin();
@@ -3744,9 +3698,9 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
foundRoot = true;
if( reg_top && (foundRoot || !is_toroidal) )
break;
-
+
}
-
+
if( is_toroidal && !reg_top && foundRoot && root != reg_source )
reg_top = root;
if( !reg_top ){
@@ -3760,7 +3714,7 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
// insert ALL fixed edges into quadToTri_edges, all forbidden edges on recombined quads
// into forbidden_edges, and insert into lat_tri_diags ALL lateral diagonal edges.
- QuadToTriGetRegionDiags( gr, quadToTri_edges, forbidden_edges, lat_tri_diags, pos );
+ QuadToTriGetRegionDiags( gr, quadToTri_edges, forbidden_edges, lat_tri_diags, pos);
/*unsigned int Rnum = gr->tag()-1;
std::vector<MVertex*> verts;
MElement *elem;
@@ -3827,21 +3781,21 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
// BRUTE FORCE diagonalization of elements with all vertices on a lateral boundary of region:
// This has to be done for all cases with such elements if
- if( !makeEdgesForElemsWithAllVertsOnBnd( gr, is_addverts, cat_src_elems, quadToTri_edges,
- lat_tri_diags, forbidden_edges, problems, pos ) ){
- Msg::Error("In QuadToTriEdgeGenerator(), failed to make edges for the elements in region %d "
- "with all vertices on a lateral boundary", gr->tag() );
- return 0;
+ if(!makeEdgesForElemsWithAllVertsOnBnd(gr, is_addverts, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, forbidden_edges, problems, pos)){
+ Msg::Error("In QuadToTriEdgeGenerator(), failed to make edges for the elements in region %d "
+ "with all vertices on a lateral boundary", gr->tag());
+ return 0;
}
// now do the "elegant" diagonalization of all the rest of the surface elements....
// Extrude source triangles that are on the source boundary edges and find any diagonals
- if( !makeEdgesForOtherBndPrisms( gr, is_addverts, cat_src_elems, quadToTri_edges,
- lat_tri_diags, forbidden_edges, problems, pos ) ){
+ if(!makeEdgesForOtherBndPrisms(gr, is_addverts, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, forbidden_edges, problems, pos)){
Msg::Error("In QuadToTriEdgeGenerator(), failed to make edges for the prism extrusions in region %d with "
- "source triangles having some but not all vertices on the boundary", gr->tag() );
+ "source triangles having some but not all vertices on the boundary", gr->tag());
return 0;
}
@@ -3854,8 +3808,8 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
// Edge creation for extruded quadrangles with some but not all vertices on a boundary.
- if( !makeEdgesForOtherBndHexa( gr, is_addverts, cat_src_elems, quadToTri_edges,
- lat_tri_diags, forbidden_edges, problems, pos ) ){
+ if(!makeEdgesForOtherBndHexa(gr, is_addverts, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, forbidden_edges, problems, pos)){
Msg::Error("In QuadToTriEdgeGenerator(), failed to make edges for the hexahedral extrusions in region %d with "
"source quads having some but not all vertices on the boundary", gr->tag() );
return 0;
@@ -3865,8 +3819,8 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
// Find diagonals for elements touching a "pivot vertex" of a hexa element that has
// a source quad with only one vertex on a lateral boundary (see inside makeEdgesForOtherBndHexa() and
// makeEdgesForElemsTouchingPivotVert() for details of "pivot vertex".
- if( !makeEdgesForElemsTouchPivotVert( gr, is_addverts, cat_src_elems, quadToTri_edges,
- lat_tri_diags, forbidden_edges, pos ) ){
+ if(!makeEdgesForElemsTouchPivotVert(gr, is_addverts, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, forbidden_edges, pos)){
Msg::Error("In QuadToTriEdgeGenerator(), failed to make edges for "
"the elements in region %d touching a \'pivot vertex\' of a "
"hexa element with source quad having one vertex on a boundary.", gr->tag() );
@@ -3875,20 +3829,20 @@ int QuadToTriEdgeGenerator(GRegion *gr, CategorizedSourceElements &cat_src_elem
// Mesh internal elements in the top layer (just add lateral diagonals for the
// Do this by lowest pointer in top surface
- if( !is_toroidal && !makeEdgesInternalTopLayer( gr, is_addverts, cat_src_elems, quadToTri_edges, pos ) ){
+ if(!is_toroidal && !makeEdgesInternalTopLayer(gr, is_addverts, cat_src_elems, quadToTri_edges, pos)){
Msg::Error("In QuadToTriEdgeGenerator(), failed to make internal edges "
- "in top extrusion layer of region %d.", gr->tag() );
+ "in top extrusion layer of region %d.", gr->tag());
return 0;
}
return 1;
-} // End of QuadToTriEdgeGenerator()
+}
// Remesh the lateral 2D faces of QuadToTri regions using edges in quadToTri_edges as contraints
// Added 2010-01-24
-static bool QuadToTriLateralRemesh( GRegion *gr, std::set<std::pair<MVertex*,MVertex*> > &quadToTri_edges )
+static bool QuadToTriLateralRemesh(GRegion *gr, std::set<std::pair<MVertex*,MVertex*> > &quadToTri_edges)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -3912,10 +3866,10 @@ static bool QuadToTriLateralRemesh( GRegion *gr, std::set<std::pair<MVertex*,MVe
// If shared laterals are all static (quad or non subdivide triangles),
// set the allStaticSharedLaterals argument to true.
// If any lateral is unstructured, error.
-
+
bool is_toroidal = IsInToroidalQuadToTri(reg_source);
GFace *root = findRootSourceFaceForFace(reg_source);
-
+
bool foundTop = false, foundRoot = false;
GFace *reg_top = NULL;
std::list<GFace *> faces = gr->faces();
@@ -3973,11 +3927,11 @@ static bool QuadToTriLateralRemesh( GRegion *gr, std::set<std::pair<MVertex*,MVe
// Adds the face- or body-center vertices needed for some QuadToTri elements
-static bool addBodyCenteredVertices( GRegion *to, CategorizedSourceElements &c,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- bool is_addverts, unsigned int lat_tri_diags_size,
- std::set<MVertex *, MVertexLessThanLexicographic> &pos )
+static bool addBodyCenteredVertices(GRegion *to, CategorizedSourceElements &c,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ bool is_addverts, unsigned int lat_tri_diags_size,
+ MVertexRTree &pos)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
@@ -4124,8 +4078,8 @@ static bool addBodyCenteredVertices( GRegion *to, CategorizedSourceElements &c,
// Meshes either a prism or a hexahedral set of mesh vertices with an internal vertex
// created here in the function.
// Added 2010-03-30
-static void MeshWithInternalVertex( GRegion *to, MElement *source, std::vector<MVertex *> v, std::vector<int> n1,
- std::vector<int> n2, std::set<MVertex *, MVertexLessThanLexicographic> &pos )
+static void MeshWithInternalVertex(GRegion *to, MElement *source, std::vector<MVertex *> v, std::vector<int> n1,
+ std::vector<int> n2, MVertexRTree &pos)
{
int v_size = v.size();
@@ -4148,24 +4102,18 @@ static void MeshWithInternalVertex( GRegion *to, MElement *source, std::vector<M
// find the internal vertex
std::vector<double> centroid = QtFindVertsCentroid(v);
- MVertex tmp(centroid[0], centroid[1], centroid[2], 0, -1);
// it's too dangerous to use the 'new' command in here even with body-centered vertices.
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()){
+ MVertex *tmp = pos.find(centroid[0], centroid[1], centroid[2]);
+ if(!tmp){
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g) in geometrical entity %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
+ centroid[0], centroid[1], centroid[2], to->tag());
Msg::Error("MeshWithInternalVertex() failed to find body-centered vertex.");
return;
}
- MVertex *v_int = (*itp);
+ MVertex *v_int = tmp;
// build all pyramids/tetra
for( int p = 0; p < n_lat; p++ ){
@@ -4213,161 +4161,15 @@ static void MeshWithInternalVertex( GRegion *to, MElement *source, std::vector<M
}
-// Meshes either a prism or a hexahedral set of mesh vertices with face centered vertices
-// Can pick top or bottom faces to have the vertex, or both, based on the top_flag and bottom_flag args.
-// created here in the code
-// Added 2010-04-05
-/*static void MeshWithFaceCenteredVertex( GRegion *to, MElement *source, std::vector<MVertex *> v, std::vector<int> n1,
- std::vector<int> n2, bool bottom_flag, bool top_flag,
- std::set<MVertex *, MVertexLessThanLexicographic> &pos )
-{
-
- int v_size = v.size();
- int n_lat_tmp;
- if( v_size == 6 )
- n_lat_tmp = 3;
- else if( v_size == 8 )
- n_lat_tmp = 4;
- else{
- Msg::Error("In MeshWithFaceCenteredVertex(), number of element vertices does not equal 6 or 8.");
- return;
- }
-
- const int n_lat = n_lat_tmp;
-
- if( !top_flag && !bottom_flag ){
- Msg::Error("In MeshWithFaceCenteredVertex(), neither the top nor the bottom face were selected for vertex creation.");
- return;
- }
-
- if( (n_lat == 3 && n1.size() != 3) || (n_lat == 4 && n2.size() != 6) ){
- Msg::Error("In MeshWithFaceCenteredVertex(), size of diagonal node vectors is not does not equal 3 or 6.");
- return;
- }
-
- std::vector<MVertex *> face_vertices;
- // create the face-centered vertices
- for( int s = 0; s < 2; s++ ){
- if( (!s && !bottom_flag) || (s && !top_flag) )
- continue;
- std::vector<MVertex *> v_face;
- v_face.assign(n_lat, (MVertex*)(NULL) );
- int add = !s ? 0 : n_lat;
- for( int t = 0; t < n_lat; t++ )
- v_face[t] = v[t+add];
-
- std::vector<double> centroid = QtFindVertsCentroid(v_face);
- MVertex tmp(centroid[0], centroid[1], centroid[2], 0, -1);
-
- // it's too dangerous to use the 'new' command in here on face-centered vertices.
-
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()){
- Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g) in geometrical entity %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
- Msg::Error("MeshWithFaceCenteredVertex() failed to find face-centered vertex.");
- return;
- }
-
- face_vertices.push_back(*itp);
-
- }
-
- // If just one face-centered vertex, make a pyramid/tetra (or base-divided pyramid )
- // with the face-centered vertex as the apex. Fill in
- // pyramids/tetra aroud it laterally.
- if( face_vertices.size() == 1 ){
- int base, add;
- if( top_flag ){
- base = n_lat;
- add = 0;
- }
- else{
- base = n_lat+1;
- add = n_lat;
- }
- if( n_lat == 4 && n1[base] >= 0 ){
- addTetrahedron( v[n1[base]], v[n2[base]], v[(n1[base]-add+1)%n_lat+add], face_vertices[0], to, source );
- addTetrahedron( v[n1[base]], v[n2[base]], v[(n1[base]-add+n_lat-1)%n_lat+add], face_vertices[0], to, source );
- }
- else if( n_lat == 4 )
- addPyramid( v[0+add], v[1+add], v[2+add], v[3+add], face_vertices[0], to, source );
- else
- addTetrahedron( v[0+add], v[1+add], v[2+add], face_vertices[0], to, source );
-
- for( int p = 0; p < n_lat; p++ ){
- int p2 = (p+1)%n_lat;
- if( v[p] == v[p+n_lat] && v[p2] == v[p2+n_lat] )
- continue;
- else if( v[p] == v[p+n_lat] || v[p2] == v[p2+n_lat] ){
- MVertex *v_dup = (v[p] == v[p+n_lat]) ? v[p] : v[p2];
- MVertex *v_non_dup = (v_dup == v[p]) ? v[p2] : v[p];
- MVertex *v_non_dup2 = (v_non_dup == v[p]) ? v[p+n_lat] : v[p2+n_lat];
- addTetrahedron( v_dup, v_non_dup, v_non_dup2, face_vertices[0], to, source );
- }
- else if( n1[p] >= 0 ){
- int add_2 = n1[p] < n_lat ? n_lat : -n_lat;
- addTetrahedron( v[n1[p]], v[n2[p]], v[n1[p]+add_2], face_vertices[0], to, source );
- addTetrahedron( v[n1[p]], v[n2[p]], v[n2[p]-add_2], face_vertices[0], to, source );
- }
- else
- addPyramid( v[p], v[p+n_lat], v[p2+n_lat], v[p2], face_vertices[0], to, source );
- }
- return;
- }
-
- // If two face-centered vertices, make prisms around the central axis formed by the face-centered vertices.
- // Divide the prisms on the two inside faces to make the validity of the subdivision independent of the outer
- // later surface states (diagonalized or not)
- else if( face_vertices.size() == 2 ){
- for( int p = 0; p < n_lat; p++ ){
- int p2 = (p+1)%n_lat;
- if( v[p] == v[p+n_lat] && v[p2] == v[p2+n_lat] ){
- addTetrahedron( v[p], v[p2], face_vertices[0], face_vertices[1], to, source );
- }
- else if( v[p] == v[p+n_lat] || v[p2] == v[p2+n_lat] ){
- MVertex *v_dup = (v[p] == v[p+n_lat]) ? v[p] : v[p2];
- MVertex *v_non_dup = (v_dup == v[p]) ? v[p2] : v[p];
- MVertex *v_non_dup2 = (v_non_dup == v[p]) ? v[p+n_lat] : v[p2+n_lat];
- addTetrahedron( v_non_dup, v_non_dup2, face_vertices[1], v_dup, to, source );
- addTetrahedron( v_non_dup, face_vertices[0], face_vertices[1], v_dup, to, source );
- }
- else{
- addTetrahedron( v[p], v[(p+1)%n_lat], face_vertices[1], face_vertices[0], to, source );
- if( n1[p] >= 0 ){
- int add_2 = n1[p] < n_lat ? n_lat : -n_lat;
- addTetrahedron( v[n1[p]], v[n2[p]], v[n1[p]+add_2], face_vertices[1], to, source );
- addTetrahedron( v[n1[p]], v[n2[p]], v[n2[p]-add_2], face_vertices[1], to, source );
- }
- else
- addPyramid( v[p], v[p+n_lat], v[(p+1)%n_lat+n_lat], v[(p+1)%n_lat], face_vertices[1], to, source );
- }
- }
- return;
- }
-
- // If execute to here, send error message
- Msg::Error("MeshWithFaceCenteredVertex() failed. on region %d", to->tag() );
-
-}
-
-*/
-
-
// Construct the elements that subdivide a prism (or degenerated prism) in a QuadToTri interface;
// Added 2010-01-24
static inline void QuadToTriPriPyrTet(std::vector<MVertex*> &v, GRegion *to, int j,
- int k, MElement* source,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
- unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool diag_search,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+ int k, MElement* source,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
+ unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool diag_search,
+ MVertexRTree &pos)
{
int dup[3];
int m = 0;
@@ -4422,7 +4224,7 @@ static inline void QuadToTriPriPyrTet(std::vector<MVertex*> &v, GRegion *to, int
// mesh with added internal body centered vertex
// is this prism part of a QuadToTri internal vertex extrusion and does it need to be extruded as such?
if( is_addverts && bnd_elem && found_diags ){
- MeshWithInternalVertex( to, source, v, n1, n2, pos );
+ MeshWithInternalVertex( to, source, v, n1, n2, pos);
return;
}
@@ -4486,32 +4288,28 @@ static inline void QuadToTriPriPyrTet(std::vector<MVertex*> &v, GRegion *to, int
" be divided without internal vertex, but was not previously detected as such. "
" This is a bug. Please Report.");
Msg::Error("j: %d, k: %d", j, k);
- QtMakeCentroidVertex( v, &(to->mesh_vertices), to, pos );
+ QtMakeCentroidVertex( v, &(to->mesh_vertices), to, pos);
std::pair<unsigned int, unsigned int> jk_pair(j,k);
problems_new[source].insert(jk_pair);
is_problem = true;
}
if( is_problem ){
- MeshWithInternalVertex( to, source, v, n1, n2, pos );
+ MeshWithInternalVertex( to, source, v, n1, n2, pos);
return;
}
-
-
}
// Construct the elements that subdivde a two-point degenerated hexahedron (prism).
-static inline bool createTwoPtDegenHexElems( std::vector<MVertex*> &v, GRegion *to, ExtrudeParams *ep, int j,
- int k, int dup[], MElement* source, std::vector<int> n1, std::vector<int> n2,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
- unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool found_diags,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static inline bool createTwoPtDegenHexElems(std::vector<MVertex*> &v, GRegion *to, ExtrudeParams *ep, int j,
+ int k, int dup[], MElement* source, std::vector<int> n1, std::vector<int> n2,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
+ unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool found_diags)
{
-
if( !ep )
return 0;
@@ -4596,18 +4394,16 @@ static inline bool createTwoPtDegenHexElems( std::vector<MVertex*> &v, GRegion *
return 0;
-} // end of createTwoPtDegenHexElems()
-
+}
// Construct the elements that subdivide a one-point degenerated hexahedron extrusion
-static inline bool createOnePtDegenHexElems( std::vector<MVertex*> &v, GRegion *to, ExtrudeParams *ep, int j,
- int k, int dup[], MElement* source, std::vector<int> n1, std::vector<int> n2,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
- unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool found_diags,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static inline bool createOnePtDegenHexElems(std::vector<MVertex*> &v, GRegion *to, ExtrudeParams *ep, int j,
+ int k, int dup[], MElement* source, std::vector<int> n1, std::vector<int> n2,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
+ unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool found_diags)
{
if( !ep )
@@ -4838,23 +4634,20 @@ static inline bool createOnePtDegenHexElems( std::vector<MVertex*> &v, GRegion *
return 0;
-} // end of createOnePtDegenHexElems()
+}
// Construct the elements that subdivide a full hexahedron extrusion.
-static inline bool createFullHexElems( std::vector<MVertex*> &v, GRegion *to, ExtrudeParams *ep, int j,
- int k, int dup[], MElement* source, std::vector<int> n1, std::vector<int> n2,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
- unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool found_diags,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+static inline bool createFullHexElems(std::vector<MVertex*> &v, GRegion *to, ExtrudeParams *ep, int j,
+ int k, int dup[], MElement* source, std::vector<int> n1, std::vector<int> n2,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
+ unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool found_diags)
{
-
if( !ep )
return 0 ;
-
// First: does this hexa have ANY dividing diagonals? If no, return
if( !found_diags ){
addHexahedron( v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], to, source );
@@ -5236,17 +5029,17 @@ static inline bool createFullHexElems( std::vector<MVertex*> &v, GRegion *to, Ex
return 0; // if exhaust possibilities, default to this
-} // end of createFullHexElems()
+}
// Overall function that creates the elements that subdivide any whole element extruded from a quadrangle.
static inline void QuadToTriHexPri(std::vector<MVertex*> &v, GRegion *to, int j,
- int k, MElement* source,
- std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
- std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
- unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool diag_search,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+ int k, MElement* source,
+ std::set<std::pair<MVertex*, MVertex*> > &quadToTri_edges,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems,
+ std::map<MElement*, std::set<std::pair<unsigned int, unsigned int> > > &problems_new,
+ unsigned int lat_tri_diags_size, bool bnd_elem, bool is_addverts, bool diag_search,
+ MVertexRTree &pos)
{
int dup[4];
@@ -5270,7 +5063,7 @@ static inline void QuadToTriHexPri(std::vector<MVertex*> &v, GRegion *to, int j,
ExtrudeParams *ep = to->meshAttributes.extrude;
-
+
// variables to hold of each faces's diagonal vertex nodes
bool found_diags = false;
std::vector<int> n1, n2;
@@ -5315,7 +5108,7 @@ static inline void QuadToTriHexPri(std::vector<MVertex*> &v, GRegion *to, int j,
// Divide by new internal vertex extrusion method?
if( is_addverts && ( found_diags || m==1 ) ){
- MeshWithInternalVertex( to, source, v, n1, n2, pos );
+ MeshWithInternalVertex( to, source, v, n1, n2, pos);
return;
}
@@ -5323,46 +5116,46 @@ static inline void QuadToTriHexPri(std::vector<MVertex*> &v, GRegion *to, int j,
// The of the possibilites are for a 'no new vertex' extrusion
// PRISM
- else if( m == 2 && !is_problem ){
- if( createTwoPtDegenHexElems( v, to, ep, j, k, dup, source, n1, n2, quadToTri_edges, problems,
- problems_new, lat_tri_diags_size, bnd_elem, is_addverts, found_diags, pos ) )
+ else if(m == 2 && !is_problem){
+ if(createTwoPtDegenHexElems(v, to, ep, j, k, dup, source, n1, n2, quadToTri_edges, problems,
+ problems_new, lat_tri_diags_size, bnd_elem, is_addverts, found_diags))
return;
}
// DEGENERATE HEXAHEDRON
- else if( m == 1 && !is_problem ){
- if( createOnePtDegenHexElems( v, to, ep, j, k, dup, source, n1, n2, quadToTri_edges, problems,
- problems_new, lat_tri_diags_size, bnd_elem, is_addverts, found_diags, pos ) )
+ else if(m == 1 && !is_problem){
+ if(createOnePtDegenHexElems(v, to, ep, j, k, dup, source, n1, n2, quadToTri_edges, problems,
+ problems_new, lat_tri_diags_size, bnd_elem, is_addverts, found_diags))
return;
}
// FULL HEXAHEDRON
- else if( !is_problem ){
- if( createFullHexElems( v, to, ep, j, k, dup, source, n1, n2, quadToTri_edges, problems,
- problems_new, lat_tri_diags_size, bnd_elem, is_addverts, found_diags, pos ) )
+ else if(!is_problem){
+ if(createFullHexElems(v, to, ep, j, k, dup, source, n1, n2, quadToTri_edges, problems,
+ problems_new, lat_tri_diags_size, bnd_elem, is_addverts, found_diags))
return;
}
// now take care of unexpected failure to divide without internal vertex
- if( !is_problem ) {
+ if(!is_problem) {
Msg::Error("In QuadToTriHexPri(), Extruded hexahedron needs subdivision, but cannot "
" be divided without internal vertex, and was not previously detected as such. "
" This is a bug. Please Report.");
Msg::Error("j: %d, k: %d", j, k);
- QtMakeCentroidVertex( v, &(to->mesh_vertices), to, pos );
+ QtMakeCentroidVertex(v, &(to->mesh_vertices), to, pos);
std::pair<unsigned int, unsigned int> jk_pair(j,k);
problems_new[source].insert(jk_pair);
is_problem = true;
}
// Mesh with internal vertex
- if( is_problem ){
- MeshWithInternalVertex( to, source, v, n1, n2, pos );
+ if(is_problem){
+ MeshWithInternalVertex(to, source, v, n1, n2, pos);
return;
}
-} // end of QuadToTriPriPyrTet()
+}
// reserves approximately the right amount of memory for quadToTri extrusions
// in the element vectors in the region 'to'
@@ -5412,10 +5205,10 @@ static void reserveQuadToTriCapacityForRegion( GRegion *to, GFace *from, bool i
*/
// displays for the user a list of the body centered vertices created for problem elements.
-static void listBodyCenteredVertices( GRegion *to, bool is_addverts,
- std::map<MElement *, std::set<std::pair<unsigned int,unsigned int> > > *problems,
- std::map<MElement *, std::set<std::pair<unsigned int,unsigned int> > > *problems_new,
- std::set<MVertex*, MVertexLessThanLexicographic> *pos )
+static void listBodyCenteredVertices(GRegion *to, bool is_addverts,
+ std::map<MElement *, std::set<std::pair<unsigned int,unsigned int> > > *problems,
+ std::map<MElement *, std::set<std::pair<unsigned int,unsigned int> > > *problems_new,
+ MVertexRTree *pos)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
@@ -5467,8 +5260,8 @@ static void listBodyCenteredVertices( GRegion *to, bool is_addverts,
std::set<std::pair<unsigned int, unsigned int> >::iterator itset;
for( itset = itmap->second.begin(); itset != itmap->second.end(); itset++ ){
std::vector<MVertex *> verts;
- getExtrudedVertices( itmap->first, ep, (*itset).first,
- (*itset).second, (*pos), verts);
+ getExtrudedVertices(itmap->first, ep, (*itset).first,
+ (*itset).second, (*pos), verts);
// find centroid
std::vector<double> centroid = QtFindVertsCentroid(verts);
int_verts_count2++;
@@ -5486,10 +5279,10 @@ static void listBodyCenteredVertices( GRegion *to, bool is_addverts,
// Function that makes all the elements in a QuadToTri region, both
// the divided elements and the whole elements, using already-created subdivision edges.
bool QuadToTriCreateElements(GRegion *to, CategorizedSourceElements &cat_src_elems,
- std::set<std::pair<MVertex*,MVertex*> > &quadToTri_edges,
- std::set<std::pair<MVertex*,MVertex*> > &lat_tri_diags,
- std::map<MElement*, std::set<std::pair<unsigned int,unsigned int> > > &problems,
- std::set<MVertex *, MVertexLessThanLexicographic> &pos)
+ std::set<std::pair<MVertex*,MVertex*> > &quadToTri_edges,
+ std::set<std::pair<MVertex*,MVertex*> > &lat_tri_diags,
+ std::map<MElement*, std::set<std::pair<unsigned int,unsigned int> > > &problems,
+ MVertexRTree &pos)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
@@ -5550,8 +5343,8 @@ bool QuadToTriCreateElements(GRegion *to, CategorizedSourceElements &cat_src_el
// Make the extra vertices needed for Some QuadToTri elements
- if( !addBodyCenteredVertices( to, cat_src_elems, quadToTri_edges, problems, is_addverts,
- lat_tri_diags_size, pos ) ){
+ if(!addBodyCenteredVertices(to, cat_src_elems, quadToTri_edges, problems, is_addverts,
+ lat_tri_diags_size, pos)){
Msg::Error("QuadToTriCreateElements() could not add face or body vertices for QuadToTri region %d.", to->tag() );
return false;
}
@@ -5590,10 +5383,10 @@ bool QuadToTriCreateElements(GRegion *to, CategorizedSourceElements &cat_src_el
for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
// keeps old allocation
verts.resize(0);
- if( getExtrudedVertices(elem, ep, j, k, pos, verts) == 6 ){
- QuadToTriPriPyrTet( verts, to, j, k, elem, quadToTri_edges,
+ if(getExtrudedVertices(elem, ep, j, k, pos, verts) == 6){
+ QuadToTriPriPyrTet(verts, to, j, k, elem, quadToTri_edges,
problems, problems_new, lat_tri_diags_size,
- bnd_elem, is_addverts, 1, pos );
+ bnd_elem, is_addverts, 1, pos);
}
}
}
@@ -5637,7 +5430,7 @@ bool QuadToTriCreateElements(GRegion *to, CategorizedSourceElements &cat_src_el
for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
// keeps old allocation
verts.resize(0);
- if(getExtrudedVertices(elem, ep, j, k, pos, verts) == 8 ){
+ if(getExtrudedVertices(elem, ep, j, k, pos, verts) == 8){
QuadToTriHexPri(verts, to, j, k, elem, quadToTri_edges,
problems, problems_new, lat_tri_diags_size,
bnd_elem, is_addverts, 1, pos);
@@ -5655,8 +5448,7 @@ bool QuadToTriCreateElements(GRegion *to, CategorizedSourceElements &cat_src_el
}
// List for the user any elements with internal vertices:
- listBodyCenteredVertices( to, is_addverts, &problems, &problems_new, &pos );
-
+ listBodyCenteredVertices(to, is_addverts, &problems, &problems_new, &pos);
// Now revert any elements that have positive volume.
// Does this even need to be done?
@@ -5699,7 +5491,7 @@ bool QuadToTriCreateElements(GRegion *to, CategorizedSourceElements &cat_src_el
// Mesh QuadToTri region from extrudeMesh() in meshGRegionExtruded.cpp
// Added 04/08/2011:
-int meshQuadToTriRegion( GRegion *gr, std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+int meshQuadToTriRegion(GRegion *gr, MVertexRTree &pos)
{
// Perform some checks to see if this is a valid QuadToTri region.
// If so, a decision has to be made: if this surface is NOT laterally adjacent to
@@ -5747,14 +5539,14 @@ int meshQuadToTriRegion( GRegion *gr, std::set<MVertex*, MVertexLessThanLexicogr
return 0;
}
- if( !QuadToTriEdgeGenerator( gr, cat_src_elems, quadToTri_edges,
- lat_tri_diags, problems, pos) ){
+ if(!QuadToTriEdgeGenerator(gr, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, problems, pos)){
Msg::Error("In meshQuadToTriRegion(), failed to create edges for QuadToTri "
"region %d.", gr->tag() );
return 0;
}
- if( !QuadToTriCreateElements( gr, cat_src_elems, quadToTri_edges,
- lat_tri_diags, problems, pos) ){
+ if(!QuadToTriCreateElements(gr, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, problems, pos)){
Msg::Error("In meshQuadToTriRegion, failed to create elements for QuadToTri "
"region %d.", gr->tag() );
return 0;
@@ -5773,8 +5565,8 @@ int meshQuadToTriRegion( GRegion *gr, std::set<MVertex*, MVertexLessThanLexicogr
// The function that is called from meshGRegionExtruded.cpp to mesh QuadToTri regions
// that are adjacent to subdivided regions, after the global Subdivide command is called.
// Added 04/08/11.
-int meshQuadToTriRegionAfterGlobalSubdivide( GRegion *gr, std::set<std::pair<MVertex*, MVertex*> > *edges,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+int meshQuadToTriRegionAfterGlobalSubdivide(GRegion *gr, std::set<std::pair<MVertex*, MVertex*> > *edges,
+ MVertexRTree &pos)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -5835,16 +5627,16 @@ int meshQuadToTriRegionAfterGlobalSubdivide( GRegion *gr, std::set<std::pair<MVe
}
// Mesh quadToTri
- if( !QuadToTriEdgeGenerator( gr, cat_src_elems, quadToTri_edges,
- lat_tri_diags, problems, pos ) ){
+ if(!QuadToTriEdgeGenerator(gr, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, problems, pos)){
Msg::Error("In meshQuadToTriRegionAfterGlobalSubdivide(), edge generation failed for QuadToTri "
"region %d.", gr->tag() );
return 0;
}
- if( !QuadToTriCreateElements( gr, cat_src_elems, quadToTri_edges,
- lat_tri_diags, problems, pos ) ){
+ if(!QuadToTriCreateElements(gr, cat_src_elems, quadToTri_edges,
+ lat_tri_diags, problems, pos)){
Msg::Error("In meshQuadToTriRegionAfterGlobalSubdivide(), element creation failed for QuadToTri "
- "region %d.", gr->tag() );
+ "region %d.", gr->tag());
return 0;
}
@@ -5853,4 +5645,4 @@ int meshQuadToTriRegionAfterGlobalSubdivide( GRegion *gr, std::set<std::pair<MVe
return 1;
-} // end of meshQuadToTriRegionAfterGlobalSubdivide()
+}
diff --git a/Mesh/QuadTriExtruded3D.h b/Mesh/QuadTriExtruded3D.h
index c7efd0e311884d29b95fb39cbd9194e657ca5e66..6bb19b6cb6e2df97231cf258015994a9d54ac7c3 100644
--- a/Mesh/QuadTriExtruded3D.h
+++ b/Mesh/QuadTriExtruded3D.h
@@ -1,36 +1,37 @@
-/**************************************************************************************************
+/********************************************************************************
QuadTriExtruded3D.h
The code in this file was written by Dr. Trevor S. Strickler.
-email: <trevor.strickler@gmail.com>
+email: <trevor.strickler@gmail.com>
-This file is part of the QuadTri contribution to Gmsh. QuadTri allows the conformal interface
-of quadrangle faces to triangle faces using pyramids and other mesh elements.
+This file is part of the QuadTri contribution to Gmsh. QuadTri allows the
+conformal interface of quadrangle faces to triangle faces using pyramids and
+other mesh elements.
-See READMEQUADTRI.txt for more information. The license information is in LICENSE.txt.
+See READMEQUADTRI.txt for more information. The license information is in
+LICENSE.txt.
-Trevor S. Strickler hereby transfers copyright of QuadTri files to
-Christophe Geuzaine and J.-F. Remacle with the understanding that
-his contribution shall be cited appropriately.
+Trevor S. Strickler hereby transfers copyright of QuadTri files to Christophe
+Geuzaine and J.-F. Remacle with the understanding that his contribution shall be
+cited appropriately.
-All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
-Gmsh is available at: www.geuz.org/gmsh
+All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine,
+J.-F. Remacle Gmsh is available at: www.geuz.org/gmsh
-For Gmsh license information, see the LICENSE.txt file for license information. Please report all
-Gmsh bugs and problems to the public mailing list <gmsh@geuz.org>.
+For Gmsh license information, see the LICENSE.txt file for license
+information. Please report all Gmsh bugs and problems to the public mailing list
+<gmsh@geuz.org>.
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License, Version 2,
-as published by the Free Software Foundation, or (at your option)
-any later version, with or without the exception given in the
-LICENSE.txt file supplied with this code and with Gmsh.
+This program is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License, Version 2, as published by the Free
+Software Foundation, or (at your option) any later version, with or without the
+exception given in the LICENSE.txt file supplied with this code and with Gmsh.
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
+This program is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE. See the GNU General Public License for more details.
-****************************************************************************************************/
+********************************************************************************/
#if !defined( _QTEXTR3D_H_ )
#define _QTEXTR3D_H_
@@ -59,21 +60,20 @@ GNU General Public License for more details.
#include <math.h>
#include "QuadTriUtils.h"
-// Determines whether the region is a valid QuadToTri region. Performs some
-// basic checks, including whether there is a valid top, valid source,
-// and that the surfaces serving as laterals are structured
-// Added 2010-12-30
+// Determines whether the region is a valid QuadToTri region. Performs some
+// basic checks, including whether there is a valid top, valid source, and that
+// the surfaces serving as laterals are structured Added 2010-12-30
bool IsValidQuadToTriRegion(GRegion *region, bool *allNonGlobalSharedLaterals);
-// Mesh QuadToTri region from extrudeMesh() from meshGRegionExtruded.cpp.
-// Added 04/08/2011:
-int meshQuadToTriRegion( GRegion *gr, std::set<MVertex*, MVertexLessThanLexicographic> &pos );
-
-// The function that is called from meshGRegionExtruded.cpp to mesh QuadToTri regions
-// that are adjacent to subdivided regions, after the global Subdivide command is called.
-// Added 04/08/11.
-int meshQuadToTriRegionAfterGlobalSubdivide( GRegion *gr, std::set<std::pair<MVertex*, MVertex*> > *edges,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos );
+// Mesh QuadToTri region from extrudeMesh() from meshGRegionExtruded.cpp. Added
+// 04/08/2011:
+int meshQuadToTriRegion(GRegion *gr, MVertexRTree &pos);
+// The function that is called from meshGRegionExtruded.cpp to mesh QuadToTri
+// regions that are adjacent to subdivided regions, after the global Subdivide
+// command is called. Added 04/08/11.
+int meshQuadToTriRegionAfterGlobalSubdivide(GRegion *gr,
+ std::set<std::pair<MVertex*, MVertex*> > *edges,
+ MVertexRTree &pos);
#endif
diff --git a/Mesh/QuadTriUtils.cpp b/Mesh/QuadTriUtils.cpp
index ff34f064fc15123371b989bf38d617b8f7f57c7d..055cdd785ef75faff5957a24d864ec150953c12f 100644
--- a/Mesh/QuadTriUtils.cpp
+++ b/Mesh/QuadTriUtils.cpp
@@ -1,65 +1,41 @@
-/**************************************************************************************************
+/*******************************************************************************
QuadTriUtils.cpp
The code in this file was written by Dr. Trevor S. Strickler.
email: <trevor.strickler@gmail.com>
-This file is part of the QuadTri contribution to Gmsh. QuadTri allows the conformal interface
-of quadrangle faces to triangle faces using pyramids and other mesh elements.
+This file is part of the QuadTri contribution to Gmsh. QuadTri allows the
+conformal interface of quadrangle faces to triangle faces using pyramids and
+other mesh elements.
-See READMEQUADTRI.txt for more information. The license information is in LICENSE.txt.
+See READMEQUADTRI.txt for more information. The license information is in
+LICENSE.txt.
-Trevor S. Strickler hereby transfers copyright of QuadTri files to
-Christophe Geuzaine and J.-F. Remacle with the understanding that
-his contribution shall be cited appropriately.
+Trevor S. Strickler hereby transfers copyright of QuadTri files to Christophe
+Geuzaine and J.-F. Remacle with the understanding that his contribution shall be
+cited appropriately.
-All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
-Gmsh is available at: www.geuz.org/gmsh
+All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine,
+J.-F. Remacle Gmsh is available at: www.geuz.org/gmsh
-For Gmsh license information, see the LICENSE.txt file for license information. Please report all
-Gmsh bugs and problems to the public mailing list <gmsh@geuz.org>.
+For Gmsh license information, see the LICENSE.txt file for license
+information. Please report all Gmsh bugs and problems to the public mailing list
+<gmsh@geuz.org>.
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License, Version 2,
-as published by the Free Software Foundation, or (at your option)
-any later version, with or without the exception given in the
-LICENSE.txt file supplied with this code and with Gmsh.
+This program is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License, Version 2, as published by the Free
+Software Foundation, or (at your option) any later version, with or without the
+exception given in the LICENSE.txt file supplied with this code and with Gmsh.
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
+This program is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE. See the GNU General Public License for more details.
-****************************************************************************************************/
+********************************************************************************/
#include <stdlib.h>
#include "QuadTriUtils.h"
-// This is a member function for the element map in ExtrudeParams.
-// This allows insertion of a whole vector at once.
-/*void ExtrudeParams::
-ExtrusionElementMap::addExtrudedElemVector(MElement* source, std::vector<MElement*> *extrudedVector )
-{
- std::map<MElement*,std::vector<MElement*> >::iterator it = _extrudedElements.find(source);
-
- if(it != _extrudedElements.end()){
- it->second.reserve(it->second.size()+extrudedVector->size());
- it->second.insert(it->second.end(), extrudedVector->begin(), extrudedVector->end());
- }
- else {
- int totalNbElems = 0;
- for (int i = 0; i <_parent->mesh.NbLayer;i++)
- totalNbElems += _parent->mesh.NbElmLayer[i];
- unsigned int new_cap = totalNbElems > extrudedVector->size() ? totalNbElems : extrudedVector->size();
- // automatically creates source key in map (single argument 'insert' is also logarithmic too, just like this).
- std::vector<MElement*> *vec = &(_extrudedElements[source]);
- vec->reserve( new_cap );
- vec->insert( vec->end(), extrudedVector->begin(), extrudedVector->end());
- }
-
-}*/
-
-
// By Geuzaine, Remacle...
static void addTriangle(MVertex* v1, MVertex* v2, MVertex* v3,
GFace *to)
@@ -68,17 +44,17 @@ static void addTriangle(MVertex* v1, MVertex* v2, MVertex* v3,
to->triangles.push_back(newTri);
}
-// this determines if a face is a non-lateral face in a structured toroidal volume extrusion with at
+// this determines if a face is a non-lateral face in a structured toroidal volume extrusion with at
// least one QuadToTri region...
int IsInToroidalQuadToTri(GFace *face)
{
- if( !face )
+ if(!face)
return false;
-
+
GModel *model = face->model();
-
+
bool is_toroidal = false, is_quadtri = false, is_noaddverts = false;
-
+
// Find the root face first...then step back through extrusions as far as can find
// another structured region. If there is a single quadtri region, and this is a torus that
// extrudes back onto the root source surface, then return true.
@@ -86,7 +62,7 @@ int IsInToroidalQuadToTri(GFace *face)
root_face = findRootSourceFaceForFace(face);
unsigned int numRegions = 0;
std::vector<GRegion*> adj_extruded_reg;
-
+
//find the two regions adjacent to the root face. If this is a structured torus, then both regions
// should be structured extrusions and BOTH should have the same root face
std::set<GRegion *, GEntityLessThan>::iterator itreg;
@@ -104,10 +80,10 @@ int IsInToroidalQuadToTri(GFace *face)
// does face belong to region and if so is it a structured extrusion?
std::list<GFace *> region_faces = std::list<GFace *>( region->faces() );
if( std::find( region_faces.begin(), region_faces.end(), root_face ) !=
- region_faces.end() && region->meshAttributes.extrude &&
+ region_faces.end() && region->meshAttributes.extrude &&
region->meshAttributes.extrude->mesh.ExtrudeMesh &&
region->meshAttributes.extrude->geo.Mode == EXTRUDED_ENTITY ){
-
+
adj_extruded_reg.push_back(region);
numRegions++;
}
@@ -115,9 +91,9 @@ int IsInToroidalQuadToTri(GFace *face)
continue;
}
// if there are two structured extruded regions adjacent to the root face,
- // then find the one that is NOT extruded from the root directly. Then follow this
+ // then find the one that is NOT extruded from the root directly. Then follow this
// face as far as possible back to the source.
-
+
GRegion *last_region = 0;
GFace *last_reg_source = 0;
bool found_first = 0, found_last = 0;
@@ -125,12 +101,12 @@ int IsInToroidalQuadToTri(GFace *face)
for( int ind = 0; ind <= 1; ind++ ){
ExtrudeParams *adj_ep = adj_extruded_reg[ind]->meshAttributes.extrude;
GFace *reg_source = 0;
-
+
if( adj_ep && adj_ep->mesh.ExtrudeMesh ){
reg_source = model->getFaceByTag(std::abs( adj_ep->geo.Source ) );
if( adj_ep->mesh.QuadToTri ){
is_quadtri = true;
- if( adj_ep->mesh.QuadToTri == QUADTRI_NOVERTS_1 ||
+ if( adj_ep->mesh.QuadToTri == QUADTRI_NOVERTS_1 ||
adj_ep->mesh.QuadToTri == QUADTRI_NOVERTS_1_RECOMB )
is_noaddverts = true;
}
@@ -153,7 +129,7 @@ int IsInToroidalQuadToTri(GFace *face)
}
}
- //walk back around to beginning if possible
+ //walk back around to beginning if possible
if( last_region && found_first && found_last ){
GFace *iter_face = last_reg_source;
GFace *iter_source_face = 0;
@@ -164,7 +140,7 @@ int IsInToroidalQuadToTri(GFace *face)
counter++;
if( iter_face ){
ExtrudeParams *iter_eps = iter_face->meshAttributes.extrude;
- if( iter_eps && iter_eps->mesh.ExtrudeMesh &&
+ if( iter_eps && iter_eps->mesh.ExtrudeMesh &&
iter_eps->geo.Mode == COPIED_ENTITY ){
if( iter_eps->mesh.QuadToTri )
is_quadtri = true;
@@ -192,7 +168,7 @@ int IsInToroidalQuadToTri(GFace *face)
}
}
}
-
+
// now return
if( is_toroidal && is_quadtri ){
if( !is_noaddverts )
@@ -209,37 +185,38 @@ void ReplaceBndQuadsInFace(GFace *face)
{
ExtrudeParams *ep = face->meshAttributes.extrude;
bool is_struct = false;
- if( (ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY) ||
- face->meshAttributes.method == MESH_TRANSFINITE )
+ if((ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY) ||
+ face->meshAttributes.method == MESH_TRANSFINITE)
is_struct = true;
- GFace *root_face = findRootSourceFaceForFace( face );
- if( root_face == face ){
- std::set<MVertex*, MVertexLessThanLexicographic> pos_src_edge;
+ GFace *root_face = findRootSourceFaceForFace(face);
+ if(root_face == face){
+ MVertexRTree pos_src_edge(1.e-12 * CTX::instance()->lc);
QuadToTriInsertFaceEdgeVertices(face, pos_src_edge);
- std::vector<MQuadrangle*> quads2;
+ std::vector<MQuadrangle*> quads2;
//loop through source quads, if on boundary, delete them
for(unsigned int i = 0; i < face->quadrangles.size(); i++){
std::vector<MVertex*> verts;
int num_verts = face->quadrangles[i]->getNumVertices();
bool on_edge = false;
-
- for( int j = 0; j < num_verts; j++ )
+
+ for(int j = 0; j < num_verts; j++)
verts.push_back(face->quadrangles[i]->getVertex(j));
-
- for( int j = 0; j < num_verts; j++ ){
- if( pos_src_edge.find( verts[j] ) != pos_src_edge.end() ){
+
+ for(int j = 0; j < num_verts; j++){
+ if(pos_src_edge.find(verts[j]->x(), verts[j]->y(), verts[j]->z())){
on_edge = true;
break;
}
}
- if( on_edge ){
+ if(on_edge){
delete face->quadrangles[i];
- if( is_struct ){
- addTriangle( verts[0],verts[2], verts[3],face);
- addTriangle( verts[0],verts[1], verts[2],face);
+ if(is_struct){
+ addTriangle(verts[0],verts[2], verts[3],face);
+ addTriangle(verts[0],verts[1], verts[2],face);
}
else{
- // draw other diagonals to minimize difference in average edge length with diagonal length, in quadrature
+ // draw other diagonals to minimize difference in average edge length
+ // with diagonal length, in quadrature
double mag_sq_ave = 0.0;
for( int p = 0; p < 4; p++ ){
@@ -267,15 +244,12 @@ void ReplaceBndQuadsInFace(GFace *face)
face->quadrangles.clear();
face->quadrangles = quads2;
}
-
-}
+}
-// Insert all vertices on a region's source edge, including corners,
-// into pos_src_edge set.
-// Added 2010-01-09
-void QuadToTriInsertSourceEdgeVertices(GRegion *gr,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos_src_edge)
+// Insert all vertices on a region's source edge, including corners, into
+// pos_src_edge set. Added 2010-01-09
+void QuadToTriInsertSourceEdgeVertices(GRegion *gr, MVertexRTree &pos_src_edge)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
if(!ep || !ep->mesh.ExtrudeMesh || ep->geo.Mode != EXTRUDED_ENTITY){
@@ -284,45 +258,35 @@ void QuadToTriInsertSourceEdgeVertices(GRegion *gr,
return;
}
-
- GFace *source_face = gr->model()->getFaceByTag( std::abs(ep->geo.Source) );
+ GFace *source_face = gr->model()->getFaceByTag(std::abs(ep->geo.Source));
std::list<GEdge*> edges = source_face->edges();
std::list<GEdge*>::iterator ite = edges.begin();
- for(ite = edges.begin(); ite != edges.end(); ite++ ){
- pos_src_edge.insert( (*ite)->mesh_vertices.begin(), (*ite)->mesh_vertices.end() );
- pos_src_edge.insert( (*ite)->getBeginVertex()->mesh_vertices.begin(),
- (*ite)->getBeginVertex()->mesh_vertices.end() );
- pos_src_edge.insert( (*ite)->getEndVertex()->mesh_vertices.begin(),
- (*ite)->getEndVertex()->mesh_vertices.end() );
+ for(ite = edges.begin(); ite != edges.end(); ite++){
+ pos_src_edge.insert((*ite)->mesh_vertices);
+ pos_src_edge.insert((*ite)->getBeginVertex()->mesh_vertices);
+ pos_src_edge.insert((*ite)->getEndVertex()->mesh_vertices);
}
}
-// Insert all vertices on a faces edges, including corners,
-// into pos_edges set.
+// Insert all vertices on a faces edges, including corners, into pos_edges set.
// Added 2010-01-18
-void QuadToTriInsertFaceEdgeVertices(GFace *face,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos_edges)
+void QuadToTriInsertFaceEdgeVertices(GFace *face, MVertexRTree &pos_edges)
{
-
-
std::list<GEdge*> edges = face->edges();
std::list<GEdge*>::iterator ite = edges.begin();
while(ite != edges.end()){
- pos_edges.insert((*ite)->mesh_vertices.begin(), (*ite)->mesh_vertices.end());
- pos_edges.insert((*ite)->getBeginVertex()->mesh_vertices.begin(),
- (*ite)->getBeginVertex()->mesh_vertices.end());
- pos_edges.insert((*ite)->getEndVertex()->mesh_vertices.begin(),
- (*ite)->getEndVertex()->mesh_vertices.end());
+ pos_edges.insert((*ite)->mesh_vertices);
+ pos_edges.insert((*ite)->getBeginVertex()->mesh_vertices);
+ pos_edges.insert((*ite)->getEndVertex()->mesh_vertices);
++ite;
}
}
-// Constructor for the CategorizedSourceElements structure.
-// See definition of CategorizedSourceElements in QuadTriUtils.h
-// file for details.
-CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
+// Constructor for the CategorizedSourceElements structure. See definition of
+// CategorizedSourceElements in QuadTriUtils.h file for details.
+CategorizedSourceElements::CategorizedSourceElements(GRegion *gr)
{
region = (GRegion*)(NULL);
source_face = (GFace*)(NULL);
@@ -330,7 +294,7 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
ExtrudeParams *ep = gr->meshAttributes.extrude;
- if( !ep || !ep->mesh.QuadToTri || !ep->mesh.ExtrudeMesh ){
+ if(!ep || !ep->mesh.QuadToTri || !ep->mesh.ExtrudeMesh){
Msg::Error("In CategorizedSourceElements constructor, invalid extrusion "
"in region %d for performing QuadToTri mesh generation.",
gr->tag() );
@@ -338,7 +302,7 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
}
GModel *model = gr->model();
- if( !model ){
+ if(!model){
Msg::Error("In CategorizedSourceElements constructor, invalid model for region "
"%d.", gr->tag() );
return;
@@ -347,7 +311,7 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
// now find and verify the source face
GFace *source_tmp = model->getFaceByTag( std::abs( ep->geo.Source ) );
- if( !source_tmp ){
+ if(!source_tmp){
Msg::Error("In CategorizedSourceElements constructor, invalid source face for region "
"%d.", gr->tag() );
return;
@@ -359,10 +323,9 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
source_face = source_tmp;
// get source face boundary verts
- std::set<MVertex*, MVertexLessThanLexicographic> bnd_verts;
+ MVertexRTree bnd_verts(1.e-12 * CTX::instance()->lc);
QuadToTriInsertSourceEdgeVertices(gr, bnd_verts);
-
unsigned int num_tri = source_face->triangles.size();
unsigned int num_quad = source_face->quadrangles.size();
@@ -370,18 +333,16 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
tri_bool.assign(num_tri*4, false);
quad_bool.assign(num_quad*5, false);
-
// keep temporary set of one boundary point quad pivot indices
std::set<MVertex*> one_pt_quad_pivots;
// now classify the source elements in the vectors
- for( int t = 0; t < 2; t++ ){ // t = 0 loop over triangles, t=1 loop over quads
-
+ for(int t = 0; t < 2; t++){ // t = 0 loop over triangles, t=1 loop over quads
int size = !t ? source_face->triangles.size() : source_face->quadrangles.size();
- for( int i = 0; i < size; i++ ){ // loop over total elements in current vector
+ for(int i = 0; i < size; i++){ // loop over total elements in current vector
std::vector<MVertex*> elem_verts;
MElement *elem;
- if( !t ){
+ if(!t){
elem = source_face->triangles[i];
elem->getVertices(elem_verts);
}
@@ -393,31 +354,31 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
int bnd_count = 0;
int bnd_vert = 0;
// get the boundary vert bool values
- for( int k = 0; k < elem_size; k++ ){
- if( bnd_verts.find(elem_verts[k]) != bnd_verts.end() ){
+ for(int k = 0; k < elem_size; k++){
+ if(bnd_verts.find(elem_verts[k]->x(), elem_verts[k]->y(), elem_verts[k]->z())){
if( !t ) tri_bool[4*i+k+1] = true;
else quad_bool[5*i+k+1] = true;
bnd_vert = k;
bnd_count++;
}
}
- if( bnd_count ){
- if( !t ) tri_bool[4*i] = true;
- else quad_bool[5*i] = true;
+ if(bnd_count){
+ if(!t) tri_bool[4*i] = true;
+ else quad_bool[5*i] = true;
}
// Place element vector indices into appropriate set
- if( !bnd_count )
+ if(!bnd_count)
(!t) ? internal_tri.insert(i) : internal_quad.insert(i);
- else if( bnd_count == 1 || bnd_count == 2 ||
- (bnd_count == 3 && t) )
+ else if(bnd_count == 1 || bnd_count == 2 ||
+ (bnd_count == 3 && t))
(!t) ? other_bnd_tri.insert(i) : other_bnd_quad.insert(i);
- else if( (bnd_count == 3 && !t) || bnd_count == 4 ){
+ else if((bnd_count == 3 && !t) || bnd_count == 4){
(!t) ? three_bnd_pt_tri.insert(i) : four_bnd_pt_quad.insert(i);
}
// if a one boundary point quad, record it in one_pt_quads set
- if( t && bnd_count == 1 ){
+ if(t && bnd_count == 1){
one_pt_quad_pivots.insert(elem_verts[(bnd_vert+2)%4]);
}
}
@@ -435,28 +396,28 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
for( int s = 0; s < 2; s++ ){
for( int q = 0; q < 2; q++ ){
std::set<unsigned int> *int_elems;
- if( !s ){
- if( !q ) int_elems = &internal_tri;
+ if(!s){
+ if(!q) int_elems = &internal_tri;
else continue;
}
else
int_elems = !q ? &internal_quad : &other_bnd_quad;
- for( it_int = (*int_elems).begin(); it_int != (*int_elems).end(); it_int++ ){
+ for(it_int = (*int_elems).begin(); it_int != (*int_elems).end(); it_int++){
std::vector<MVertex*> verts;
- if( !s )
+ if(!s)
source_face->triangles[(*it_int)]->getVertices(verts);
else
source_face->quadrangles[(*it_int)]->getVertices(verts);
// for t == s == 1 (other_bnd_quads), only retain those with two boundary verts
- if( s && q ){
+ if(s && q){
int bnd_count = 0;
- for(unsigned int k = 0; k < verts.size(); k++ ){
- if( quad_bool[5*(*it_int)+k+1] )
+ for(unsigned int k = 0; k < verts.size(); k++){
+ if(quad_bool[5*(*it_int)+k+1])
bnd_count++;
}
- if( bnd_count != 2 )
+ if(bnd_count != 2)
continue;
}
@@ -467,17 +428,17 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
// will sometimes choose the wrong pivot.
std::set<MVertex*>::iterator it_piv;
bool found = false;
- for( it_piv = one_pt_quad_pivots.begin(); it_piv != one_pt_quad_pivots.end(); it_piv++ ){
- for(unsigned int t = 0; t < verts.size(); t++ ){
- if( (*it_piv) == verts[t] ) {
+ for(it_piv = one_pt_quad_pivots.begin(); it_piv != one_pt_quad_pivots.end(); it_piv++){
+ for(unsigned int t = 0; t < verts.size(); t++){
+ if((*it_piv) == verts[t]) {
found = true;
- if( !s ){
+ if(!s){
internal_tri_touch_one_bnd_pt_quad.insert(*it_int);
tri_bool[4*(*it_int)] = true;
tri_bool[4*(*it_int)+t+1] = true;
}
else{
- if( !q ){
+ if(!q){
internal_quad_touch_one_bnd_pt_quad.insert(*it_int);
quad_bool[5*(*it_int)] = true;
quad_bool[5*(*it_int)+t+1] = true;
@@ -488,19 +449,17 @@ CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
break;
}
}
- if( found ) break;
+ if(found) break;
}
}
}
}
-} // end of CategorizedSourceElements::CategorizedSourceElements( GRegion *gr )
-
+}
// Find centroid of vertices in vector v, return in vector
-std::vector<double> QtFindVertsCentroid( std::vector<MVertex*> v )
+std::vector<double> QtFindVertsCentroid(std::vector<MVertex*> v)
{
-
std::vector<double> v_return;
const int v_size = v.size();
@@ -538,15 +497,13 @@ std::vector<double> QtFindVertsCentroid( std::vector<MVertex*> v )
v_return.push_back(z);
return v_return;
-
}
// Add a new vertex at the centroid of a vector of vertices (this goes into a region
// Added 2010-02-06
-MVertex* QtMakeCentroidVertex( std::vector<MVertex*> v, std::vector<MVertex*> *target,
- GEntity *entity, std::set<MVertex*, MVertexLessThanLexicographic> &pos )
+MVertex* QtMakeCentroidVertex(std::vector<MVertex*> v, std::vector<MVertex*> *target,
+ GEntity *entity, MVertexRTree &pos)
{
-
int v_size = v.size();
if( v_size != 6 && v_size != 8 && v_size != 3 && v_size != 4){
Msg::Error("In makeCentroidVertex(), number of vertices does not equal 3, 4, 6, or 8.");
@@ -563,28 +520,24 @@ MVertex* QtMakeCentroidVertex( std::vector<MVertex*> v, std::vector<MVertex*> *t
}
// make new vertex
- MVertex tmp(x, y, z, 0, -1);
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp = pos.find(&tmp);
+ MVertex *tmp = pos.find(x, y, z);
MVertex *v_int;
// simple check if it exists
- if(itp == pos.end()){
+ if(!tmp){
v_int = new MVertex(x, y, z, entity);
target->push_back(v_int);
pos.insert(v_int);
}
else
- v_int = (*itp);
+ v_int = tmp;
return (v_int);
-
}
-
// Finds the index of the lowest valued pointer in a vector of MVertex pointers
// Added 2011-03-10
-int getIndexForLowestVertexPointer( std::vector<MVertex*> v )
+int getIndexForLowestVertexPointer(std::vector<MVertex*> v)
{
-
int ind_low = 0;
int v_size = v.size();
for( int i = 1; i < v_size; i++ ){
@@ -594,7 +547,6 @@ int getIndexForLowestVertexPointer( std::vector<MVertex*> v )
return ind_low;
}
-
// Given 4 verts on a face, find an existent diagonal, if any.
// Two possible methods: If the 'index_guess' argument is the index of the correct triangle,
// finding it is simple. If not, have to do a complete pedantic search.
@@ -653,7 +605,6 @@ std::pair<int, int> FindDiagonalEdgeIndices( std::vector<MVertex*> verts,
}
-
// Get number of regions neighboring a face
int GetNeighborRegionsOfFace(GFace *face, std::vector<GRegion *> &neighbors)
{
@@ -698,7 +649,6 @@ int GetNeighborRegionsOfFace(GFace *face, std::vector<GRegion *> &neighbors)
// Trevor Strickler 12/09/10
int IsSurfaceALateralForRegion(GRegion *region, GFace *face)
{
-
// NOTE: don't necessarily require the face to be extruded! just in the position
// of lateral is all we care about here.
ExtrudeParams *ep = face->meshAttributes.extrude;
@@ -731,7 +681,6 @@ int IsSurfaceALateralForRegion(GRegion *region, GFace *face)
region_faces.end() )
return 0;
-
// if this face is a COPIED_ENTITY with source = region source face, this is the top. Exit.
if( ep && ep->geo.Mode == COPIED_ENTITY &&
reg_source == model->getFaceByTag(std::abs(ep->geo.Source)))
@@ -899,6 +848,3 @@ void fill_touch_bnd( int touch_bnd[], std::vector<bool> vert_bnd, int n_lat )
}
}
-
-
-
diff --git a/Mesh/QuadTriUtils.h b/Mesh/QuadTriUtils.h
index d3205b1fa8b9b31352e84059a07a4c7ae777f429..33f1c29cc444eb9b2f5474c3edda678fb098c5fa 100644
--- a/Mesh/QuadTriUtils.h
+++ b/Mesh/QuadTriUtils.h
@@ -1,16 +1,16 @@
-/**************************************************************************************************
+/**************************************************************************************************
QuadTriUtils.h
The code in this file was written by Dr. Trevor S. Strickler.
-email: <trevor.strickler@gmail.com>
+email: <trevor.strickler@gmail.com>
This file is part of the QuadTri contribution to Gmsh. QuadTri allows the conformal interface
-of quadrangle faces to triangle faces using pyramids and other mesh elements.
+of quadrangle faces to triangle faces using pyramids and other mesh elements.
See READMEQUADTRI.txt for more information. The license information is in LICENSE.txt.
-Trevor S. Strickler hereby transfers copyright of QuadTri files to
-Christophe Geuzaine and J.-F. Remacle with the understanding that
+Trevor S. Strickler hereby transfers copyright of QuadTri files to
+Christophe Geuzaine and J.-F. Remacle with the understanding that
his contribution shall be cited appropriately.
All reused or original Gmsh code is Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
@@ -21,8 +21,8 @@ Gmsh bugs and problems to the public mailing list <gmsh@geuz.org>.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License, Version 2,
-as published by the Free Software Foundation, or (at your option)
-any later version, with or without the exception given in the
+as published by the Free Software Foundation, or (at your option)
+any later version, with or without the exception given in the
LICENSE.txt file supplied with this code and with Gmsh.
This program is distributed in the hope that it will be useful,
@@ -45,6 +45,7 @@ GNU General Public License for more details.
#include "GModel.h"
#include "GmshDefines.h"
#include "MVertex.h"
+#include "MVertexRTree.h"
#include "Context.h"
#include "GModel.h"
#include "meshGFace.h"
@@ -68,9 +69,9 @@ GNU General Public License for more details.
// four_bnd_pt_quad : quadrangle vector indices for quads with four vertices on a surface boundary edge
// other_bnd_pt_tri, other_bnd_pt_quad : vector indices of the respective elements that have some
// but not all vertices on a boundary edge.
-// internal_quad_touch_one_bnd_pt_quad : quad vector indices of quads that share a vertex with a
+// internal_quad_touch_one_bnd_pt_quad : quad vector indices of quads that share a vertex with a
// one boundary point quad.
-// internal_tri_touch_one_bnd_pt_quad : triangle vector indices of triangles that share a vertex with a
+// internal_tri_touch_one_bnd_pt_quad : triangle vector indices of triangles that share a vertex with a
// one boundary point quad.
// two_pt_bnd_quad_touch_one_bnd_pt_quad : set of indices to quads with two points on a boundary that also touch quads
// with one point on a boundary.
@@ -82,22 +83,22 @@ GNU General Public License for more details.
// tri_bool, quad_bool : These vectors hold bits to tell which vertices in which elements are on the boundary, OR
// in the case of the [*]_touch_one_bnd_pt_quad elements, what vertex is the non-boundary 'pivot' vertex shared
// with a quad having one point on a boundary (diagonals are drawn to the pivot vertex in the single layer quadToTri
-// method).
-// Format: In tri_bool, there are 4*(number of triangles) elements. Each source triangle of index i has
+// method).
+// Format: In tri_bool, there are 4*(number of triangles) elements. Each source triangle of index i has
// four consecutive bool in tri_bool, beginning at i*4. The first bool is true if the triangle touches an edge
-// boundary (or touches a one boundary point quad). The other bools correspond to triangle vertices, in order
-// of appearance in the triangle's own vertex vector, and are true if the corresponding vertex is on the boundary.
+// boundary (or touches a one boundary point quad). The other bools correspond to triangle vertices, in order
+// of appearance in the triangle's own vertex vector, and are true if the corresponding vertex is on the boundary.
// Everything about tri_bool applies to quad_bool, but there are 5 bool per quad, accessed starting at i*5.
-//
+//
struct CategorizedSourceElements{
public:
- GRegion *region;
+ GRegion *region;
GFace *source_face;
bool valid;
- std::set<unsigned int> three_bnd_pt_tri, four_bnd_pt_quad,
+ std::set<unsigned int> three_bnd_pt_tri, four_bnd_pt_quad,
other_bnd_tri, other_bnd_quad;
- std::set<unsigned int> internal_tri_touch_one_bnd_pt_quad,
+ std::set<unsigned int> internal_tri_touch_one_bnd_pt_quad,
internal_quad_touch_one_bnd_pt_quad,
two_bnd_pt_quad_touch_one_bnd_pt_quad;
@@ -107,7 +108,7 @@ struct CategorizedSourceElements{
CategorizedSourceElements( GRegion *gr );
};
-// this determines if a face is a non-lateral face in a structured toroidal volume extrusion with at
+// this determines if a face is a non-lateral face in a structured toroidal volume extrusion with at
// least one QuadToTri region...
int IsInToroidalQuadToTri(GFace *face);
@@ -124,34 +125,32 @@ ExtrusionElementMap::addExtrudedElemVector(MElement* source, std::vector<MElemen
// Insert all vertices on a region's source edge, including corners,
// into pos_src_edge set.
// Added 2010-01-09
-void QuadToTriInsertSourceEdgeVertices(GRegion *gr,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos_src_edge);
+void QuadToTriInsertSourceEdgeVertices(GRegion *gr, MVertexRTree &pos_src_edge);
// Insert all vertices on a faces edges, including corners,
// into pos_edges set.
// Added 2010-01-18
-void QuadToTriInsertFaceEdgeVertices(GFace *face,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos_edges);
+void QuadToTriInsertFaceEdgeVertices(GFace *face, MVertexRTree &pos_edges);
// Find centroid of vertices in vector v, return in vector
-std::vector<double> QtFindVertsCentroid( std::vector<MVertex*> v );
+std::vector<double> QtFindVertsCentroid(std::vector<MVertex*> v);
// Add a new vertex at the centroid of a vector of vertices (this goes into a region
// Added 2010-02-06
-MVertex* QtMakeCentroidVertex( std::vector<MVertex*> v, std::vector<MVertex*> *target,
- GEntity *entity, std::set<MVertex*, MVertexLessThanLexicographic> &pos );
+MVertex* QtMakeCentroidVertex(std::vector<MVertex*> v, std::vector<MVertex*> *target,
+ GEntity *entity, MVertexRTree &pos);
// Finds the index of the lowest valued pointer in a vector of MVertex pointers
// Added 2011-03-10
-int getIndexForLowestVertexPointer( std::vector<MVertex*> v );
+int getIndexForLowestVertexPointer(std::vector<MVertex*> v);
// Given 4 verts on a face, find an existent diagonal, if any.
// Two possible methods: If the 'index_guess' argument is the index of the correct triangle,
// finding it is simple. If not, have to do a complete pedantic search.
// Added 2010-01-26
-std::pair<int, int> FindDiagonalEdgeIndices( std::vector<MVertex*> verts, GFace *face, bool lateral,
- unsigned int index_guess = 0 );
+std::pair<int, int> FindDiagonalEdgeIndices(std::vector<MVertex*> verts, GFace *face, bool lateral,
+ unsigned int index_guess = 0);
// Get number of regions neighboring a face
int GetNeighborRegionsOfFace(GFace *face, std::vector<GRegion *> &neighbors);
@@ -162,7 +161,7 @@ int IsSurfaceALateralForRegion(GRegion *region, GFace *face);
// Function to determine if a face is a top surface for a region. It returns 1
-// if the face is COPIED_ENTITY with source = region's source and if face belongs to region.
+// if the face is COPIED_ENTITY with source = region's source and if face belongs to region.
// Otherwise, return 0 (NOTE: ReplaceDuplicateSurfaces() can remove a top surface
// and replace it. If that happens, this will return 0. That is INTENDED for THIS function.
// Added 2010-12-13
@@ -173,12 +172,12 @@ int IsSurfaceATopForRegion(GRegion *region, GFace *face);
GFace* findRootSourceFaceForFace(GFace *face);
-// Input is vert_bnd[], which describes some 2D element: vert_bnd[i] is true if
-// the ith vertex the element touches a lateral edge boundary of the surface the
-// element is in.
-// Output is touch_bnd[]: Each element of touch_bnd[] corresponds to an edge of
+// Input is vert_bnd[], which describes some 2D element: vert_bnd[i] is true if
+// the ith vertex the element touches a lateral edge boundary of the surface the
+// element is in.
+// Output is touch_bnd[]: Each element of touch_bnd[] corresponds to an edge of
// the element described by vert_bnd[]. Edge i of touch_bnd[] is formed by
-// vertices i and (i+1)%element_size of the element. The value of touch_bnd[i] is non-zero
+// vertices i and (i+1)%element_size of the element. The value of touch_bnd[i] is non-zero
// if that edge touches a boundary edge of the surface that the element is in.
// Added 2011-03-10
void fill_touch_bnd( int touch_bnd[], std::vector<bool> vert_bnd, int n_lat );
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index 86a734c8c5b12f9a4030dc11dd6dd8161be79954..ca73da285093e62da0adb987fc3cc146c63520a0 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -9,6 +9,7 @@
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "ExtrudeParams.h"
+#include "MVertexRTree.h"
#include "Context.h"
#include "GmshMessage.h"
#include "QuadTriExtruded2D.h"
@@ -75,8 +76,7 @@ static void createQuaTri(std::vector<MVertex*> &v, GFace *to,
}
}
-static void extrudeMesh(GEdge *from, GFace *to,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+static void extrudeMesh(GEdge *from, GFace *to, MVertexRTree &pos,
std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
@@ -120,7 +120,6 @@ static void extrudeMesh(GEdge *from, GFace *to,
// create elements (note that it would be faster to access the
// *interior* nodes by direct indexing, but it's just simpler to
// query everything by position)
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
for(unsigned int i = 0; i < from->lines.size(); i++){
MVertex *v0 = from->lines[i]->getVertex(0);
MVertex *v1 = from->lines[i]->getVertex(1);
@@ -135,27 +134,21 @@ static void extrudeMesh(GEdge *from, GFace *to,
ep->Extrude(j, k + 1, x[p + 2], y[p + 2], z[p + 2]);
}
for(int p = 0; p < 4; p++){
- MVertex tmp(x[p], y[p], z[p], 0, -1);
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()){
+ MVertex *tmp = pos.find(x[p], y[p], z[p]);
+ if(!tmp){
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
+ x[p], y[p], z[p], to->tag());
return;
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
- createQuaTri(verts, to, constrainedEdges,from->lines[i], tri_quad_flag);
+ createQuaTri(verts, to, constrainedEdges, from->lines[i], tri_quad_flag);
}
}
}
}
-static void copyMesh(GFace *from, GFace *to,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+static void copyMesh(GFace *from, GFace *to, MVertexRTree &pos)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
@@ -198,58 +191,48 @@ static void copyMesh(GFace *from, GFace *to,
}
// create triangle elements
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
for(unsigned int i = 0; i < from->triangles.size(); i++){
std::vector<MVertex*> verts;
for(int j = 0; j < 3; j++){
MVertex *v = from->triangles[i]->getVertex(j);
- MVertex tmp(v->x(), v->y(), v->z(), 0, -1);
+ double x = v->x(), y = v->y(), z = v->z();
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
- tmp.x(), tmp.y(), tmp.z());
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()) {
+ x, y, z);
+ MVertex *tmp = pos.find(x, y, z);
+ if(!tmp) {
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
+ x, y, z, to->tag());
return;
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
addTriangle(verts[0], verts[1], verts[2], to);
}
// Add triangles for divided quads for QuadTri -- Trevor Strickler
// if quadtotri and not part of a toroidal extrusion, mesh the top surface accordingly
- if( detectQuadToTriTop && !is_toroidal ){
- if( !MeshQuadToTriTopSurface(from, to, pos))
+ if(detectQuadToTriTop && !is_toroidal){
+ if(!MeshQuadToTriTopSurface(from, to, pos))
Msg::Error("In MeshExtrudedSurface()::copyMesh(), mesh of QuadToTri top "
- "surface %d failed.", to->tag() );
+ "surface %d failed.", to->tag() );
return;
}
-
// create quadrangle elements if NOT QuadToTri and NOT toroidal
for(unsigned int i = 0; i < from->quadrangles.size(); i++){
std::vector<MVertex*> verts;
for(int j = 0; j < 4; j++){
MVertex *v = from->quadrangles[i]->getVertex(j);
- MVertex tmp(v->x(), v->y(), v->z(), 0, -1);
+ double x = v->x(), y = v->y(), z = v->z();
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
- tmp.x(), tmp.y(), tmp.z());
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end()) {
+ x, y, z);
+ MVertex *tmp = pos.find(x, y, z);
+ if(!tmp) {
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
- tmp.x(), tmp.y(), tmp.z(), to->tag());
+ x, y, z, to->tag());
return;
}
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
addQuadrangle(verts[0], verts[1], verts[2], verts[3], to);
}
@@ -265,26 +248,21 @@ int MeshExtrudedSurface(GFace *gf,
Msg::Info("Meshing surface %d (extruded)", gf->tag());
- // build a set with all the vertices on the boundary of the face gf
- double old_tol = MVertexLessThanLexicographic::tolerance;
- MVertexLessThanLexicographic::tolerance = 1.e-12 * CTX::instance()->lc;
-
- std::set<MVertex*, MVertexLessThanLexicographic> pos;
+ // build an rtree with all the vertices on the boundary of the face gf
+ MVertexRTree pos(1.e-12 * CTX::instance()->lc);
std::list<GEdge*> edges = gf->edges();
std::list<GEdge*>::iterator it = edges.begin();
while(it != edges.end()){
- pos.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
- pos.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
- (*it)->getBeginVertex()->mesh_vertices.end());
- pos.insert((*it)->getEndVertex()->mesh_vertices.begin(),
- (*it)->getEndVertex()->mesh_vertices.end());
+ pos.insert((*it)->mesh_vertices);
+ pos.insert((*it)->getBeginVertex()->mesh_vertices);
+ pos.insert((*it)->getEndVertex()->mesh_vertices);
++it;
}
// if the edges of the mesh are constrained, the vertices already
// exist on the face--so we add them to the set
if(constrainedEdges)
- pos.insert(gf->mesh_vertices.begin(), gf->mesh_vertices.end());
+ pos.insert(gf->mesh_vertices);
if(ep->geo.Mode == EXTRUDED_ENTITY) {
// surface is extruded from a curve
@@ -311,8 +289,6 @@ int MeshExtrudedSurface(GFace *gf,
copyMesh(from, gf, pos);
}
- MVertexLessThanLexicographic::tolerance = old_tol;
-
gf->meshStatistics.status = GFace::DONE;
return 1;
}
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index f4df794766b768b811fcf5b1b8a62c8599a1e137..e5317f1adb372037a23c1b132e1396833b907923 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -107,10 +107,8 @@ static void createTet(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, GRegio
}
static int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
- std::vector<MVertex*> &verts)
+ MVertexRTree &pos, std::vector<MVertex*> &verts)
{
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
double x[8], y[8], z[8];
int n = ele->getNumVertices();
for(int p = 0; p < n; p++){
@@ -124,23 +122,17 @@ static int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
ep->Extrude(j, k + 1, x[p + n], y[p + n], z[p + n]);
}
for(int p = 0; p < 2 * n; p++){
- MVertex tmp(x[p], y[p], z[p], 0, -1);
- itp = pos.find(&tmp);
- if(itp == pos.end()){ // FIXME: workaround
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
- itp = tmp.linearSearch(pos);
- }
- if(itp == pos.end())
+ MVertex *tmp = pos.find(x[p], y[p], z[p]);
+ if(!tmp)
Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g)",
- tmp.x(), tmp.y(), tmp.z());
+ x[p], y[p], z[p]);
else
- verts.push_back(*itp);
+ verts.push_back(tmp);
}
return verts.size();
}
-static void extrudeMesh(GFace *from, GRegion *to,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+static void extrudeMesh(GFace *from, GRegion *to, MVertexRTree &pos)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
@@ -165,9 +157,9 @@ static void extrudeMesh(GFace *from, GRegion *to,
return;
}
- // create elements (note that it would be faster to access the
- // *interior* nodes by direct indexing, but it's just simpler to
- // query everything by position)
+ // create elements (note that it would be faster to access the *interior*
+ // nodes by direct indexing, but it's just simpler to query everything by
+ // position)
for(unsigned int i = 0; i < from->triangles.size(); i++){
for(int j = 0; j < ep->mesh.NbLayer; j++) {
for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
@@ -195,22 +187,19 @@ static void extrudeMesh(GFace *from, GRegion *to,
}
}
-static void insertAllVertices(GRegion *gr,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+static void insertAllVertices(GRegion *gr, MVertexRTree &pos)
{
- pos.insert(gr->mesh_vertices.begin(), gr->mesh_vertices.end());
+ pos.insert(gr->mesh_vertices);
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator itf = faces.begin();
while(itf != faces.end()){
- pos.insert((*itf)->mesh_vertices.begin(), (*itf)->mesh_vertices.end());
+ pos.insert((*itf)->mesh_vertices);
std::list<GEdge*> edges = (*itf)->edges();
std::list<GEdge*>::iterator ite = edges.begin();
while(ite != edges.end()){
- pos.insert((*ite)->mesh_vertices.begin(), (*ite)->mesh_vertices.end());
- pos.insert((*ite)->getBeginVertex()->mesh_vertices.begin(),
- (*ite)->getBeginVertex()->mesh_vertices.end());
- pos.insert((*ite)->getEndVertex()->mesh_vertices.begin(),
- (*ite)->getEndVertex()->mesh_vertices.end());
+ pos.insert((*ite)->mesh_vertices);
+ pos.insert((*ite)->getBeginVertex()->mesh_vertices);
+ pos.insert((*ite)->getEndVertex()->mesh_vertices);
++ite;
}
++itf;
@@ -233,10 +222,8 @@ void meshGRegionExtruded::operator() (GRegion *gr)
deMeshGRegion dem;
dem(gr);
- // build a set with all the vertices on the boundary of gr
- double old_tol = MVertexLessThanLexicographic::tolerance;
- MVertexLessThanLexicographic::tolerance = 1.e-12 * CTX::instance()->lc;
- std::set<MVertex*, MVertexLessThanLexicographic> pos;
+ // build an rtree with all the vertices on the boundary of gr
+ MVertexRTree pos(1.e-12 * CTX::instance()->lc);
insertAllVertices(gr, pos);
// volume is extruded from a surface
@@ -255,8 +242,6 @@ void meshGRegionExtruded::operator() (GRegion *gr)
for(it = ep->mesh.Holes.begin(); it != ep->mesh.Holes.end(); it++)
carveHole(gr, it->first, it->second.first, it->second.second);
}
-
- MVertexLessThanLexicographic::tolerance = old_tol;
}
static int edgeExists(MVertex *v1, MVertex *v2,
@@ -281,8 +266,7 @@ static void deleteEdge(MVertex *v1, MVertex *v2,
}
// subdivide the 3 lateral faces of each prism
-static void phase1(GRegion *gr,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+static void phase1(GRegion *gr, MVertexRTree &pos,
std::set<std::pair<MVertex*, MVertex*> > &edges)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -316,8 +300,7 @@ static void phase1(GRegion *gr,
}
// modify lateral edges to make them "tet-compatible"
-static void phase2(GRegion *gr,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+static void phase2(GRegion *gr, MVertexRTree &pos,
std::set<std::pair<MVertex*, MVertex*> > &edges,
std::set<std::pair<MVertex*, MVertex*> > &edges_swap,
int &swap)
@@ -384,8 +367,7 @@ static void phase2(GRegion *gr,
}
// create tets
-static void phase3(GRegion *gr,
- std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+static void phase3(GRegion *gr, MVertexRTree &pos,
std::set<std::pair<MVertex*, MVertex*> > &edges)
{
ExtrudeParams *ep = gr->meshAttributes.extrude;
@@ -452,9 +434,7 @@ int SubdivideExtrudedMesh(GModel *m)
// create a vector of quadToTri regions that have NOT been meshed
// yet
std::vector<GRegion*> regions, regions_quadToTri;
- double old_tol = MVertexLessThanLexicographic::tolerance;
- MVertexLessThanLexicographic::tolerance = 1.e-12 * CTX::instance()->lc;
- std::set<MVertex*, MVertexLessThanLexicographic> pos;
+ MVertexRTree pos(1.e-12 * CTX::instance()->lc);
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
ExtrudeParams *ep = (*it)->meshAttributes.extrude;
if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY &&
@@ -539,7 +519,7 @@ int SubdivideExtrudedMesh(GModel *m)
// mesh the region (should already be done in ExtrudeMesh).
for(unsigned int i = 0; i < regions_quadToTri.size(); i++){
GRegion *gr = regions_quadToTri[i];
- std::set<MVertex*, MVertexLessThanLexicographic> pos_local;
+ MVertexRTree pos_local(1.e-12 * CTX::instance()->lc);
insertAllVertices(gr, pos_local);
meshQuadToTriRegionAfterGlobalSubdivide(gr, &edges, pos_local);
}
@@ -565,6 +545,5 @@ int SubdivideExtrudedMesh(GModel *m)
}
}
- MVertexLessThanLexicographic::tolerance = old_tol;
return 1;
}
diff --git a/Post/PViewDataListIO.cpp b/Post/PViewDataListIO.cpp
index 3228133e4c966436e71f273f56bef202828d75ee..6292094e40cf8f0e62232371b6bbf029dda955d9 100644
--- a/Post/PViewDataListIO.cpp
+++ b/Post/PViewDataListIO.cpp
@@ -11,7 +11,7 @@
#include "StringUtils.h"
#include "GmshMessage.h"
#include "GmshDefines.h"
-#include "MVertexPositionSet.h"
+#include "MVertexRTree.h"
#include "Context.h"
#include "adaptiveData.h"
#include "OS.h"
@@ -463,9 +463,8 @@ public:
};
static void createElements(std::vector<double> &list, int nbelm, int nbnod,
- MVertexPositionSet &pos, std::vector<MElement*> &elements,
- double eps, int type,
- std::map<MVertex*, nodeData> *vertexData)
+ MVertexRTree &pos, std::vector<MElement*> &elements,
+ int type, std::map<MVertex*, nodeData> *vertexData)
{
if(!nbelm) return;
int t = 0;
@@ -533,7 +532,7 @@ static void createElements(std::vector<double> &list, int nbelm, int nbnod,
double *z = &list[i + 2 * nbnod];
std::vector<MVertex*> verts(nbnod);
for(int j = 0; j < nbnod; j++){
- verts[j] = pos.find(x[j], y[j], z[j], eps);
+ verts[j] = pos.find(x[j], y[j], z[j]);
if(vertexData)
(*vertexData)[verts[j]] = nodeData(nbnod, j, &list[i + 3 * nbnod]);
}
@@ -572,7 +571,8 @@ bool PViewDataList::writeMSH(const std::string &fileName, double version, bool b
if(*numEle) numComponents = std::min(numComponents, numComp);
createVertices(*list, *numEle, numNodes, vertices);
}
- MVertexPositionSet pos(vertices);
+ MVertexRTree pos(eps);
+ pos.insert(vertices);
std::map<MVertex *, nodeData> vertexData;
@@ -580,7 +580,7 @@ bool PViewDataList::writeMSH(const std::string &fileName, double version, bool b
std::vector<double> *list = 0;
int *numEle = 0, numComp, numNodes;
int typ = _getRawData(i, &list, &numEle, &numComp, &numNodes);
- createElements(*list, *numEle, numNodes, pos, elements, eps, typ,
+ createElements(*list, *numEle, numNodes, pos, elements, typ,
forceNodeData ? &vertexData : 0);
}
diff --git a/benchmarks/2d/world.geo b/benchmarks/2d/world.geo
index a7d51c1dcdd01c9bfdaf0168cf808a5b7767675f..f7f16ff25936553343e15483235c6865a0a6f6a7 100644
--- a/benchmarks/2d/world.geo
+++ b/benchmarks/2d/world.geo
@@ -2,7 +2,6 @@ lc = 0.1;
Point(1) = {0,0,0,lc};
Point(2) = {1,0,0,lc};
-// Fabien - enleve ceci:
PolarSphere (1) = {1,2};
Point(3) = {0.335046,-1.96526,0,lc};