/* T O D O T O D O T O D O T O D O T O D O T O D O T O D O T O D O T O D O Dès que Find_P_in_elt sera en quadtree virer tous les 'last_find' d'accélération relative Dérécursivation de insertPoint et surtout scanNeigh Optimiser la détection d'une frontière dans le coloriage Mettre les lignes de contraintes à NULL quand les Quadtrees seront fonctionnels T O D O T O D O T O D O T O D O T O D O T O D O T O D O T O D O T O D O */ #include #include #include "triangulation.hpp" #define ALTERNATE_INCIRCLE #define DEBUG_FIND //#define STEP_DETAILS #define STEP_WRITE #define STANDALONE #ifdef STANDALONE template inline FLT sqr(const FLT &x) { return x * x; } #else /* STANDALONE */ //#include "timer.hpp" //#include "communicator.hpp" //#include "local_math.hpp" //#include "utils.hpp" #endif /* STANDALONE */ /************************************************************/ /******************** Fonctions locales *********************/ /************************************************************/ //! Produit vectoriel spécialisé 2D inline double SDet(const Vertex & V1, const Vertex & V2, const Vertex & V3) { return (V2[0]-V1[0]) * (V3[1]-V1[1]) - (V2[1]-V1[1]) * (V3[0]-V1[0]); } // inline double SDet(const Vertex & V1, const Vertex & V2, const Vertex & V3) { // Triangle tr(&V1,&V2,&V3); // return Vertex((tr[1]-tr[0]) ^ (tr[2]-tr[0]))[2]; // } inline bool checkInter(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) { // retourne true si P1P2 inter Q1Q2 != vide assert(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2)); // const double // s1 = SDet(P1,P2,Q1)*SDet(P1,P2,Q2), // s2 = SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2); // if (s1 == 0 || s2 == 0) ffout(4) << "warning s=0\n"; // return (s1 <= 0 && s2 <= 0); const double eps = 1e-17; if (std::abs(SDet(P1,P2,Q1)) < eps || std::abs(SDet(P1,P2,Q2)) < eps || std::abs(SDet(Q1,Q2,P1)) < eps || std::abs(SDet(Q1,Q2,P2)) < eps) { //ffout(0) << " !!!!! NO !!!!!!!" << std::endl; return false; } return (SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0); // version optimisée } inline double Angle(const Vertex &V1, const Vertex &V2, const Vertex &V3) { // cos de l'angle * ||V2-V1|| (V1 point pivot) const double distance = sqrt(sqr(V3[0]-V1[0])+sqr(V3[1]-V1[1])); return (((V2[0]-V1[0])*(V3[0]-V1[0]))+((V2[1]-V1[1])*(V3[1]-V1[1])))/distance; } /************************************************************/ /********************** Triangulation ***********************/ /************************************************************/ void Triangulation::setBox(const PointNumType xmin, const PointNumType ymin, const PointNumType xmax, const PointNumType ymax) { triangles.clear(); // pour des bases saines corners[0] = Vertex(xmin,ymax,0); corners[1] = Vertex(xmin,ymin,0); corners[2] = Vertex(xmax,ymin,0); corners[3] = Vertex(xmax,ymax,0); triangles.push_back(ConnectedTriangle(&corners[0],&corners[2],&corners[3])); triangles.push_back(ConnectedTriangle(&corners[0],&corners[1],&corners[2])); triangles.front().neigh(2) = (&triangles.back()); triangles.back().neigh(1) = (&triangles.front()); #ifdef STEP_WRITE { assert(checkAll()); std::ofstream o("init.mesh"); export_mesh(o); } #endif /* STEP_WRITE */ } #ifndef ALTERNATE_INCIRCLE bool Triangulation::isInCircle(const Vertex & position, const Triangle & tri) const { double x1,x2,x3,y1,y2,y3; double a,b,c,d,e,f; double centre_cercle_x,centre_cercle_y; double rayon_cercle,distance; double tolerance = 1.e-13; x1=tri[0][0]; x2=tri[1][0]; x3=tri[2][0]; y1=tri[0][1]; y2=tri[1][1]; y3=tri[2][1]; a=x2-x3; b=y2-y3; c=((y2+y3)*(y2-y3)+(x2+x3)*(x2-x3))/2.; d=x2-x1; e=y2-y1; f=((y2+y1)*(y2-y1)+(x2+x1)*(x2-x1))/2.; centre_cercle_y = (d*c-a*f)/(b*d-a*e); centre_cercle_x = (c*e-f*b)/(a*e-b*d); //printf("centre du cercle :%g %g\n",(d*c-a*f)/(b*d-a*e),(c*e-f*b)/(a*e-b*d)); distance = sqr(centre_cercle_x-position[0]) + sqr(centre_cercle_y-position[1]); rayon_cercle = sqr(centre_cercle_x-x1) + sqr(centre_cercle_y-y1); return (distance<=(rayon_cercle+tolerance)); } #else /* ALTERNATE_INCIRCLE */ bool Triangulation::isInCircle(const Vertex & position, const Triangle & tri) const { const double eps = 1e-13; // On peut éviter quelques indirections par des variables locales /* 9 affectations : 3+3+3*3 (+) et 0+0+3*2 (*) */ const double BxmAx=tri(1)[0]-tri(0)[0], CxmAx=tri(2)[0]-tri(0)[0], DxmAx=position[0]-tri(0)[0]; const double BymAy=tri(1)[1]-tri(0)[1], CymAy=tri(2)[1]-tri(0)[1], DymAy=position[1]-tri(0)[1]; const double BA2=BxmAx*(tri(1)[0]+tri(0)[0])+BymAy*(tri(1)[1]+tri(0)[1]), CA2=CxmAx*(tri(2)[0]+tri(0)[0])+CymAy*(tri(2)[1]+tri(0)[1]), DA2=DxmAx*(position[0]+tri(0)[0])+DymAy*(position[1]+tri(0)[1]); const double r= // 5 (+) , 3*3 (*) (BxmAx*CymAy-CxmAx*BymAy)*DA2+(CxmAx*DymAy-DxmAx*CymAy)*BA2+(DxmAx*BymAy-BxmAx*DymAy)*CA2; // Total 20 (+) et 12 (*) return (r modification de V1 qui deivient T2 et V3 qui devient T1 const TriangleIndex num_tri2 = num_tri1->neigh(k); assert(num_tri2 != NULL); PointIndex S1=num_tri1->base()[k], S2=num_tri1->base()[(k+1)%3], S4=num_tri1->base()[(k+2)%3]; TriangleIndex V1=num_tri1->neigh((k+1)%3), V2=num_tri1->neigh((k+2)%3); PointIndex S3; TriangleIndex V3,V4; // recherche du sommet dans le triangle num_tri2 opposé au coté commun (ie opposé à k dans num_tri1) for(unsigned k=0;k<3;++k) if(num_tri2->neigh(k) == num_tri1) { S3=num_tri2->base()[k]; V3=num_tri2->neigh((k+1)%3); // Ici, un bug est mort dans d'attroces souffrances V4=num_tri2->neigh((k+2)%3); break; } for(unsigned k=0;k<3;++k) { if(V1 != NULL) if(V1->neigh(k) == num_tri1) V1->neigh(k) = num_tri2; if(V3 != NULL) if(V3->neigh(k) == num_tri2) V3->neigh(k) = num_tri1; } num_tri1->setVertices(S1,S2,S3); num_tri2->setVertices(S1,S3,S4); num_tri1->setNeighs(V3,num_tri2,V2); num_tri2->setNeighs(V4,V1,num_tri1); } void Triangulation::elimination_tri(TriangleIndex num_tri) { // elimination fu triangle pointe par neigh_index for(unsigned k=0;k<3;++k) { TriangleIndex neigh = num_tri->neigh(k); if(neigh != NULL) { for(unsigned k2=0;k2<3;++k2) if(neigh->neigh(k2) == num_tri) { neigh->neigh(k2) = NULL; break; } } num_tri->neigh(k) = NULL; } // num_tri->tag = -1; deleteTriangle(num_tri); } // recherche du triangle comportant l'element P TriangleIndex Triangulation::find_P_in_elt(const Vertex & P) { #ifdef DEBUG_FIND // Version non optimisée de debug (en cas de trous de maillages supposé convexe) for(Triangles::iterator i = triangles.begin();true;++i) { assert(i != triangles.end()); const Vertex & P1 = (*i)(0); const Vertex & P2 = (*i)(1); const Vertex & P3 = (*i)(2); assert(SDet(P1,P2,P3)>=0); // ffout(4)<=-1e-15) && (SDet(P2,P3,P)>=-1e-15) && (SDet(P3,P1,P)>=-1e-15) ) { return &*i; } } #else /* DEBUG_FIND */ // static TriangleIndex lastFind = NULL; // if (lastFind == NULL) lastFind = &*triangles.begin(); for(TriangleIndex i = &*triangles.begin();;) { assert(i != NULL); const Vertex & P2 = (*i)(1); const Vertex & P3 = (*i)(2); if(SDet(P2,P3,P)<0) { i=i->neigh(0); } else { const Vertex & P1 = (*i)(0); // on en a pas besoin avant if(SDet(P3,P1,P)<0) { i=i->neigh(1); } else { if(SDet(P1,P2,P)<0) { i=i->neigh(2); } else { // lastFind = i; // Lié à Find_P_in_elt return i; } } } } #endif /* DEBUG_FIND */ } void Triangulation::insertPoint(Vertex * newPoint, const TriangleIndex & curTriangle) { // triangle forcant cette position; évite les problèmes dus à l'éffacement de triangles newTriangle(ConnectedTriangle(&corners[0],&corners[0],&corners[0])); const Triangles::iterator first = --triangles.end(); Vertex * last = curTriangle->base()[1]; for(unsigned i=0;i<3;++i) { if (scanNeigh(curTriangle, i, newPoint, last) ) { // au bord de la bulle TriangleIndex tmpNeigh = curTriangle->neigh(i); newTriangle(ConnectedTriangle(last,curTriangle->base()[(i+2)%3],newPoint,NULL,NULL,tmpNeigh)); if (tmpNeigh != NULL) tmpNeigh->neigh(tmpNeigh->localizeNeigh(curTriangle)) = &triangles.back(); last = curTriangle->base()[(i+2)%3]; } } deleteTriangle(curTriangle); Triangles::iterator LTi = first, LTf = triangles.end(), LTaux; ++LTi; LTi->neigh(1) = &*(--LTf); // Cas particulier LTf->neigh(0) = &*LTi; LTaux = LTi; ++LTaux; while(LTi != LTf) { // ici LTf = --triangle.end() LTi->neigh(0) = &*LTaux; LTaux->neigh(1) = &*LTi; ++LTi; ++LTaux; } triangles.erase(first); // plus utile } bool Triangulation::scanNeigh(const TriangleIndex & start, const unsigned in, Vertex * newPoint, Vertex * &last) { // Scan à partir du triangle Start par le voisin in, à la recherche de triangles voisins contenant dans // leur cercle circonscrit le point newPoint; si OK met à jour la liste L entre LP1 et LP2 const TriangleIndex curNeigh = start->neigh(in); if (curNeigh != NULL && isInCircle(*newPoint,*curNeigh)) { const unsigned k = curNeigh->localizeNeigh(&*start); // sommet opposé au coté d'entrée if (scanNeigh(curNeigh, (k+1)%3, newPoint, last) ) { // au bord de la bulle TriangleIndex tmpNeigh = curNeigh->neigh((k+1)%3); newTriangle(ConnectedTriangle(last,curNeigh->base()[k],newPoint,NULL,NULL,tmpNeigh)); if (tmpNeigh != NULL) tmpNeigh->neigh(tmpNeigh->localizeNeigh(curNeigh)) = &triangles.back(); last = curNeigh->base()[k]; } // Idem que précédemment, mais pour pas de boucle, code dupliqué if (scanNeigh(curNeigh, (k+2)%3, newPoint, last)) { // au bord de la bulle TriangleIndex tmpNeigh = curNeigh->neigh((k+2)%3); newTriangle(ConnectedTriangle(last,curNeigh->base()[(k+1)%3],newPoint,NULL,NULL,tmpNeigh)); if (tmpNeigh != NULL) tmpNeigh->neigh(tmpNeigh->localizeNeigh(curNeigh)) = &triangles.back(); last = curNeigh->base()[(k+1)%3]; } deleteTriangle(curNeigh); return false; } return true; } void Triangulation::resetTags() { for(Triangles::iterator i = triangles.begin(); i != triangles.end(); ++i) i->tag = 0; } bool Triangulation::checkLine(const CurveVertex &L, const bool closed) { assert(L.size() > 1); if (closed) assert(L.front() == L.back()); bool change = false; CurveVertex::const_iterator position = L.begin(); const CurveVertex::const_iterator end = L.end(); PointIndex curVertex = *position++, curTarget = *position; // Etape initiale : recherche d'un triangle possédant en sommet le premier point de L // On peut y mettre une recherche quadtree ou par adjacente optimale (cf find_P_in_elt) TriangleIndex curTriangle = findVertex(curVertex); bool found = false; while (!found) { assert(curTriangle->base()[0] == curVertex || curTriangle->base()[1] == curVertex || curTriangle->base()[2] == curVertex); unsigned localVertex=curTriangle->localizeVertex(curVertex); // position local de curVertex if (curTriangle->base()[(localVertex+1)%3] == curTarget) { // On a trouvé curTarget dans le triangle courant, on passe au suivant. // On marque le coté comme frontière si closed == true ie c'est pas une ligne interne if (closed) { // coté entre localVertex et LocalVertex+1 => neigh[localVertex+2] curTriangle->tag |= HSide[(localVertex+2)%3]; TriangleIndex NTri = curTriangle->neigh((localVertex+2)%3); // forcement non null assert(NTri != NULL); // Car tout est inclus dans la boite NTri->tag |= HSide[NTri->localizeNeigh(curTriangle)]; } curVertex = curTarget; found = (++position == end); curTarget = *position; } else if (curTriangle->base()[(localVertex+2)%3] == curTarget) { // On a trouvé curTarget dans le triangle courant, on passe au suivant. // On marque le coté comme frontière si closed == true ie c'est pas une ligne interne if (closed) { // coté entre localVertex et LocalVertex+1 => neigh[localVertex+2] curTriangle->tag |= HSide[(localVertex+1)%3]; TriangleIndex NTri = curTriangle->neigh((localVertex+1)%3); // forcement non null assert(NTri != NULL); // Car tout est inclus dans la boite NTri->tag |= HSide[NTri->localizeNeigh(curTriangle)]; } curVertex = curTarget; found = (++position == end); curTarget = *position; } else if(checkInter(*curVertex,*curTarget, (*curTriangle)((localVertex+1)%3), (*curTriangle)((localVertex+2)%3))) { // Ok : found a problem on the border // gestion spécifique du cnflit arète-frontière bool atteint = false; while(!atteint) { assert(curTriangle != NULL); // description globale du coté opposée à localvertex dans curTriangle ie l'entrée const PointIndex g1 = curTriangle->base()[(localVertex+1)%3], g2 = curTriangle->base()[(localVertex+2)%3]; TriangleIndex new_t = curTriangle->neigh(localVertex); assert(new_t != NULL); unsigned l=new_t->localizeNeigh(curTriangle); const PointIndex from = curTriangle->base()[localVertex]; const PointIndex to = new_t->base()[l]; // Test de convexité if (checkInter(*from,*to,*g1,*g2)) { // c'est convexe ou fusinnable #ifndef STEP_DETAILS // Test de fusion // g1 sur [from,to] if (SDet(*from,*to,*g1) == 0) ffout(4) << "G1 is on the border\n"; // g2 sur [from,to] if (SDet(*from,*to,*g2) == 0) ffout(4) << "G2 is on the border\n"; // from sur [g1,g2] if (SDet(*g1,*g2,*from) == 0) ffout(4) << "from is on the border\n"; // tp sur [g1,g2] if (SDet(*g1,*g2,*to) == 0) ffout(4) << "to is on the border\n"; #endif /* STEP_DETAILS */ permutation(curTriangle,localVertex); change = true; // changement // le cas où g3 serait curTarget localVertex=curTriangle->localizeVertex(from); l=new_t->localizeVertex(from); if (to == curTarget) { // passer au point suivant. #warning ! // if (closed && from == curVertex) { // // on a l'arête complète // curTriangle->tag |= HSide[curTriangle->localizeNeigh(new_t)]; // new_t->tag |= HSide[new_t->localizeNeigh(curTriangle)]; // } #warning end of ! curVertex = curTarget; found = (++position == end); curTarget = *position; atteint = true; } else { // on pourrait y mettre les noms globaux que l'on y connait if (checkInter(*curVertex,*curTarget, (*curTriangle)((localVertex+1)%3), (*curTriangle)((localVertex+2)%3))) { // curTriangle : OK // localVertex : OK } else { curTriangle = new_t; localVertex = l; assert(checkInter(*curVertex,*curTarget, (*curTriangle)((l+1)%3), (*curTriangle)((l+2)%3))); } } } else { // c'est pas convexe et ça coupe encore l'arète curTriangle = new_t; // test de continuation if (checkInter(*curVertex,*curTarget, (*curTriangle)(l), (*curTriangle)((l+1)%3))) { // sinon c'est l'autre // c'est sur la droite [l, l+1] opp à l+2 localVertex = (l+2)%3; } else { // c'est sur la droite [l,l+2] opp à l+1 assert(checkInter(*curVertex,*curTarget, (*curTriangle)(l), (*curTriangle)((l+2)%3))); localVertex = (l+1)%3; } } } } else { // A fortment optimiser, peut coupable avec le calcul d'intersection const double theta1 = Angle(*curVertex,*curTarget, (*curTriangle)((localVertex+1)%3)); const double theta2 = Angle(*curVertex,*curTarget, (*curTriangle)((localVertex+2)%3)); const int sens = (theta1 > theta2)?2:1; // ffout(4) << "angles : " << theta1 << " " << theta2 << " sens : " << sens << "\n"; curTriangle = curTriangle->neigh((localVertex+sens)%3); } } // end of while(!found) return change; } bool Triangulation::colorize(const Points & points) { // les variables pour les différents parcours std::stack Stack; // On parcours les triangles pour propager la couleur de l'extérieur : 1 // Recherche de l'élément initial : ie un triangle de coin de la box englobante Stack.push(find_P_in_elt(corners[0])); while(!Stack.empty()) { TriangleIndex itr = Stack.top(); Stack.pop(); if (!(itr->tag & HMask)) { // pas de couleur itr->tag |= 1; // 1 : Couleur de remplissage // parcours des voisins (il faudrait peut être pas de boucle pour être plus rapide) for(unsigned j=0;j<3;j++) { TriangleIndex ktr = itr->neigh(j); if (ktr != NULL && // hors domaine !(itr->tag & HSide[j]) && !(ktr->tag & HMask) // Déjà colorié ) { Stack.push(ktr); } } } } Triangles::iterator llf = triangles.begin(); const Triangles::iterator lle = triangles.end(); while(llf != lle && !((llf->tag & HUMask) && !(llf->tag & HMask) && (llf->neigh(0)->tag & HMask || llf->neigh(1)->tag & HMask || llf->neigh(2)->tag & HMask))) { ++llf; } if (llf == lle) { #ifndef STEP_DETAILS ffout(4) << "Bad Domain\n"; // c'est pour cela qu'il faut commencer par colorier l'extérieur #endif /* STEP_DETAILS */ return false; } // On parcours les triangles pour propager une couleur // Initialisation Stack.push(&*llf); while(!Stack.empty()) { TriangleIndex itr = Stack.top(); Stack.pop(); if (!(itr->tag & HMask)) { // pas de couleur itr->tag |= 2; // colorisation du domaine avec la couleur 2 // parcours des voisins (il faudrait peut être pas de boucle pour être plus rapide) for(unsigned j=0;j<3;j++) { TriangleIndex ktr = itr->neigh(j); if (ktr != NULL && !(itr->tag & HSide[j]) && !(ktr->tag & HMask)) { Stack.push(ktr); } } } } return true; } TriangleIndex Triangulation::findVertex(const Vertex * const p) { for(Triangles::iterator i = triangles.begin(); i != triangles.end(); ++i) { for(unsigned j=0;j<3;++j) if (i->base()[j] == p) return &*i; } assert(false); return NULL; } bool Triangulation::triangulize(Points & points, const CurveVertex & envVertex, const std::vector & holes, const std::vector & curves) { assert(points.size() > 0); _points_ = &points; const Vertex & firstPoint = points[0]; PointNumType xmin=firstPoint[0],xmax=firstPoint[0], ymin=firstPoint[1],ymax=firstPoint[1]; for(unsigned i=0;i::iterator lf = holes.begin(); lf != holes.end(); ++lf) // for(unsigned i=1;checkLine(*lf,true);) // ffout(4) << "\tPasse " << ++i << "\n"; //#else /* STEP_DETAILS */ for(std::vector::const_iterator lf = holes.begin(); lf != holes.end(); ++lf) while(checkLine(*lf,true)); //#endif /* STEP_DETAILS */ #ifdef STEP_DETAILS ffout(4) << "Check Courbe(s)" << "\n"; for(std::vector::const_iterator lf = curves.begin(); lf != curves.end(); ++lf) for(unsigned i=1;checkLine(*lf,false);) ffout(4) << "\tPasse " << ++i << "\n"; #else /* STEP_DETAILS */ for(std::vector::const_iterator lf = curves.begin(); lf != curves.end(); ++lf) while(checkLine(*lf,false)); #endif /* STEP_DETAILS */ #ifdef STEP_WRITE { assert(checkAll()); std::ofstream o("check.mesh"); export_mesh(o); } #endif /* STEP_WRITE */ #ifdef STEP_DETAILS ffout(4) << "Colorize" << "\n"; #endif /* STEP_DETAILS */ colorize(points); #ifdef STEP_WRITE { assert(checkAll()); std::ofstream o("color.mesh"); export_mesh(o); } #endif /* STEP_WRITE */ unsigned count_tri = 0; const Triangles::iterator e = triangles.end(); for(Triangles::iterator ff,f = triangles.begin(); f != e;) { ff = f++; if ((ff->tag & HMask) != 2) { elimination_tri(&*ff); ++count_tri; } } #ifdef STEP_DETAILS ffout(4) << count_tri << " Triangles Deleted " << triangles.size() << " triangles left" << "\n"; #endif /* STEP_DETAILS */ #ifdef STEP_WRITE { assert(checkAll()); std::ofstream o("final.mesh"); export_mesh(o); } #endif /* STEP_WRITE */ return true; } void Triangulation::export_mesh(std::ostream &o) const { const Points & points = *_points_; o << "MeshVersionFormatted 1\n" "Dimension\n" "\t2\n" "Vertices\n" "\t" << points.size()+4 << "\n"; for(unsigned i=0;i<4;++i) o << corners[i][0] << " " << corners[i][1] << " 1\n"; for(Points::const_iterator i = points.begin(); i != points.end(); ++i) { o << i->x() << " " << i->y() << " 1\n"; } class BestRef { private: public: const Vertex *points, *corners; public: // BestRef(const Points * _points, const Point * _corners) // : points(_points), // corners(_corners) { } unsigned operator()(const Vertex & p) { const unsigned n = &p-corners; if (n < 4) return n+1; else return &p-points+5; } } ref; // (firstPointer(points),firstPointer(corners)); ref.points = &points[0]; ref.corners = corners; o << "Triangles\n\t" << triangles.size() << "\n"; for(Triangles::const_iterator i = triangles.begin(); i != triangles.end(); ++i) { o << " " << ref((*i)(0)) << " " << ref((*i)(1)) << " " << ref((*i)(2)) << " " << (i->tag & HMask) << "\n"; } } bool Triangulation::checkAll() { bool ok = true; unsigned jj = 0; for(Triangles::iterator i = triangles.begin(); i != triangles.end(); ++i,++jj) { // Check Orientation const double S = SDet((*i)(0),(*i)(1),(*i)(2)); const bool check_point = ((i->base())[0] != (i->base())[1]) && ((i->base())[0] != (i->base())[2]) && ((i->base())[1] != (i->base())[2]); if (!check_point) { ffout(4) << "Warning : Triangle " << jj << " Plat : " /*<< *i */<< "\n"; ok = false; } assert(check_point); if (S <= 0) { ffout(4) << "< triangle " << jj << " " /*<< *i */<< " : " << S << "\n"; ok = false; } // Check Neigh for(unsigned j=0;j<3;j++) { const TriangleIndex k = i->neigh(j); if (k != NULL) { const unsigned l = k->localizeNeigh(&*i); if (!(k->base()[(l+2)%3] == i->base()[(j+1)%3] && k->base()[(l+1)%3] == i->base()[(j+2)%3])) { ffout(4) << "Bad edge " << j << " for triangle : " /*<< *i */<< "\n"; ok = false; } } } } return ok; }