From 8c33ad3eedb6b987ea2a6d116d65ad67ed6c7a47 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 29 Aug 2009 11:04:08 +0000
Subject: [PATCH] never ever change the immutable vertex numbers: otherwise we
cannot use our algorithms from outside with prescribed boundary meshes
-> use MVertex::_index instead
(renamed setNum to forceNum to make sure we don't use setNum by mistake)
---
Geo/GModelIO_OCC.cpp | 9 +-
Geo/MVertex.cpp | 2 +-
Geo/MVertex.h | 18 +--
Geo/discreteEdge.cpp | 4 +-
Mesh/BDS.cpp | 12 +-
Mesh/Generator.cpp | 13 --
Mesh/meshGFaceDelaunayInsertion.cpp | 188 +++++++++++++-------------
Mesh/meshGFaceOptimize.cpp | 159 +++++++++++-----------
Mesh/meshGRegion.cpp | 67 ++++-----
Mesh/meshGRegionDelaunayInsertion.cpp | 12 +-
Mesh/meshGRegionDelaunayInsertion.h | 16 +--
Mesh/qualityMeasures.cpp | 36 ++---
12 files changed, 258 insertions(+), 278 deletions(-)
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 53793f621b..066171a1e1 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1015,7 +1015,7 @@ static void _applyOCCMeshConstraintsOnVertices
TopoDS_Shape *shape2 = (TopoDS_Shape*)gf->getNativePtr();
if(shape.IsSame(*shape2)){
Msg::Debug("... embedding vertex %d in face %d", nodeNum, gf->tag());
- gv->mesh_vertices[0]->setNum(nodeNum);
+ gv->mesh_vertices[0]->forceNum(nodeNum);
gf->addEmbeddedVertex(gv);
}
}
@@ -1057,8 +1057,8 @@ static void _applyOCCMeshConstraintsOnEdges
int numbeg = nodeNum.Value(invert ? n : 1);
int numend = nodeNum.Value(invert ? 1 : n);
Msg::Debug("... beg=%d end=%d", numbeg, numend);
- ge->getBeginVertex()->mesh_vertices[0]->setNum(numbeg);
- ge->getEndVertex()->mesh_vertices[0]->setNum(numend);
+ ge->getBeginVertex()->mesh_vertices[0]->forceNum(numbeg);
+ ge->getEndVertex()->mesh_vertices[0]->forceNum(numend);
// set the mesh on the edge
for(int i = 2; i < n; i++){
int num = nodeNum.Value(invert ? n - i + 1 : i);
@@ -1066,8 +1066,7 @@ static void _applyOCCMeshConstraintsOnEdges
GPoint p = ge->point(u);
Msg::Debug("... adding mesh vertex num=%d u=%g xyz=(%g,%g,%g)",
num, u, p.x(), p.y(), p.z());
- MEdgeVertex *v = new MEdgeVertex(p.x(), p.y(), p.z(), ge, u);
- v->setNum(num);
+ MEdgeVertex *v = new MEdgeVertex(p.x(), p.y(), p.z(), ge, u, -1., num);
ge->mesh_vertices.push_back(v);
}
for(unsigned int i = 0; i < ge->mesh_vertices.size() + 1; i++){
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 7c7c5dab2e..3f9b6e9759 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -42,7 +42,7 @@ MVertex::MVertex(double x, double y, double z, GEntity *ge, int num)
}
}
-void MVertex::setNum(int num)
+void MVertex::forceNum(int num)
{
#pragma omp critical
{
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 1cf422768f..78cb939917 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -27,13 +27,13 @@ class MVertex{
protected:
// the maximum vertex id number in the mesh
static int _globalNum;
- // the id number of the vertex (this number is unique and is
- // guaranteed never to change once a mesh has been generated)
+ // the immutable id number of the vertex (this number is unique and
+ // is guaranteed never to change once a vertex has been created)
int _num;
- // a vertex index, used for example when saving a mesh (this index
- // is not necessarily unique, can change after mesh renumbering,
- // ...). By convention, vertices with negative indices are not
- // saved
+ // a vertex index, used for example during mesh generation or when
+ // saving a mesh (this index is not necessarily unique, can change
+ // after mesh renumbering, etc.). By convention, vertices with
+ // negative indices are not saved
int _index;
// a visibility and polynomial order flags
char _visible, _order;
@@ -71,9 +71,11 @@ class MVertex{
inline GEntity* onWhat() const { return _ge; }
inline void setEntity(GEntity *ge) { _ge = ge; }
- // get/set the number
+ // get the immutable vertex number
inline int getNum() const { return _num; }
- void setNum(int num);
+
+ // force the immutable number (this should normally never be used)
+ void forceNum(int num);
// get/set the index
inline int getIndex() const { return _index; }
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index c995316213..4caeca1900 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -218,10 +218,10 @@ void discreteEdge::parametrize()
if (_orientation[i] == 1 ) vR = lines[i]->getVertex(1);
else vR = lines[i]->getVertex(0);
int param = i+1;
- MVertex *vNEW = new MEdgeVertex(vR->x(),vR->y(),vR->z(), this, param);
+ MVertex *vNEW = new MEdgeVertex(vR->x(),vR->y(),vR->z(), this, param,
+ -1., vR->getNum());
old2new.insert(std::make_pair(vR,vNEW));
newVertices.push_back(vNEW);
- vNEW->setNum(vR->getNum());
newLines.push_back(new MLine(vL, vNEW));
_orientation[i] = 1;
vL = vNEW;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 0731fea655..a341dcd827 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -170,7 +170,7 @@ BDS_Point *BDS_Mesh::add_point(int num, double u, double v, GFace *gf)
BDS_Point *BDS_Mesh::find_point(int p)
{
BDS_Point P(p);
- std::set < BDS_Point *, PointLessThan >::iterator it = points.find(&P);
+ std::set<BDS_Point *, PointLessThan>::iterator it = points.find(&P);
if(it != points.end())
return (BDS_Point*)(*it);
else
@@ -331,7 +331,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
outputScalarField(triangles, "debugp.pos", 1);
outputScalarField(triangles, "debugr.pos", 0);
Msg::Debug("edge %d %d cannot be recovered at all, look at debugp.pos "
- "and debugr.pos", num1, num2);
+ "and debugr.pos", num1, num2);
return 0;
}
return eee;
@@ -503,7 +503,7 @@ void BDS_Edge::oppositeof(BDS_Point *oface[2]) const
BDS_GeomEntity *BDS_Mesh::get_geom(int p1, int p2)
{
BDS_GeomEntity ge(p1, p2);
- std::set < BDS_GeomEntity *, GeomLessThan >::iterator it = geom.find(&ge);
+ std::set<BDS_GeomEntity *, GeomLessThan >::iterator it = geom.find(&ge);
if(it == geom.end()) return 0;
return (BDS_GeomEntity*)(*it);
}
@@ -528,7 +528,7 @@ void recur_tag(BDS_Face *t, BDS_GeomEntity *g)
double PointLessThanLexicographic::t = 0;
double BDS_Vector::t = 0;
-template < class IT > void DESTROOOY(IT beg, IT end)
+template <class IT> void DESTROOOY(IT beg, IT end)
{
while(beg != end) {
delete *beg;
@@ -1177,8 +1177,8 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
std::list<BDS_Face*> ts;
p->getTriangles(ts);
- std::list < BDS_Edge * >::iterator ited = p->edges.begin();
- std::list < BDS_Edge * >::iterator itede = p->edges.end();
+ std::list<BDS_Edge *>::iterator ited = p->edges.begin();
+ std::list<BDS_Edge *>::iterator itede = p->edges.end();
double sTot = 0;
while(ited != itede) {
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 67d4409472..1782257c86 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -409,20 +409,7 @@ static void Mesh2D(GModel *m)
// and curve meshes) is global as it depends on a smooth normal
// field generated from the surface mesh of the source surfaces
if(!Mesh2DWithBoundaryLayers(m)){
-
-#if 1
std::for_each(m->firstFace(), m->lastFace(), meshGFace());
-#else
- // test openmp
- std::vector<GFace*> faces;
- faces.insert(faces.begin(), m->firstFace(), m->lastFace());
-#pragma omp parallel for schedule(dynamic)
- for(unsigned int i = 0; i < faces.size(); i++){
- meshGFace mesher;
- mesher(faces[i]);
- }
-#endif
-
int nIter = 0;
while(1){
meshGFace mesher;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 28ea0a71f0..6379ccd2dc 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -117,12 +117,12 @@ void circumCenterMetric(MTriangle *base, const double *metric,
double *x, double &Radius2)
{
// d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
- double pa[2] = {Us[base->getVertex(0)->getNum()],
- Vs[base->getVertex(0)->getNum()]};
- double pb[2] = {Us[base->getVertex(1)->getNum()],
- Vs[base->getVertex(1)->getNum()]};
- double pc[2] = {Us[base->getVertex(2)->getNum()],
- Vs[base->getVertex(2)->getNum()]};
+ double pa[2] = {Us[base->getVertex(0)->getIndex()],
+ Vs[base->getVertex(0)->getIndex()]};
+ double pb[2] = {Us[base->getVertex(1)->getIndex()],
+ Vs[base->getVertex(1)->getIndex()]};
+ double pc[2] = {Us[base->getVertex(2)->getIndex()],
+ Vs[base->getVertex(2)->getIndex()]};
circumCenterMetric(pa, pb, pc, metric, x, Radius2);
}
@@ -186,10 +186,12 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base,
const std::vector<double> &Vs)
{
double x[2], Radius2, metric[3];
- double pa[2] = {(Us[base->getVertex(0)->getNum()] + Us[base->getVertex(1)->getNum()] +
- Us[base->getVertex(2)->getNum()]) / 3.,
- (Vs[base->getVertex(0)->getNum()] + Vs[base->getVertex(1)->getNum()] +
- Vs[base->getVertex(2)->getNum()]) / 3.};
+ double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
+ Us[base->getVertex(1)->getIndex()] +
+ Us[base->getVertex(2)->getIndex()]) / 3.,
+ (Vs[base->getVertex(0)->getIndex()] +
+ Vs[base->getVertex(1)->getIndex()] +
+ Vs[base->getVertex(2)->getIndex()]) / 3.};
buildMetric(gf, pa, metric);
circumCenterMetric(base, metric, Us, Vs, x, Radius2);
@@ -253,12 +255,12 @@ int inCircumCircle(MTriangle *base,
const double *p, const double *param ,
std::vector<double> &Us, std::vector<double> &Vs)
{
- double pa[2] = {Us[base->getVertex(0)->getNum()],
- Vs[base->getVertex(0)->getNum()]};
- double pb[2] = {Us[base->getVertex(1)->getNum()],
- Vs[base->getVertex(1)->getNum()]};
- double pc[2] = {Us[base->getVertex(2)->getNum()],
- Vs[base->getVertex(2)->getNum()]};
+ double pa[2] = {Us[base->getVertex(0)->getIndex()],
+ Vs[base->getVertex(0)->getIndex()]};
+ double pb[2] = {Us[base->getVertex(1)->getIndex()],
+ Vs[base->getVertex(1)->getIndex()]};
+ double pc[2] = {Us[base->getVertex(2)->getIndex()],
+ Vs[base->getVertex(2)->getIndex()]};
double result =
gmsh::incircle(pa, pb, pc, (double*)param) * gmsh::orient2d(pa, pb, pc);
@@ -351,14 +353,14 @@ void recurFindCavityAniso(GFace *gf,
bool circUV(MTriangle *t, std::vector<double> & Us, std::vector<double> &Vs,
double *res, GFace *gf)
{
- double u1 [3] = {Us[t->getVertex(0)->getNum()], Vs[t->getVertex(0)->getNum()], 0};
- double u2 [3] = {Us[t->getVertex(1)->getNum()], Vs[t->getVertex(1)->getNum()], 0};
- double u3 [3] = {Us[t->getVertex(2)->getNum()], Vs[t->getVertex(2)->getNum()], 0};
+ double u1[3] = {Us[t->getVertex(0)->getIndex()], Vs[t->getVertex(0)->getIndex()], 0};
+ double u2[3] = {Us[t->getVertex(1)->getIndex()], Vs[t->getVertex(1)->getIndex()], 0};
+ double u3[3] = {Us[t->getVertex(2)->getIndex()], Vs[t->getVertex(2)->getIndex()], 0};
circumCenterXY(u1, u2, u3, res);
return true;
- double p1 [3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
- double p2 [3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
- double p3 [3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
+ double p1[3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
+ double p2[3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
+ double p3[3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
double resxy[3], uv[2];
circumCenterXYZ(p1, p2, p3, resxy,uv);
return true;
@@ -371,12 +373,12 @@ bool invMapUV(MTriangle *t, double *p,
double mat[2][2];
double b[2];
- double u0 = Us[t->getVertex(0)->getNum()];
- double v0 = Vs[t->getVertex(0)->getNum()];
- double u1 = Us[t->getVertex(1)->getNum()];
- double v1 = Vs[t->getVertex(1)->getNum()];
- double u2 = Us[t->getVertex(2)->getNum()];
- double v2 = Vs[t->getVertex(2)->getNum()];
+ double u0 = Us[t->getVertex(0)->getIndex()];
+ double v0 = Vs[t->getVertex(0)->getIndex()];
+ double u1 = Us[t->getVertex(1)->getIndex()];
+ double v1 = Vs[t->getVertex(1)->getIndex()];
+ double u2 = Us[t->getVertex(2)->getIndex()];
+ double v2 = Vs[t->getVertex(2)->getIndex()];
mat[0][0] = u1 - u0;
mat[0][1] = u2 - u0;
@@ -399,14 +401,14 @@ bool invMapUV(MTriangle *t, double *p,
double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
{
- double u1 = Us[t->getVertex(0)->getNum()];
- double v1 = Vs[t->getVertex(0)->getNum()];
- double u2 = Us[t->getVertex(1)->getNum()];
- double v2 = Vs[t->getVertex(1)->getNum()];
- double u3 = Us[t->getVertex(2)->getNum()];
- double v3 = Vs[t->getVertex(2)->getNum()];
- const double vv1 [2] = {u2 - u1, v2 - v1};
- const double vv2 [2] = {u3 - u1, v3 - v1};
+ double u1 = Us[t->getVertex(0)->getIndex()];
+ double v1 = Vs[t->getVertex(0)->getIndex()];
+ double u2 = Us[t->getVertex(1)->getIndex()];
+ double v2 = Vs[t->getVertex(1)->getIndex()];
+ double u3 = Us[t->getVertex(2)->getIndex()];
+ double v3 = Vs[t->getVertex(2)->getIndex()];
+ const double vv1[2] = {u2 - u1, v2 - v1};
+ const double vv2[2] = {u3 - u1, v3 - v1};
double s = vv1[0] * vv2[1] - vv1[1] * vv2[0];
return s * 0.5;
}
@@ -452,20 +454,20 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
while (it != shell.end()){
MTriangle *t = new MTriangle(it->v[0], it->v[1], v);
const double ONE_THIRD = 1./3.;
- double lc = ONE_THIRD * (vSizes[t->getVertex(0)->getNum()] +
- vSizes[t->getVertex(1)->getNum()] +
- vSizes[t->getVertex(2)->getNum()]);
- double lcBGM = ONE_THIRD * (vSizesBGM[t->getVertex(0)->getNum()] +
- vSizesBGM[t->getVertex(1)->getNum()] +
- vSizesBGM[t->getVertex(2)->getNum()]);
+ double lc = ONE_THIRD * (vSizes[t->getVertex(0)->getIndex()] +
+ vSizes[t->getVertex(1)->getIndex()] +
+ vSizes[t->getVertex(2)->getIndex()]);
+ double lcBGM = ONE_THIRD * (vSizesBGM[t->getVertex(0)->getIndex()] +
+ vSizesBGM[t->getVertex(1)->getIndex()] +
+ vSizesBGM[t->getVertex(2)->getIndex()]);
double LL = Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM;
MTri3 *t4;
if (_experimental_anisotropic_blues_band_){
- SMetric3 metrBGM = interpolation ( vMetricsBGM[t->getVertex(0)->getNum()],
- vMetricsBGM[t->getVertex(1)->getNum()],
- vMetricsBGM[t->getVertex(2)->getNum()],
- ONE_THIRD,ONE_THIRD);
+ SMetric3 metrBGM = interpolation(vMetricsBGM[t->getVertex(0)->getIndex()],
+ vMetricsBGM[t->getVertex(1)->getIndex()],
+ vMetricsBGM[t->getVertex(2)->getIndex()],
+ ONE_THIRD, ONE_THIRD);
t4 = new MTri3(t, LL,&metrBGM);
}
@@ -537,14 +539,14 @@ void _printTris(char *name, std::set<MTri3*, compareTri3Ptr> &AllTris,
if (!worst->isDeleted()){
if (param)
fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- Us [(worst)->tri()->getVertex(0)->getNum()],
- Vs [(worst)->tri()->getVertex(0)->getNum()],
+ Us[(worst)->tri()->getVertex(0)->getIndex()],
+ Vs[(worst)->tri()->getVertex(0)->getIndex()],
0.0,
- Us [(worst)->tri()->getVertex(1)->getNum()],
- Vs [(worst)->tri()->getVertex(1)->getNum()],
+ Us[(worst)->tri()->getVertex(1)->getIndex()],
+ Vs[(worst)->tri()->getVertex(1)->getIndex()],
0.0,
- Us [(worst)->tri()->getVertex(2)->getNum()],
- Vs [(worst)->tri()->getVertex(2)->getNum()],
+ Us[(worst)->tri()->getVertex(2)->getIndex()],
+ Vs[(worst)->tri()->getVertex(2)->getIndex()],
0.0,
(worst)->getRadius(),
(worst)->getRadius(),
@@ -616,10 +618,10 @@ static void insertAPoint(GFace *gf,
GPoint p = gf->point(center[0], center[1]);
// printf("the new point is %g %g\n",p.x(),p.y());
MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]);
- v->setNum(Us.size());
- double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getNum()] +
- uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] +
- uv[1] * vSizes [ptin->tri()->getVertex(2)->getNum()]);
+ v->setIndex(Us.size());
+ double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getIndex()] +
+ uv[0] * vSizes[ptin->tri()->getVertex(1)->getIndex()] +
+ uv[1] * vSizes[ptin->tri()->getVertex(2)->getIndex()]);
double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
//SMetric3 metr = BGM_MeshMetric(gf,center[0],center[1],p.x(),p.y(),p.z());
// vMetricsBGM.push_back(metr);
@@ -642,14 +644,14 @@ static void insertAPoint(GFace *gf,
}
}
else {
- /* Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",
+ /* Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",
center[0], center[1],
- Us[base->getVertex(0)->getNum()],
- Vs[base->getVertex(0)->getNum()],
- Us[base->getVertex(1)->getNum()],
- Vs[base->getVertex(1)->getNum()],
- Us[base->getVertex(2)->getNum()],
- Vs[base->getVertex(2)->getNum()],
+ Us[base->getVertex(0)->getIndex()],
+ Vs[base->getVertex(0)->getIndex()],
+ Us[base->getVertex(1)->getIndex()],
+ Vs[base->getVertex(1)->getIndex()],
+ Us[base->getVertex(2)->getIndex()],
+ Vs[base->getVertex(2)->getIndex()],
metric[0], metric[1], metric[2]);*/
AllTris.erase(it);
worst->forceRadius(0);
@@ -682,21 +684,21 @@ void gmshBowyerWatson(GFace *gf)
else{
if(ITER++ % 5000 == 0)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
- vSizes.size(), worst->getRadius());
+ vSizes.size(), worst->getRadius());
double center[2],metric[3],r2;
if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
circUV(worst->tri(), Us, Vs, center, gf);
MTriangle *base = worst->tri();
- double pa[2] = {(Us[base->getVertex(0)->getNum()] +
- Us[base->getVertex(1)->getNum()] +
- Us[base->getVertex(2)->getNum()]) / 3.,
- (Vs[base->getVertex(0)->getNum()] +
- Vs[base->getVertex(1)->getNum()] +
- Vs[base->getVertex(2)->getNum()]) / 3.};
+ double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
+ Us[base->getVertex(1)->getIndex()] +
+ Us[base->getVertex(2)->getIndex()]) / 3.,
+ (Vs[base->getVertex(0)->getIndex()] +
+ Vs[base->getVertex(1)->getIndex()] +
+ Vs[base->getVertex(2)->getIndex()]) / 3.};
if (_experimental_anisotropic_blues_band_){
- SMetric3 m = interpolation (vMetricsBGM[base->getVertex(0)->getNum()],
- vMetricsBGM[base->getVertex(1)->getNum()],
- vMetricsBGM[base->getVertex(2)->getNum()],
+ SMetric3 m = interpolation (vMetricsBGM[base->getVertex(0)->getIndex()],
+ vMetricsBGM[base->getVertex(1)->getIndex()],
+ vMetricsBGM[base->getVertex(2)->getIndex()],
pa[0],pa[1]);
buildMetric(gf, pa, m, metric);
}
@@ -726,9 +728,9 @@ void gmshBowyerWatson(GFace *gf)
static double lengthMetric(const double p[2], const double q[2],
const double metric[3])
{
- return sqrt ( (p[0]-q[0]) * metric [0] * (p[0]-q[0]) +
- 2 * (p[0]-q[0]) * metric [1] * (p[1]-q[1]) +
- (p[1]-q[1]) * metric [2] * (p[1]-q[1]) );
+ return sqrt ( (p[0] - q[0]) * metric[0] * (p[0] - q[0]) +
+ 2 * (p[0] - q[0]) * metric[1] * (p[1] - q[1]) +
+ (p[1] - q[1]) * metric[2] * (p[1] - q[1]) );
}
/*
@@ -807,17 +809,17 @@ void gmshBowyerWatsonFrontal(GFace *gf)
worst->getRadius() > LIMIT_){
if(ITER++ % 5000 == 0)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
- vSizes.size(), worst->getRadius());
+ vSizes.size(), worst->getRadius());
// compute circum center of that guy
double center[2],metric[3],r2;
MTriangle *base = worst->tri();
circUV(base, Us, Vs, center, gf);
- double pa[2] = {(Us[base->getVertex(0)->getNum()] +
- Us[base->getVertex(1)->getNum()] +
- Us[base->getVertex(2)->getNum()]) / 3.,
- (Vs[base->getVertex(0)->getNum()] +
- Vs[base->getVertex(1)->getNum()] +
- Vs[base->getVertex(2)->getNum()]) / 3.};
+ double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
+ Us[base->getVertex(1)->getIndex()] +
+ Us[base->getVertex(2)->getIndex()]) / 3.,
+ (Vs[base->getVertex(0)->getIndex()] +
+ Vs[base->getVertex(1)->getIndex()] +
+ Vs[base->getVertex(2)->getIndex()]) / 3.};
buildMetric(gf, pa, metric);
circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
// compute the middle point of the edge
@@ -826,10 +828,10 @@ void gmshBowyerWatsonFrontal(GFace *gf)
// printf("the active edge is %d : %g %g -> %g %g\n",
// active_edge,base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),
// base->getVertex(ip2)->x(),base->getVertex(ip2)->y());
- double P[2] = {Us[base->getVertex(ip1)->getNum()],
- Vs[base->getVertex(ip1)->getNum()]};
- double Q[2] = {Us[base->getVertex(ip2)->getNum()],
- Vs[base->getVertex(ip2)->getNum()]};
+ double P[2] = {Us[base->getVertex(ip1)->getIndex()],
+ Vs[base->getVertex(ip1)->getIndex()]};
+ double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
+ Vs[base->getVertex(ip2)->getIndex()]};
double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
// now we have the edge center and the center of the circumcircle,
@@ -847,15 +849,15 @@ void gmshBowyerWatsonFrontal(GFace *gf)
const double p = 0.5 * lengthMetric(P, Q, metric); // / RATIO;
const double q = lengthMetric(center, midpoint, metric);
const double rhoM1 = 0.5 *
- (vSizes[base->getVertex(ip1)->getNum()] +
- vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
+ (vSizes[base->getVertex(ip1)->getIndex()] +
+ vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
const double rhoM2 = 0.5 *
- (vSizesBGM[base->getVertex(ip1)->getNum()] +
- vSizesBGM[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ vSizesBGM[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
// const double rhoM = 0.5 *
- // (vSizes[base->getVertex(ip1)->getNum()] +
- // vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
+ // (vSizes[base->getVertex(ip1)->getIndex()] +
+ // vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
const double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) / RATIO;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 88ce8456c6..6cd2547914 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -34,9 +34,9 @@ static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes)
for (int j = i + 1; j < 3; j++){
MVertex *vi = t->getVertex(i);
MVertex *vj = t->getVertex(j);
- double dx = vi->x()-vj->x();
- double dy = vi->y()-vj->y();
- double dz = vi->z()-vj->z();
+ double dx = vi->x() - vj->x();
+ double dy = vi->y() - vj->y();
+ double dz = vi->z() - vj->z();
double l = sqrt(dx * dx + dy * dy + dz * dz);
std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
@@ -77,7 +77,10 @@ void buildMeshGenerationDataStructures(GFace *gf,
int NUM = 0;
for (std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
- it->first->setNum(NUM++);
+ // FIXME: this vertex-stored indexing makes the optimization
+ // routines not thread-safe (we cannot concurrently optimize two
+ // surfaces that share an edge)
+ it->first->setIndex(NUM++);
vSizes.push_back(it->second);
vSizesBGM.push_back(it->second);
vMetricsBGM.push_back(SMetric3(it->second));
@@ -87,9 +90,9 @@ void buildMeshGenerationDataStructures(GFace *gf,
Vs.push_back(param[1]);
}
for(unsigned int i = 0; i < gf->triangles.size(); i++){
- double lc = 0.3333333333 * (vSizes[gf->triangles[i]->getVertex(0)->getNum()] +
- vSizes[gf->triangles[i]->getVertex(1)->getNum()] +
- vSizes[gf->triangles[i]->getVertex(2)->getNum()]);
+ double lc = 0.3333333333 * (vSizes[gf->triangles[i]->getVertex(0)->getIndex()] +
+ vSizes[gf->triangles[i]->getVertex(1)->getIndex()] +
+ vSizes[gf->triangles[i]->getVertex(2)->getIndex()]);
AllTris.insert(new MTri3(gf->triangles[i], lc));
}
gf->triangles.clear();
@@ -118,15 +121,15 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
double n1[3], n2[3];
MTriangle *t = gf->triangles[0];
MVertex *v0 = t->getVertex(0), *v1 = t->getVertex(1), *v2 = t->getVertex(2);
- normal3points(Us[v0->getNum()], Vs[v0->getNum()], 0.,
- Us[v1->getNum()], Vs[v1->getNum()], 0.,
- Us[v2->getNum()], Vs[v2->getNum()], 0., n1);
+ normal3points(Us[v0->getIndex()], Vs[v0->getIndex()], 0.,
+ Us[v1->getIndex()], Vs[v1->getIndex()], 0.,
+ Us[v2->getIndex()], Vs[v2->getIndex()], 0., n1);
for(unsigned int j = 1; j < gf->triangles.size(); j++){
t = gf->triangles[j];
v0 = t->getVertex(0); v1 = t->getVertex(1); v2 = t->getVertex(2);
- normal3points(Us[v0->getNum()], Vs[v0->getNum()], 0.,
- Us[v1->getNum()], Vs[v1->getNum()], 0.,
- Us[v2->getNum()], Vs[v2->getNum()], 0., n2);
+ normal3points(Us[v0->getIndex()], Vs[v0->getIndex()], 0.,
+ Us[v1->getIndex()], Vs[v1->getIndex()], 0.,
+ Us[v2->getIndex()], Vs[v2->getIndex()], 0., n2);
double pp; prosca(n1, n2, &pp);
if(pp < 0) t->revert();
}
@@ -253,10 +256,10 @@ double surfaceTriangleUV(MVertex *v1, MVertex *v2, MVertex *v3,
const std::vector<double> &Us,
const std::vector<double> &Vs)
{
- const double v12[2] = {Us[v2->getNum()] - Us[v1->getNum()],
- Vs[v2->getNum()] - Vs[v1->getNum()]};
- const double v13[2] = {Us[v3->getNum()] - Us[v1->getNum()],
- Vs[v3->getNum()] - Vs[v1->getNum()]};
+ const double v12[2] = {Us[v2->getIndex()] - Us[v1->getIndex()],
+ Vs[v2->getIndex()] - Vs[v1->getIndex()]};
+ const double v13[2] = {Us[v3->getIndex()] - Us[v1->getIndex()],
+ Vs[v3->getIndex()] - Vs[v1->getIndex()]};
return 0.5 * fabs (v12[0] * v13[1] - v12[1] * v13[0]);
}
@@ -312,11 +315,11 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
}
case SWCR_DEL:
{
- double edgeCenter[2] ={(Us[v1->getNum()] + Us[v2->getNum()] + Us[v3->getNum()] +
- Us[v4->getNum()]) * .25,
- (Vs[v1->getNum()] + Vs[v2->getNum()] + Vs[v3->getNum()] +
- Vs[v4->getNum()]) * .25};
- double uv4[2] ={Us[v4->getNum()], Vs[v4->getNum()]};
+ double edgeCenter[2] ={(Us[v1->getIndex()] + Us[v2->getIndex()] + Us[v3->getIndex()] +
+ Us[v4->getIndex()]) * .25,
+ (Vs[v1->getIndex()] + Vs[v2->getIndex()] + Vs[v3->getIndex()] +
+ Vs[v4->getIndex()]) * .25};
+ double uv4[2] ={Us[v4->getIndex()], Vs[v4->getIndex()]};
double metric[3];
buildMetric(gf, edgeCenter, metric);
if (!inCircumCircleAniso(gf, t1->tri(), uv4, metric, Us, Vs)){
@@ -333,10 +336,10 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
double avg2[3] = {(v3->x() + v4->x()) *.5,(v3->y() + v4->y()) *.5,
(v3->z() + v4->z()) *.5};
- GPoint gp1 = gf->point(SPoint2((Us[v1->getNum()] + Us[v2->getNum()]) * .5,
- (Vs[v1->getNum()] + Vs[v2->getNum()]) * .5));
- GPoint gp2 = gf->point(SPoint2((Us[v3->getNum()] + Us[v4->getNum()]) * .5,
- (Vs[v3->getNum()] + Vs[v4->getNum()]) * .5));
+ GPoint gp1 = gf->point(SPoint2((Us[v1->getIndex()] + Us[v2->getIndex()]) * .5,
+ (Vs[v1->getIndex()] + Vs[v2->getIndex()]) * .5));
+ GPoint gp2 = gf->point(SPoint2((Us[v3->getIndex()] + Us[v4->getIndex()]) * .5,
+ (Vs[v3->getIndex()] + Vs[v4->getIndex()]) * .5));
double d1 = (avg1[0] - gp1.x()) * (avg1[0] - gp1.x()) + (avg1[1] - gp1.y()) *
(avg1[1]-gp1.y()) + (avg1[2] - gp1.z()) * (avg1[2] - gp1.z());
double d2 = (avg2[0] - gp2.x()) * (avg2[0] - gp2.x()) + (avg2[1] - gp2.y()) *
@@ -372,18 +375,18 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
if (!found)cavity.push_back(t2->getNeigh(i));
}
}
- double lc1 = 0.3333333333 * (vSizes[t1b->getVertex(0)->getNum()] +
- vSizes[t1b->getVertex(1)->getNum()] +
- vSizes[t1b->getVertex(2)->getNum()]);
- double lcBGM1 = 0.3333333333 * (vSizesBGM[t1b->getVertex(0)->getNum()] +
- vSizesBGM[t1b->getVertex(1)->getNum()] +
- vSizesBGM[t1b->getVertex(2)->getNum()]);
- double lc2 = 0.3333333333 * (vSizes[t2b->getVertex(0)->getNum()] +
- vSizes[t2b->getVertex(1)->getNum()] +
- vSizes[t2b->getVertex(2)->getNum()]);
- double lcBGM2 = 0.3333333333 * (vSizesBGM[t2b->getVertex(0)->getNum()] +
- vSizesBGM[t2b->getVertex(1)->getNum()] +
- vSizesBGM[t2b->getVertex(2)->getNum()]);
+ double lc1 = 0.3333333333 * (vSizes[t1b->getVertex(0)->getIndex()] +
+ vSizes[t1b->getVertex(1)->getIndex()] +
+ vSizes[t1b->getVertex(2)->getIndex()]);
+ double lcBGM1 = 0.3333333333 * (vSizesBGM[t1b->getVertex(0)->getIndex()] +
+ vSizesBGM[t1b->getVertex(1)->getIndex()] +
+ vSizesBGM[t1b->getVertex(2)->getIndex()]);
+ double lc2 = 0.3333333333 * (vSizes[t2b->getVertex(0)->getIndex()] +
+ vSizes[t2b->getVertex(1)->getIndex()] +
+ vSizes[t2b->getVertex(2)->getIndex()]);
+ double lcBGM2 = 0.3333333333 * (vSizesBGM[t2b->getVertex(0)->getIndex()] +
+ vSizesBGM[t2b->getVertex(1)->getIndex()] +
+ vSizesBGM[t2b->getVertex(2)->getIndex()]);
MTri3 *t1b3 = new MTri3(t1b, Extend1dMeshIn2dSurfaces() ?
std::min(lc1, lcBGM1) : lcBGM1);
MTri3 *t2b3 = new MTri3(t2b, Extend1dMeshIn2dSurfaces() ?
@@ -406,8 +409,8 @@ inline double computeEdgeAdimLength(MVertex *v1, MVertex *v2, GFace *f,
const std::vector<double> &vSizes ,
const std::vector<double> &vSizesBGM)
{
- const double edgeCenter[2] ={(Us[v1->getNum()] + Us[v2->getNum()]) * .5,
- (Vs[v1->getNum()] + Vs[v2->getNum()]) * .5};
+ const double edgeCenter[2] ={(Us[v1->getIndex()] + Us[v2->getIndex()]) * .5,
+ (Vs[v1->getIndex()] + Vs[v2->getIndex()]) * .5};
GPoint GP = f->point (edgeCenter[0], edgeCenter[1]);
const double dx1 = v1->x() - GP.x();
@@ -419,9 +422,9 @@ inline double computeEdgeAdimLength(MVertex *v1, MVertex *v2, GFace *f,
const double dz2 = v2->z() - GP.z();
const double l2 = sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
if (Extend1dMeshIn2dSurfaces())
- return 2 * (l1 + l2) / (std::min(vSizes[v1->getNum()], vSizesBGM[v1->getNum()]) +
- std::min(vSizes[v2->getNum()], vSizesBGM[v2->getNum()]));
- return 2 * (l1 + l2) / (vSizesBGM[v1->getNum()] + vSizesBGM[v2->getNum()]);
+ return 2 * (l1 + l2) / (std::min(vSizes[v1->getIndex()], vSizesBGM[v1->getIndex()]) +
+ std::min(vSizes[v2->getIndex()], vSizesBGM[v2->getIndex()]));
+ return 2 * (l1 + l2) / (vSizesBGM[v1->getIndex()] + vSizesBGM[v2->getIndex()]);
}
bool gmshEdgeSplit(const double lMax, MTri3 *t1, GFace *gf, int iLocalEdge,
@@ -434,8 +437,8 @@ bool gmshEdgeSplit(const double lMax, MTri3 *t1, GFace *gf, int iLocalEdge,
MVertex *v1 = t1->tri()->getVertex(iLocalEdge == 0 ? 2 : iLocalEdge - 1);
MVertex *v2 = t1->tri()->getVertex(iLocalEdge % 3);
- double edgeCenter[2] ={(Us[v1->getNum()] + Us[v2->getNum()]) * .5,
- (Vs[v1->getNum()] + Vs[v2->getNum()]) * .5};
+ double edgeCenter[2] ={(Us[v1->getIndex()] + Us[v2->getIndex()]) * .5,
+ (Vs[v1->getIndex()] + Vs[v2->getIndex()]) * .5};
double al = computeEdgeAdimLength(v1, v2, gf,Us, Vs, vSizes, vSizesBGM);
if (al <= lMax) return false;
@@ -483,8 +486,8 @@ bool gmshEdgeSplit(const double lMax, MTri3 *t1, GFace *gf, int iLocalEdge,
}
gf->mesh_vertices.push_back(vnew);
- vnew->setNum(Us.size());
- double lcN = 0.5 * (vSizes [v1->getNum()] + vSizes [v2->getNum()] );
+ vnew->setIndex(Us.size());
+ double lcN = 0.5 * (vSizes[v1->getIndex()] + vSizes[v2->getIndex()] );
double lcBGM = BGM_MeshSize(gf, edgeCenter[0], edgeCenter[1], p.x(), p.y(), p.z());
vSizesBGM.push_back(lcBGM);
@@ -511,30 +514,30 @@ bool gmshEdgeSplit(const double lMax, MTri3 *t1, GFace *gf, int iLocalEdge,
if (!found) cavity.push_back(t2->getNeigh(i));
}
}
- double lc1 = 0.3333333333 * (vSizes[t1b->getVertex(0)->getNum()] +
- vSizes[t1b->getVertex(1)->getNum()] +
- vSizes[t1b->getVertex(2)->getNum()]);
- double lcBGM1 = 0.3333333333 * (vSizesBGM[t1b->getVertex(0)->getNum()] +
- vSizesBGM[t1b->getVertex(1)->getNum()] +
- vSizesBGM[t1b->getVertex(2)->getNum()]);
- double lc2 = 0.3333333333 * (vSizes[t2b->getVertex(0)->getNum()] +
- vSizes[t2b->getVertex(1)->getNum()] +
- vSizes[t2b->getVertex(2)->getNum()]);
- double lcBGM2 = 0.3333333333 * (vSizesBGM[t2b->getVertex(0)->getNum()] +
- vSizesBGM[t2b->getVertex(1)->getNum()] +
- vSizesBGM[t2b->getVertex(2)->getNum()]);
- double lc3 = 0.3333333333 * (vSizes[t3b->getVertex(0)->getNum()] +
- vSizes[t3b->getVertex(1)->getNum()] +
- vSizes[t3b->getVertex(2)->getNum()]);
- double lcBGM3 = 0.3333333333 * (vSizesBGM[t3b->getVertex(0)->getNum()] +
- vSizesBGM[t3b->getVertex(1)->getNum()] +
- vSizesBGM[t3b->getVertex(2)->getNum()]);
- double lc4 = 0.3333333333 * (vSizes[t4b->getVertex(0)->getNum()] +
- vSizes[t4b->getVertex(1)->getNum()] +
- vSizes[t4b->getVertex(2)->getNum()]);
- double lcBGM4 = 0.3333333333 * (vSizesBGM[t4b->getVertex(0)->getNum()] +
- vSizesBGM[t4b->getVertex(1)->getNum()] +
- vSizesBGM[t4b->getVertex(2)->getNum()]);
+ double lc1 = 0.3333333333 * (vSizes[t1b->getVertex(0)->getIndex()] +
+ vSizes[t1b->getVertex(1)->getIndex()] +
+ vSizes[t1b->getVertex(2)->getIndex()]);
+ double lcBGM1 = 0.3333333333 * (vSizesBGM[t1b->getVertex(0)->getIndex()] +
+ vSizesBGM[t1b->getVertex(1)->getIndex()] +
+ vSizesBGM[t1b->getVertex(2)->getIndex()]);
+ double lc2 = 0.3333333333 * (vSizes[t2b->getVertex(0)->getIndex()] +
+ vSizes[t2b->getVertex(1)->getIndex()] +
+ vSizes[t2b->getVertex(2)->getIndex()]);
+ double lcBGM2 = 0.3333333333 * (vSizesBGM[t2b->getVertex(0)->getIndex()] +
+ vSizesBGM[t2b->getVertex(1)->getIndex()] +
+ vSizesBGM[t2b->getVertex(2)->getIndex()]);
+ double lc3 = 0.3333333333 * (vSizes[t3b->getVertex(0)->getIndex()] +
+ vSizes[t3b->getVertex(1)->getIndex()] +
+ vSizes[t3b->getVertex(2)->getIndex()]);
+ double lcBGM3 = 0.3333333333 * (vSizesBGM[t3b->getVertex(0)->getIndex()] +
+ vSizesBGM[t3b->getVertex(1)->getIndex()] +
+ vSizesBGM[t3b->getVertex(2)->getIndex()]);
+ double lc4 = 0.3333333333 * (vSizes[t4b->getVertex(0)->getIndex()] +
+ vSizes[t4b->getVertex(1)->getIndex()] +
+ vSizes[t4b->getVertex(2)->getIndex()]);
+ double lcBGM4 = 0.3333333333 * (vSizesBGM[t4b->getVertex(0)->getIndex()] +
+ vSizesBGM[t4b->getVertex(1)->getIndex()] +
+ vSizesBGM[t4b->getVertex(2)->getIndex()]);
MTri3 *t1b3 = new MTri3(t1b, Extend1dMeshIn2dSurfaces() ? std::min(lc1, lcBGM1) : lcBGM1);
MTri3 *t2b3 = new MTri3(t2b, Extend1dMeshIn2dSurfaces() ? std::min(lc2, lcBGM2) : lcBGM2);
MTri3 *t3b3 = new MTri3(t3b, Extend1dMeshIn2dSurfaces() ? std::min(lc3, lcBGM3) : lcBGM3);
@@ -678,14 +681,14 @@ bool gmshVertexCollapse(const double lMin, MTri3 *t1, GFace *gf,
MVertex *v1 = ring[(iMin + 1 + i) % ring.size()];
MVertex *v2 = ring[(iMin + 1 + i + 1) % ring.size()];
MTriangle *t = new MTriangle(v1,v2,ring[iMin]);
- double lc = 0.3333333333 * (vSizes[t->getVertex(0)->getNum()] +
- vSizes[t->getVertex(1)->getNum()] +
- vSizes[t->getVertex(2)->getNum()]);
- double lcBGM = 0.3333333333 * (vSizesBGM[t->getVertex(0)->getNum()] +
- vSizesBGM[t->getVertex(1)->getNum()] +
- vSizesBGM[t->getVertex(2)->getNum()]);
+ double lc = 0.3333333333 * (vSizes[t->getVertex(0)->getIndex()] +
+ vSizes[t->getVertex(1)->getIndex()] +
+ vSizes[t->getVertex(2)->getIndex()]);
+ double lcBGM = 0.3333333333 * (vSizesBGM[t->getVertex(0)->getIndex()] +
+ vSizesBGM[t->getVertex(1)->getIndex()] +
+ vSizesBGM[t->getVertex(2)->getIndex()]);
MTri3 *t3 = new MTri3(t, Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM);
- // printf("Creation %p = %d %d %d\n",t3,v1->getNum(),v2->getNum(),ring[iMin]->getNum());
+ // printf("Creation %p = %d %d %d\n",t3,v1->getIndex(),v2->getIndex(),ring[iMin]->getIndex());
outside.push_back(t3);
newTris.push_back(t3);
}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 53b1a4fc0b..4e3b8fa0b6 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -60,7 +60,7 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
in.pointlist[(I-1)*3 + 0] = (*itv)->x();
in.pointlist[(I-1)*3 + 1] = (*itv)->y();
in.pointlist[(I-1)*3 + 2] = (*itv)->z();
- (*itv)->setNum(I++);
+ (*itv)->setIndex(I++);
numberedV.push_back(*itv);
++itv;
}
@@ -93,9 +93,9 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
tetgenio::init(p);
p->numberofvertices = 3;
p->vertexlist = new int[p->numberofvertices];
- p->vertexlist[0] = t->getVertex(0)->getNum();
- p->vertexlist[1] = t->getVertex(1)->getNum();
- p->vertexlist[2] = t->getVertex(2)->getNum();
+ p->vertexlist[0] = t->getVertex(0)->getIndex();
+ p->vertexlist[1] = t->getVertex(1)->getIndex();
+ p->vertexlist[2] = t->getVertex(2)->getIndex();
in.facetmarkerlist[I] = gf->tag();
++I;
}
@@ -155,7 +155,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
v[2] = numberedV[out.trifacelist[i * 3 + 2] - 1];
GFace *gf = gr->model()->getFaceByTag(out.trifacemarkerlist[i]);
- double guess[2] = {0,0};
+ double guess[2] = {0, 0};
int Count = 0;
for(int j = 0; j < 3; j++){
if(!v[j]->onWhat()){
@@ -163,8 +163,8 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
}
else if(v[j]->onWhat()->dim() == 2){
double uu,vv;
- v[j]->getParameter(0,uu);
- v[j]->getParameter(1,vv);
+ v[j]->getParameter(0, uu);
+ v[j]->getParameter(1, vv);
guess[0] += uu;
guess[1] += vv;
Count++;
@@ -175,22 +175,22 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
v[j]->getParameter(0, UU);
SPoint2 param;
param = ge->reparamOnFace(gf, UU, 1);
- guess[0]+=param.x();
- guess[1]+=param.y();
+ guess[0] += param.x();
+ guess[1] += param.y();
Count++;
}
else if (v[j]->onWhat()->dim() == 0){
SPoint2 param;
GVertex *gv = (GVertex*)v[j]->onWhat();
param = gv->reparamOnFace(gf,1);
- guess[0]+=param.x();
- guess[1]+=param.y();
+ guess[0] += param.x();
+ guess[1] += param.y();
Count++;
}
}
if (Count != 0){
- guess[0]/=Count;
- guess[1]/=Count;
+ guess[0] /= Count;
+ guess[1] /= Count;
}
for(int j = 0; j < 3; j++){
@@ -202,28 +202,21 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
MVertex *v1b;
if (CTX::instance()->mesh.order > 1 &&
CTX::instance()->mesh.secondOrderExperimental){
- // PARAMETRIC COORDINATES SHOULD BE SET for the vertex !!!!!!!!!!!!!
- // This is not 100 % safe yet, so we reserve that operation for high order
- // meshes.
- // Msg::Debug("A new point has been inserted in mesh face %d by the 3D mesher, guess %g %g",gf->tag(),guess[0],guess[1]);
- GPoint gp = gf->closestPoint (SPoint3(v[j]->x(), v[j]->y(), v[j]->z()),guess);
-
- // To be safe, we should ensure that this mesh motion does not lead to an invalid mesh !!!!
- //v1b = new MVertex(gp.x(),gp.y(),gp.z(),gf);
- if (gp.g() ){
- //v1b = new MFaceVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf,gp.u(),gp.v());
- v1b = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
- // Msg::Info("The point has been projected back to the surface (%g %g %g) -> (%g %g %g)",
- // v[j]->x(), v[j]->y(), v[j]->z(),gp.x(),gp.y(),gp.z());
+ // parametric coordinates should be set for the vertex !
+ // (this is not 100 % safe yet, so we reserve that operation
+ // for high order meshes)
+ GPoint gp = gf->closestPoint(SPoint3(v[j]->x(), v[j]->y(), v[j]->z()),guess);
+ if (gp.g()){
+ v1b = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
}
else{
- v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
+ v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
Msg::Warning("The point was not projected back to the surface (%g %g %g)",
v[j]->x(), v[j]->y(), v[j]->z());
}
}
else{
- v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
+ v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
}
gf->mesh_vertices.push_back(v1b);
@@ -348,7 +341,7 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
tmp[0] = (*itv)->x();
tmp[1] = (*itv)->y();
tmp[2] = (*itv)->z();
- (*itv)->setNum(I++);
+ (*itv)->setIndex(I++);
numberedV.push_back(*itv);
Ng_AddPoint(ngmesh, tmp);
++itv;
@@ -360,7 +353,7 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
tmp[0] = gr->mesh_vertices[i]->x();
tmp[1] = gr->mesh_vertices[i]->y();
tmp[2] = gr->mesh_vertices[i]->z();
- gr->mesh_vertices[i]->setNum(I++);
+ gr->mesh_vertices[i]->setIndex(I++);
Ng_AddPoint(ngmesh, tmp);
}
}
@@ -372,9 +365,9 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
for(unsigned int i = 0; i< gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
int tmp[3];
- tmp[0] = t->getVertex(0)->getNum();
- tmp[1] = t->getVertex(1)->getNum();
- tmp[2] = t->getVertex(2)->getNum();
+ tmp[0] = t->getVertex(0)->getIndex();
+ tmp[1] = t->getVertex(1)->getIndex();
+ tmp[2] = t->getVertex(2)->getIndex();
Ng_AddSurfaceElement(ngmesh, NG_TRIG, tmp);
}
++it;
@@ -386,10 +379,10 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
// netgen expects tet with negative volume
if(t->getVolumeSign() > 0) t->revert();
int tmp[4];
- tmp[0] = t->getVertex(0)->getNum();
- tmp[1] = t->getVertex(1)->getNum();
- tmp[2] = t->getVertex(2)->getNum();
- tmp[3] = t->getVertex(3)->getNum();
+ tmp[0] = t->getVertex(0)->getIndex();
+ tmp[1] = t->getVertex(1)->getIndex();
+ tmp[2] = t->getVertex(2)->getIndex();
+ tmp[3] = t->getVertex(3)->getIndex();
Ng_AddVolumeElement(ngmesh, NG_TET, tmp);
}
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 51bdb2843c..96fd06270a 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -722,7 +722,7 @@ void insertVerticesInRegion (GRegion *gr)
setLcs(gr->tetrahedra[i], vSizesMap);
for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
- it->first->setNum(NUM++);
+ it->first->setIndex(NUM++);
vSizes.push_back(it->second);
vSizesBGM.push_back(it->second);
}
@@ -799,12 +799,12 @@ void insertVerticesInRegion (GRegion *gr)
worst->tet()->xyz2uvw(center, uvw);
if(worst->tet()->isInside(uvw[0], uvw[1], uvw[2])){
MVertex *v = new MVertex(center[0], center[1], center[2], worst->onWhat());
- v->setNum(NUM++);
+ v->setIndex(NUM++);
double lc1 =
- (1 - uvw[0] - uvw[1] - uvw[2]) * vSizes[worst->tet()->getVertex(0)->getNum()] +
- uvw[0] * vSizes[worst->tet()->getVertex(1)->getNum()] +
- uvw[1] * vSizes[worst->tet()->getVertex(2)->getNum()] +
- uvw[2] * vSizes[worst->tet()->getVertex(3)->getNum()];
+ (1 - uvw[0] - uvw[1] - uvw[2]) * vSizes[worst->tet()->getVertex(0)->getIndex()] +
+ uvw[0] * vSizes[worst->tet()->getVertex(1)->getIndex()] +
+ uvw[1] * vSizes[worst->tet()->getVertex(2)->getIndex()] +
+ uvw[2] * vSizes[worst->tet()->getVertex(3)->getIndex()];
double lc = BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]);
// double lc = std::min(lc1, BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]));
vSizes.push_back(lc1);
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index f65fdc3bff..9d3e4c10ac 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -114,14 +114,14 @@ class MTet4
const double dy = base->getVertex(0)->y() - center[1];
const double dz = base->getVertex(0)->z() - center[2];
circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
- double lc1 = 0.25*(sizes[base->getVertex(0)->getNum()]+
- sizes[base->getVertex(1)->getNum()]+
- sizes[base->getVertex(2)->getNum()]+
- sizes[base->getVertex(3)->getNum()]);
- double lcBGM = 0.25*(sizesBGM[base->getVertex(0)->getNum()]+
- sizesBGM[base->getVertex(1)->getNum()]+
- sizesBGM[base->getVertex(2)->getNum()]+
- sizesBGM[base->getVertex(3)->getNum()]);
+ double lc1 = 0.25*(sizes[base->getVertex(0)->getIndex()]+
+ sizes[base->getVertex(1)->getIndex()]+
+ sizes[base->getVertex(2)->getIndex()]+
+ sizes[base->getVertex(3)->getIndex()]);
+ double lcBGM = 0.25*(sizesBGM[base->getVertex(0)->getIndex()]+
+ sizesBGM[base->getVertex(1)->getIndex()]+
+ sizesBGM[base->getVertex(2)->getIndex()]+
+ sizesBGM[base->getVertex(3)->getIndex()]);
double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
circum_radius /= lc;
deleted = false;
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 5cc0acef97..a7eb3c3297 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -175,7 +175,6 @@ double qmTet(const double &x1, const double &y1, const double &z1,
}
}
-
double mesh_functional_distorsion(MTriangle *t, double u, double v)
{
// compute uncurved element jacobian d_u x and d_v x
@@ -208,7 +207,6 @@ double mesh_functional_distorsion(MTriangle *t, double u, double v)
return det;
}
-
double mesh_functional_distorsion_p2(MTriangle *t)
{
double d = mesh_functional_distorsion(t,0.,0.);
@@ -220,8 +218,6 @@ double mesh_functional_distorsion_p2(MTriangle *t)
return d;
}
-
-
double qmDistorsionOfMapping (MTriangle *e)
{
// return 1.0;
@@ -253,49 +249,47 @@ static double mesh_functional_distorsion(MTetrahedron *t, double u, double v, do
{
// compute uncurved element jacobian d_u x and d_v x
double mat[3][3];
- t->getPrimaryJacobian(u,v,w, mat);
+ t->getPrimaryJacobian(u, v, w, mat);
const double det1 = det3x3(mat);
//const double det1 = t->getJacobian(u,v,w,mat);
// const double det1 = det3x3(mat);
- t->getJacobian(u,v,w,mat);
+ t->getJacobian(u, v, w, mat);
const double detN = det3x3(mat);
// const double detN = det3x3(mat);
// printf("%g %g %g = %g %g\n",u,v,w,det1,detN);
if (det1 == 0 || detN == 0) return 0;
- double dist = std::min(detN/det1, det1/detN);
+ double dist = std::min(detN / det1, det1 / detN);
return dist;
}
-
-double qmDistorsionOfMapping (MTetrahedron *e)
+double qmDistorsionOfMapping(MTetrahedron *e)
{
- if (e->getPolynomialOrder() == 1)return 1.0;
+ if (e->getPolynomialOrder() == 1) return 1.0;
IntPt *pts;
int npts;
- e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
+ e->getIntegrationPoints(e->getPolynomialOrder(), &npts, &pts);
double dmin;
- for (int i=0;i<npts;i++){
+ for (int i = 0; i < npts; i++){
const double u = pts[i].pt[0];
const double v = pts[i].pt[1];
const double w = pts[i].pt[2];
- const double di = mesh_functional_distorsion (e,u,v,w);
- dmin = (i==0)? di : std::min(dmin,di);
+ const double di = mesh_functional_distorsion(e, u, v, w);
+ dmin = (i == 0) ? di : std::min(dmin, di);
}
const gmshMatrix<double>& points = e->getFunctionSpace()->points;
- for (int i=0;i<e->getNumPrimaryVertices();i++) {
- const double u = points(i,0);
- const double v = points(i,1);
- const double w = points(i,2);
- const double di = mesh_functional_distorsion (e,u,v,w);
- dmin = std::min(dmin,di);
+ for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
+ const double u = points(i, 0);
+ const double v = points(i, 1);
+ const double w = points(i, 2);
+ const double di = mesh_functional_distorsion(e, u, v, w);
+ dmin = std::min(dmin, di);
}
- // printf("DMIN = %g\n\n",dmin);
return dmin;
}
--
GitLab