From f5a5020298b465b7c203ed16de13cec0297addc0 Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Wed, 12 Jul 2017 07:47:43 +0200 Subject: [PATCH] fix warnings --- Mesh/delaunay3d.cpp | 50 ++++++++-------- Mesh/delaunay3d_private.h | 13 ++-- Mesh/delaunay_refinement.cpp | 88 +++++++++++++++------------- Mesh/meshGRegionBoundaryRecovery.cpp | 74 +++++++++++++---------- 4 files changed, 123 insertions(+), 102 deletions(-) diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp index f47a7eb60a..c26d282d4a 100644 --- a/Mesh/delaunay3d.cpp +++ b/Mesh/delaunay3d.cpp @@ -370,14 +370,14 @@ static void computeNeighboringTetsOfACavity(const std::vector<Tet*> &cavity, void printtet (const char *c, Tet *t){ printf("%s ",c); - + if (t->V[0]) printf("%7d %7d %7d %7d \n", t->V[0]->getNum(), t->V[1]->getNum(), t->V[2]->getNum(), t->V[3]->getNum()); - + } static bool buildVertexCavity(Vert *v, std::vector<Tet*> &cavity, std::stack<Tet*> &_stack) @@ -421,12 +421,12 @@ static bool buildVertexCavity(Vert *v, std::vector<Tet*> &cavity, std::stack<Tet return false; } } - } + } } } return true; } - + static bool buildEdgeCavity(Tet *t, int iLocalEdge, Vert **v1, Vert **v2, @@ -439,16 +439,16 @@ static bool buildEdgeCavity(Tet *t, int iLocalEdge, *v1 = t->V[edges[iLocalEdge][0]]; *v2 = t->V[edges[iLocalEdge][1]]; - + // the 5 - i th edge contains the other 2 points of the tet Vert *lastinring = t->V[edges[5 - iLocalEdge][0]]; - + ring.push_back(lastinring); cavity.push_back(t); // printf("edge %d %d \n",(*v1)->getNum(),(*v2)->getNum()); /// printtet ("first", t); - + while (1){ Vert *ov1 = t->V[edges[5 - iLocalEdge][0]]; Vert *ov2 = t->V[edges[5 - iLocalEdge][1]]; @@ -543,9 +543,9 @@ static bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread) int p2 = sp.triangles[i][1]; int p3 = sp.triangles[i][2]; tetQuality1[i] = tetQuality (ring[p1],ring[p2],ring[p3],v1,&volume1[i]); - tetQuality2[i] = tetQuality (ring[p3],ring[p2],ring[p1],v2,&volume2[i]); + tetQuality2[i] = tetQuality (ring[p3],ring[p2],ring[p1],v2,&volume2[i]); } - + // look for the best triangulation, i.e. the one that maximize the // minimum element quality double minQuality[100]; @@ -612,7 +612,7 @@ static bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread) computeAdjacencies (outside[i],j,ctnr); } } - + return true; } @@ -651,7 +651,7 @@ static bool relocateVertex (Vert* v, std::vector<Tet*> &cavity){ } v->x()=oldX, v->y()=oldY, v->z()=oldZ; return false; - + } static bool relocateVertex (Vert* v, std::stack<Tet*> &_work, std::vector<Tet*> &cavity){ @@ -676,7 +676,7 @@ void vertexRelocationPass (int numThreads, std::vector<Vert*> &v){ if (!N)break; if (iter++ >= 0)break; } - + } @@ -687,14 +687,14 @@ void edgeSwapPass (int numThreads, tetContainer &allocator, edgeContainer &embed int N = 0; for (int myThread=0; myThread < numThreads; myThread++) { int NNN = allocator.size(myThread); - for (unsigned int i=0;i<NNN;i++){ - // printf("%d %d\n",NNN,i); + for (int i=0;i<NNN;i++){ + // printf("%d %d\n",NNN,i); Tet *t = allocator(myThread,i); if (t->V[0]){ double vol; if (tetQuality (t,&vol) < 0.3){ for (int j=0;j<6;j++){ - int iLocalEdge = j; + int iLocalEdge = j; Edge e (t->V[edges[iLocalEdge][0]], t->V[edges[iLocalEdge][1]]); if (! embeddedEdges.find(e) && !ec.find(e)){ ec.addNewEdge(e); @@ -706,7 +706,7 @@ void edgeSwapPass (int numThreads, tetContainer &allocator, edgeContainer &embed } } else{ - // printf("cannot swop\n"); + // printf("cannot swap\n"); } } } @@ -869,7 +869,7 @@ static bool fixDelaunayCavity (Vert *v, starShapeness (v, bndK, _negatives); if (_negatives.empty())return false; - + // unset all tets of the cavity for (unsigned int i=0; i< cavity.size(); i++)cavity[i]->unset(myThread,K); for (unsigned int i=0; i<bndK.size(); i++)if(bndK[i].t)bndK[i].t->unset(myThread,K); @@ -882,7 +882,7 @@ static bool fixDelaunayCavity (Vert *v, Tet *containsV = tetContainsV (v,cavity); if (! containsV) return true; - + while (!_negatives.empty()) { for (unsigned int i=0;i<_negatives.size();i++){ conn &c = bndK[_negatives[i] ]; @@ -1091,11 +1091,11 @@ static Tet* randomTet (int thread, tetContainer &allocator) } -int isCavityCompatibleWithEmbeddedEdges(cavityContainer &cavity, - connContainer &bndK, +int isCavityCompatibleWithEmbeddedEdges(cavityContainer &cavity, + connContainer &bndK, edgeContainer &allEmbeddedEdges){ // return true; - + const unsigned int bSize = bndK.size(); std::vector<Edge> ed; for (unsigned int i=0; i<bSize; i++) { @@ -1105,7 +1105,7 @@ int isCavityCompatibleWithEmbeddedEdges(cavityContainer &cavity, } } } - + for (unsigned int i=0; i<cavity.size(); i++){ for (int j=0;j<6;j++){ Edge e = cavity[i]->getEdge(j); @@ -1115,7 +1115,7 @@ int isCavityCompatibleWithEmbeddedEdges(cavityContainer &cavity, } } return 1; -} +} //#define _VERBOSE 1 @@ -1186,7 +1186,7 @@ void delaunayTrgl (const unsigned int numThreads, std::vector<Vert*> vToAdd(NPTS_AT_ONCE); // printf("reaching parallel section\n"); - + #if defined(_OPENMP) #pragma omp barrier #endif @@ -1325,7 +1325,7 @@ void delaunayTrgl (const unsigned int numThreads, for (unsigned int myThread=0; myThread < numThreads;myThread++) for (unsigned int i=0;i<allocator.size(myThread);i++)allocator(myThread,i)->setAllDeleted(); - + } static void initialCube (std::vector<Vert*> &v, diff --git a/Mesh/delaunay3d_private.h b/Mesh/delaunay3d_private.h index 0c7871bb8f..fd247ff4d1 100644 --- a/Mesh/delaunay3d_private.h +++ b/Mesh/delaunay3d_private.h @@ -155,15 +155,15 @@ struct edgeContainer { _size = 0; _hash.resize(N); - _size_obj = sizeof(Edge); + _size_obj = sizeof(Edge); } inline size_t H (const Edge &e) const { const size_t h = ((size_t)e.first) ; // printf("%lu %lu %lu %lu\n",h,(h/2)%_hash.size(),h/64,h>>6); - return (h/_size_obj) %_hash.size(); + return (h/_size_obj) %_hash.size(); } - + inline bool find (const Edge &e) const { const std::vector<Edge> &v = _hash[H(e)]; for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return true;} @@ -171,7 +171,7 @@ struct edgeContainer } bool empty () const {return _size == 0;} - + bool addNewEdge (const Edge &e) { std::vector<Edge> &v = _hash[H(e)]; @@ -360,8 +360,9 @@ public: class tetContainer { std::vector<aBunchOfStuff<Tet> *> _perThread; public: - unsigned int size(int thread) const { - if (_perThread.size() <= thread)return 0; + unsigned int size(int thread) const + { + if ((int)_perThread.size() <= thread)return 0; return _perThread[thread]->size(); } inline Tet * operator () (int thread, int j) const diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp index 2be3126128..cb61906f72 100644 --- a/Mesh/delaunay_refinement.cpp +++ b/Mesh/delaunay_refinement.cpp @@ -47,7 +47,7 @@ static double GMSHSIZE (const SPoint3 &p, Field *f, double lc_nodal) { if (!CTX::instance()->mesh.lcExtendFromBoundary) lc_nodal = 1.e22; // double _a = Cpu(); if (f)f_field = (*f)(p.x(), p.y(), p.z()); - // _C3 += Cpu() - _a; + // _C3 += Cpu() - _a; double lc = std::min (lc_nodal,f_field); lc = std::max(lc, CTX::instance()->mesh.lcMin); @@ -74,13 +74,13 @@ double adaptiveTrapezoidalRule(SPoint3 p1 , SPoint3 p2 , // value of f on both sides double f1 = lc1;//GMSHSIZE(p1, bgm, lc1); //f(p1 + dp*t1,data); - double f2 = lc2;//GMSHSIZE(p2, bgm, lc2); //f(p1 + dp*t2,data); - + double f2 = lc2;//GMSHSIZE(p2, bgm, lc2); //f(p1 + dp*t2,data); + dl = p1.distance(p2); // adim_lenght of half the edge should be bigger than EXCLUSION_FACTOR // +------x--------+ - + if (dl / (2*std::min(lc1,lc2)) <= EXCLUSION_FACTOR){ // printf ("edge length %g lc %g %g\n",dl,f1,f2); return EXCLUSION_FACTOR; @@ -140,15 +140,15 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp) Field *bgm = NULL; if(fields->getBackgroundField() > 0){ bgm = fields->get(fields->getBackgroundField()); - } + } // double _a = Cpu(); - + // printf("%g %g \n",e.first->lc(), e.second->lc()); - + const double dN = adaptiveTrapezoidalRule(p1,p2,e.first->lc(), e.second->lc(), _result, dl, temp, bgm); - // _C1 += Cpu() - _a; + // _C1 += Cpu() - _a; // _a = Cpu(); const int N = (int) (dN+0.1); @@ -191,7 +191,7 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp) } } } - // _C2 += Cpu() - _a; + // _C2 += Cpu() - _a; // printf(" press enter\n"); // getchar(); // printf("%d points added\n",S.size()); @@ -199,12 +199,13 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp) // exit(1); } -inline bool distributeEdgeThroughThreads ( int nbThreads , int myThread , - const Edge &ed){ +inline bool distributeEdgeThroughThreads(int nbThreads, int myThread, + const Edge &ed) +{ // const size_t h = ((size_t)ed.first->getNum()); const size_t h = ((size_t)ed.first) >> 3; // printf("%lu %d %d %d\n",h,nbThreads,myThread,h % nbThreads); - return h % nbThreads == myThread; + return (int)(h % nbThreads) == myThread; } void saturateEdges(edgeContainer &ec, @@ -212,8 +213,9 @@ void saturateEdges(edgeContainer &ec, int nbThreads, std::vector<Vert*> &S) { - -#pragma omp parallel +#ifdef _OPENMP +#pragma omp parallel +#endif { std::vector<Vert*> Sloc; std::stack<IPT> temp; @@ -222,7 +224,7 @@ void saturateEdges(edgeContainer &ec, // nbThreads = omp_get_num_threads(); #else int myThread = 0; -#endif +#endif for (int iThread = 0; iThread<nbThreads;iThread++){ const int N = T.size(iThread); for (int i=0;i<N;i++){ @@ -241,9 +243,13 @@ void saturateEdges(edgeContainer &ec, } } } +#ifdef _OPENMP #pragma omp critical - S.insert (S.end(),Sloc.begin(), Sloc.end()); +#endif + S.insert (S.end(),Sloc.begin(), Sloc.end()); +#ifdef _OPENMP #pragma omp barrier +#endif } for (int iThread = 0; iThread<nbThreads;iThread++){ @@ -439,7 +445,7 @@ void disconnectEmbeddedFaces (GRegion *gr, connSet &faceToTet, std::vector<Vert Tet *t1 = itf->t; Tet *t2 = itf->t->T[itf->i]; // if (!(f == t1->getFace(itf->i)))printf("aie\n"); - for (int k=0;k<4;k++){ + for (int k=0;k<4;k++){ if (t1 && t1->T[k] == t2){t1->T[k]=NULL;ok1=true;} if (t2 && t2->T[k] == t1){t2->T[k]=NULL;ok2=true;} } @@ -455,8 +461,8 @@ void computeMeshSizes (GRegion *gr, std::map<MVertex*, double> &s){ for (std::list<GEdge*>::iterator it = e.begin() ; it != e.end(); ++it){ for (unsigned int i = 0; i < (*it)->lines.size(); i++){ double d = distance((*it)->lines[i]->getVertex(0), (*it)->lines[i]->getVertex(1)); - updateSize ((*it)->lines[i]->getVertex(0), s, d); - updateSize ((*it)->lines[i]->getVertex(1), s, d); + updateSize ((*it)->lines[i]->getVertex(0), s, d); + updateSize ((*it)->lines[i]->getVertex(1), s, d); } } std::list<GFace*> f ; @@ -485,31 +491,31 @@ void computeMeshSizes (GRegion *gr, std::map<MVertex*, double> &s){ } -void parallelDelaunay (int NT, std::vector<Vert*> &S, - tetContainer &allocator, - int iter, - bool explicitFiltering, +void parallelDelaunay (int NT, std::vector<Vert*> &S, + tetContainer &allocator, + int iter, + bool explicitFiltering, std::vector<int> &indices, - edgeContainer *embeddedEdges){ + edgeContainer *embeddedEdges){ // Msg::Info ("Parallel Delaunay with %d threads",NT); int N = S.size(); NT = std::min(NT, MAX_NUM_THREADS_); std::vector<Vert*> assignTo0[1]; std::vector<std::vector<Vert*> > assignTo (NT); - - for (int i=1;i<indices.size();i++){ + + for (unsigned int i = 1; i < indices.size(); i++){ int start = indices[i-1]; int end = indices[i]; int sizePerBlock = (NT*((end-start) / NT))/NT; int currentBlock = 0; int localCounter = 0; - // printf("sizePerBlock[%d] = %d (%d,%d)\n",i,sizePerBlock, start, end); + // printf("sizePerBlock[%d] = %d (%d,%d)\n",i,sizePerBlock, start, end); if (i < 2){ for (int jPt=start;jPt<end;jPt++){ assignTo0[0].push_back(S[jPt]); - } + } } else { for (int jPt=start;jPt<end;jPt++){ @@ -534,10 +540,10 @@ void edgeBasedRefinement(const int numThreads, { // fill up old Datastructures edgeContainer embeddedEdges (10000); - + std::map<MVertex*, double> sizes; - computeMeshSizes (gr, sizes); - + computeMeshSizes (gr, sizes); + tetContainer allocator (numThreads,1000000); SBoundingBox3d bb; @@ -655,7 +661,7 @@ void edgeBasedRefinement(const int numThreads, add_all.insert (add_all.end(), add.begin(), add.end()); size_t sss = 0; for (int myThread=0; myThread < numThreads; myThread++)sss+= allocator.size(myThread); - + Msg::Info("IT %3d %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d", iter,add.size(), (t2-t1), @@ -669,21 +675,21 @@ void edgeBasedRefinement(const int numThreads, } std::vector<Vert*> vv; - for (int myThread=0; myThread < numThreads; myThread++) + for (int myThread=0; myThread < numThreads; myThread++) for (unsigned int i=0;i<allocator.size(myThread);i++) for (unsigned int j=0;j<4;j++) if (allocator(myThread,i)->V[j]) - allocator(myThread,i)->V[j]->setNum(0); - - for (int myThread=0; myThread < numThreads; myThread++) + allocator(myThread,i)->V[j]->setNum(0); + + for (int myThread=0; myThread < numThreads; myThread++) for (unsigned int i=0;i<allocator.size(myThread);i++) for (unsigned int j=0;j<4;j++) - if (allocator(myThread,i)->V[j] && allocator(myThread,i)->V[j]->getNum() == 0){ + if (allocator(myThread,i)->V[j] && allocator(myThread,i)->V[j]->getNum() == 0){ allocator(myThread,i)->V[j]->setNum(1); std::map<Vert*,MVertex*>::iterator it = _ma.find(allocator(myThread,i)->V[j]); if (it == _ma.end())vv.push_back(allocator(myThread,i)->V[j]); } - + double t6 = Cpu(); Msg::Info("Optimizing"); edgeSwapPass (numThreads, allocator, embeddedEdges); @@ -698,7 +704,7 @@ void edgeBasedRefinement(const int numThreads, int cat [10] = {0,0,0,0,0,0,0,0,0,0}; double MIN = 1.0; - + for (int myThread=0; myThread < numThreads; myThread++) { for (unsigned int i=0; i< allocator.size(myThread);i++){ Tet *tt = allocator (myThread,i); @@ -728,7 +734,7 @@ void edgeBasedRefinement(const int numThreads, for (int i=0;i<10;i++){ Msg::Info ("Tet Quality [%4.3f,%4.3f] %8d",0.1*i,0.1*(i+1),cat[i]); } - + if (Msg::GetVerbosity() == 99) { std::map<Edge,double> _sizes; for (unsigned int i=0; i< allocator.size(0);i++){ @@ -759,5 +765,5 @@ void edgeBasedRefinement(const int numThreads, exp (sum / _sizes.size()),nbBad,_sizes.size()); } - + } diff --git a/Mesh/meshGRegionBoundaryRecovery.cpp b/Mesh/meshGRegionBoundaryRecovery.cpp index 4e800df262..fdbb78dbb7 100644 --- a/Mesh/meshGRegionBoundaryRecovery.cpp +++ b/Mesh/meshGRegionBoundaryRecovery.cpp @@ -742,7 +742,7 @@ bool tetgenmesh::reconstructmesh(void *p) if (!_extras.empty())Msg::Info("We add %d steiner points...",_extras.size()); - + if (l_edges.size() > 0) { @@ -761,7 +761,7 @@ bool tetgenmesh::reconstructmesh(void *p) } } assert(ge != NULL); - // Msg::Info("Steiner points exist on GEdge %d", ge->tag()); + // Msg::Info("Steiner points exist on GEdge %d", ge->tag()); // Delete the old triangles. for(unsigned int i = 0; i < ge->lines.size(); i++) delete ge->lines[i]; @@ -775,11 +775,13 @@ bool tetgenmesh::reconstructmesh(void *p) if (shellmark(segloop) == etag) { p[0] = sorg(segloop); p[1] = sdest(segloop); - MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])]; - MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])]; - // MVertex *v2 = _vertices[pointmark(p[1])]; - // printf("%d %d %d\n",pointmark(p[0]),pointmark(p[1]),_vertices.size()); - // printf("%g %g %g %g\n",v1->x(),v1->y(),v2->x(),v2->y()); + MVertex *v1 = pointmark(p[0]) >= (int)_vertices.size() ? + _extras[pointmark(p[0])] : _vertices[pointmark(p[0])]; + MVertex *v2 = pointmark(p[1]) >= (int)_vertices.size() ? + _extras[pointmark(p[1])] : _vertices[pointmark(p[1])]; + // MVertex *v2 = _vertices[pointmark(p[1])]; + // printf("%d %d %d\n",pointmark(p[0]),pointmark(p[1]),_vertices.size()); + // printf("%g %g %g %g\n",v1->x(),v1->y(),v2->x(),v2->y()); MLine *t = new MLine(v1, v2); ge->lines.push_back(t); } @@ -822,15 +824,21 @@ bool tetgenmesh::reconstructmesh(void *p) p[0] = sorg(subloop); p[1] = sdest(subloop); p[2] = sapex(subloop); - MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])]; - MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])]; - MVertex *v3 = pointmark(p[2]) >= _vertices.size() ? _extras[pointmark(p[2])] : _vertices[pointmark(p[2])]; - // if (pointmark(p[0]) >= _vertices.size())printf("F %d %d\n",v1->getIndex(),pointmark(p[0])); - // if (pointmark(p[1]) >= _vertices.size())printf("F %d %d\n",v2->getIndex(),pointmark(p[1])); - // if (pointmark(p[2]) >= _vertices.size())printf("F %d %d\n",v3->getIndex(),pointmark(p[2])); - // MVertex *v1 = _vertices[pointmark(p[0])]; - // MVertex *v2 = _vertices[pointmark(p[1])]; - // MVertex *v3 = _vertices[pointmark(p[2])]; + MVertex *v1 = pointmark(p[0]) >= (int)_vertices.size() ? + _extras[pointmark(p[0])] : _vertices[pointmark(p[0])]; + MVertex *v2 = pointmark(p[1]) >= (int)_vertices.size() ? + _extras[pointmark(p[1])] : _vertices[pointmark(p[1])]; + MVertex *v3 = pointmark(p[2]) >= (int)_vertices.size() ? + _extras[pointmark(p[2])] : _vertices[pointmark(p[2])]; + // if (pointmark(p[0]) >= _vertices.size()) + // printf("F %d %d\n",v1->getIndex(),pointmark(p[0])); + // if (pointmark(p[1]) >= _vertices.size()) + // printf("F %d %d\n",v2->getIndex(),pointmark(p[1])); + // if (pointmark(p[2]) >= _vertices.size()) + // printf("F %d %d\n",v3->getIndex(),pointmark(p[2])); + // MVertex *v1 = _vertices[pointmark(p[0])]; + // MVertex *v2 = _vertices[pointmark(p[1])]; + // MVertex *v3 = _vertices[pointmark(p[2])]; MTriangle *t = new MTriangle(v1, v2, v3); gf->triangles.push_back(t); } @@ -851,20 +859,26 @@ bool tetgenmesh::reconstructmesh(void *p) p[2] = apex(tetloop); p[3] = oppo(tetloop); - MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])]; - MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])]; - MVertex *v3 = pointmark(p[2]) >= _vertices.size() ? _extras[pointmark(p[2])] : _vertices[pointmark(p[2])]; - MVertex *v4 = pointmark(p[3]) >= _vertices.size() ? _extras[pointmark(p[3])] : _vertices[pointmark(p[3])]; - - // if (pointmark(p[0]) >= _vertices.size())printf("%d %d\n",v1->getIndex(),pointmark(p[0])); - // if (pointmark(p[1]) >= _vertices.size())printf("%d %d\n",v2->getIndex(),pointmark(p[1])); - // if (pointmark(p[2]) >= _vertices.size())printf("%d %d\n",v3->getIndex(),pointmark(p[2])); - // if (pointmark(p[3]) >= _vertices.size())printf("%d %d\n",v4->getIndex(),pointmark(p[3])); - - // MVertex *v1 = _vertices[pointmark(p[0])]; - // MVertex *v2 = _vertices[pointmark(p[1])]; - // MVertex *v3 = _vertices[pointmark(p[2])]; - // MVertex *v4 = _vertices[pointmark(p[3])]; + MVertex *v1 = pointmark(p[0]) >= (int)_vertices.size() ? + _extras[pointmark(p[0])] : _vertices[pointmark(p[0])]; + MVertex *v2 = pointmark(p[1]) >= (int)_vertices.size() ? + _extras[pointmark(p[1])] : _vertices[pointmark(p[1])]; + MVertex *v3 = pointmark(p[2]) >= (int)_vertices.size() ? + _extras[pointmark(p[2])] : _vertices[pointmark(p[2])]; + MVertex *v4 = pointmark(p[3]) >= (int)_vertices.size() ? + _extras[pointmark(p[3])] : _vertices[pointmark(p[3])]; + // if (pointmark(p[0]) >= _vertices.size()) + // printf("%d %d\n",v1->getIndex(),pointmark(p[0])); + // if (pointmark(p[1]) >= _vertices.size()) + // printf("%d %d\n",v2->getIndex(),pointmark(p[1])); + // if (pointmark(p[2]) >= _vertices.size()) + // printf("%d %d\n",v3->getIndex(),pointmark(p[2])); + // if (pointmark(p[3]) >= _vertices.size()) + /// printf("%d %d\n",v4->getIndex(),pointmark(p[3])); + // MVertex *v1 = _vertices[pointmark(p[0])]; + // MVertex *v2 = _vertices[pointmark(p[1])]; + // MVertex *v3 = _vertices[pointmark(p[2])]; + // MVertex *v4 = _vertices[pointmark(p[3])]; MTetrahedron *t = new MTetrahedron(v1, v2, v3, v4); _gr->tetrahedra.push_back(t); tetloop.tet = tetrahedrontraverse(); -- GitLab