Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2969) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2970) @@ -296,7 +296,7 @@ private: // les arete sont opposes a un sommet - Vertex* ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer - Triangle* at [3]; // 3 adjacent triangles (nu) - Int1 aa[3]; // les nu des arete dans le triangles (mod 3) + Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer + Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles + Int1 TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1 public: Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) @@ -307,5 +307,5 @@ void Echo(); void SetDet() { - if(ns[0] && ns[1] && ns[2]) det = bamg::det(*ns[0],*ns[1],*ns[2]); + if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); else det = -1; } Triangle() {} @@ -313,32 +313,32 @@ Triangle(Vertex *v0,Vertex *v1,Vertex *v2); inline void Set(const Triangle &,const Triangles &,Triangles &); - inline int In(Vertex *v) const { return ns[0]==v || ns[1]==v || ns[2]==v ;} + inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;} TriangleAdjacent FindBoundaryEdge(int i) const; void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){ if (link >=tb && link =tb && at[0] =tb && at[1] =tb && at[2] =tb && TriaAdjTriangles[0] =tb && TriaAdjTriangles[1] =tb && TriaAdjTriangles[2] =vb && ns[0] =vb && ns[1] =vb && ns[2] =vb && TriaVertices[0] =vb && TriaVertices[1] =vb && TriaVertices[2] at[aatt]=this; - tt->aa[aatt]=a + (aa[a] & 60 ) ;}// Copy all the mark + tt->TriaAdjTriangles[aatt]=this; + tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark } void SetAdj2(Int1 a,Triangle *t,Int1 aat){ - //at = pointer toward the adjacent triangles of this - //aa = position of the edge in the triangle (mod 3) - at[a]=t; //the adjacent triangle to the edge a is t - aa[a]=aat; //position of the edge in the adjacent triangle + //TriaAdjTriangles = pointer toward the adjacent triangles of this + //TriaAdjSharedEdge = position of the edge in the triangle (mod 3) + TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t + TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle //if t!=NULL add adjacent triangle to t (this) if(t) { - t->at[aat]=this; - t->aa[aat]=a; + t->TriaAdjTriangles[aat]=this; + t->TriaAdjSharedEdge[aat]=a; } } void SetTriangleContainingTheVertex() { - if (ns[0]) (ns[0]->t=this,ns[0]->vint=0); - if (ns[1]) (ns[1]->t=this,ns[1]->vint=1); - if (ns[2]) (ns[2]->t=this,ns[2]->vint=2); + if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0); + if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1); + if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2); } @@ -375,25 +375,25 @@ Int4 Optim(Int2 a,int =0); - int Locked(int a)const { return aa[a]&4;} - int Hidden(int a)const { return aa[a]&16;} - int Cracked(int a) const { return aa[a] & 32;} + int Locked(int a)const { return TriaAdjSharedEdge[a]&4;} + int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} + int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;} // for optimisation - int GetAllflag(int a){return aa[a] & 1020;} - void SetAllFlag(int a,int f){aa[a] = (aa[a] &3) + (1020 & f);} + int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;} + void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);} void SetHidden(int a){ //Get Adjacent Triangle number a - register Triangle* t = at[a]; + register Triangle* t = TriaAdjTriangles[a]; //if it exist //C|=D -> C=(C|D) bitwise inclusive OR - if(t) t->aa[aa[a] & 3] |=16; - aa[a] |= 16; + if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16; + TriaAdjSharedEdge[a] |= 16; } void SetCracked(int a){ - register Triangle * t = at[a]; - if(t) t->aa[aa[a] & 3] |=32; - aa[a] |= 32; + register Triangle * t = TriaAdjTriangles[a]; + if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32; + TriaAdjSharedEdge[a] |= 32; } @@ -403,19 +403,19 @@ void SetLocked(int a){ //mark the edge as on Boundary - register Triangle * t = at[a]; - t->aa[aa[a] & 3] |=4; - aa[a] |= 4; + register Triangle * t = TriaAdjTriangles[a]; + t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4; + TriaAdjSharedEdge[a] |= 4; } void SetMarkUnSwap(int a){ - register Triangle * t = at[a]; - t->aa[aa[a] & 3] |=8; - aa[a] |=8 ; + register Triangle * t = TriaAdjTriangles[a]; + t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8; + TriaAdjSharedEdge[a] |=8 ; } void SetUnMarkUnSwap(int a){ - register Triangle * t = at[a]; - t->aa[aa[a] & 3] &=55; // 23 + 32 - aa[a] &=55 ; + register Triangle * t = TriaAdjTriangles[a]; + t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 + TriaAdjSharedEdge[a] &=55 ; } @@ -913,14 +913,14 @@ if (link) { int a=-1; - if (aa[0] & 16 ) a=0; - if (aa[1] & 16 ) a=1; - if (aa[2] & 16 ) a=2; + if (TriaAdjSharedEdge[0] & 16 ) a=0; + if (TriaAdjSharedEdge[1] & 16 ) a=1; + if (TriaAdjSharedEdge[2] & 16 ) a=2; if (a>=0) { - t = at[a]; + t = TriaAdjTriangles[a]; // if (t-this<0) return 0; - v2 = ns[VerticesOfTriangularEdge[a][0]]; - v0 = ns[VerticesOfTriangularEdge[a][1]]; - v1 = ns[OppositeEdge[a]]; - v3 = t->ns[OppositeEdge[aa[a]&3]]; + v2 = TriaVertices[VerticesOfTriangularEdge[a][0]]; + v0 = TriaVertices[VerticesOfTriangularEdge[a][1]]; + v1 = TriaVertices[OppositeEdge[a]]; + v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]]; } } @@ -931,18 +931,18 @@ { // first do the logique part double q; - if (!link || aa[a] &4) + if (!link || TriaAdjSharedEdge[a] &4) q= -1; else { - Triangle * t = at[a]; + Triangle * t = TriaAdjTriangles[a]; if (t-this<0) q= -1;// because we do 2 times else if (!t->link ) q= -1; - else if (aa[0] & 16 || aa[1] & 16 || aa[2] & 16 || t->aa[0] & 16 || t->aa[1] & 16 || t->aa[2] & 16 ) + else if (TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16 || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 ) q= -1; else if(option) { - const Vertex & v2 = *ns[VerticesOfTriangularEdge[a][0]]; - const Vertex & v0 = *ns[VerticesOfTriangularEdge[a][1]]; - const Vertex & v1 = *ns[OppositeEdge[a]]; - const Vertex & v3 = * t->ns[OppositeEdge[aa[a]&3]]; + const Vertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]]; + const Vertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]]; + const Vertex & v1 = *TriaVertices[OppositeEdge[a]]; + const Vertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]]; q = QuadQuality(v0,v1,v2,v3); // do the float part } @@ -991,10 +991,10 @@ { *this = rec; - if ( ns[0] ) ns[0] = ThNew.vertices + Th.Number(ns[0]); - if ( ns[1] ) ns[1] = ThNew.vertices + Th.Number(ns[1]); - if ( ns[2] ) ns[2] = ThNew.vertices + Th.Number(ns[2]); - if(at[0]) at[0] = ThNew.triangles + Th.Number(at[0]); - if(at[1]) at[1] = ThNew.triangles + Th.Number(at[1]); - if(at[2]) at[2] = ThNew.triangles + Th.Number(at[2]); + if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices + Th.Number(TriaVertices[0]); + if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices + Th.Number(TriaVertices[1]); + if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices + Th.Number(TriaVertices[2]); + if(TriaAdjTriangles[0]) TriaAdjTriangles[0] = ThNew.triangles + Th.Number(TriaAdjTriangles[0]); + if(TriaAdjTriangles[1]) TriaAdjTriangles[1] = ThNew.triangles + Th.Number(TriaAdjTriangles[1]); + if(TriaAdjTriangles[2]) TriaAdjTriangles[2] = ThNew.triangles + Th.Number(TriaAdjTriangles[2]); if (link >= Th.triangles && link < Th.triangles + Th.nbt) link = ThNew.triangles + Th.Number(link); @@ -1086,10 +1086,10 @@ //set Adjacent Triangle of a triangle if(t) { - t->at[a]=ta.t; - t->aa[a]=ta.a|l; + t->TriaAdjTriangles[a]=ta.t; + t->TriaAdjSharedEdge[a]=ta.a|l; } if(ta.t) { - ta.t->at[ta.a] = t ; - ta.t->aa[ta.a] = a| l ; + ta.t->TriaAdjTriangles[ta.a] = t ; + ta.t->TriaAdjSharedEdge[ta.a] = a| l ; } } @@ -1097,12 +1097,12 @@ inline int TriangleAdjacent::Locked() const - { return t->aa[a] &4;} + { return t->TriaAdjSharedEdge[a] &4;} inline int TriangleAdjacent::Cracked() const - { return t->aa[a] &32;} + { return t->TriaAdjSharedEdge[a] &32;} inline int TriangleAdjacent::GetAllFlag_UnSwap() const - { return t->aa[a] & 1012;} // take all flag except MarkUnSwap + { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap inline int TriangleAdjacent::MarkUnSwap() const - { return t->aa[a] &8;} + { return t->TriaAdjSharedEdge[a] &8;} inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);} @@ -1114,7 +1114,7 @@ inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const - {return t->ns[VerticesOfTriangularEdge[a][i]]; } + {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; } inline Vertex * TriangleAdjacent::OppositeVertex() const - {return t->ns[bamg::OppositeVertex[a]]; } + {return t->TriaVertices[bamg::OppositeVertex[a]]; } inline Icoor2 & TriangleAdjacent::det() const { return t->det;} @@ -1158,18 +1158,18 @@ throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv")); } - ns[0]=v+i; - ns[1]=v+j; - ns[2]=v+k; - at[0]=at[1]=at[2]=0; - aa[0]=aa[1]=aa[2]=0; + TriaVertices[0]=v+i; + TriaVertices[1]=v+j; + TriaVertices[2]=v+k; + TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0; + TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0; det=0; } inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){ - ns[0]=v0; - ns[1]=v1; - ns[2]=v2; - at[0]=at[1]=at[2]=0; - aa[0]=aa[1]=aa[2]=0; + TriaVertices[0]=v0; + TriaVertices[1]=v1; + TriaVertices[2]=v2; + TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0; + TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0; if (v0) det=0; else { @@ -1180,5 +1180,5 @@ inline Real4 Triangle::qualite() { - return det < 0 ? -1 : bamg::qualite(*ns[0],*ns[1],*ns[2]); + return det < 0 ? -1 : bamg::qualite(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); } Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp (revision 2969) +++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp (revision 2970) @@ -43,5 +43,5 @@ k++; //Get ttc, adjacent triangle of t with respect to vertex j - ttc = t->at[j]; + ttc = t->TriaAdjTriangles[j]; //is the current triangle inside or outside? outside = !ttc->link; @@ -52,5 +52,5 @@ t = ttc; //NextEdge[3] = {1,2,0}; - jc = NextEdge[t->aa[j]&3]; + jc = NextEdge[t->TriaAdjSharedEdge[j]&3]; j = NextEdge[jc]; @@ -71,8 +71,8 @@ printf("Triangle:\n"); - printf(" ns pointer towards three vertices\n"); - printf(" ns[0] ns[1] ns[2] = %p %p %p\n",ns[0],ns[1],ns[2]); - printf(" at pointer towards three adjacent triangles\n"); - printf(" at[0] at[1] at[2] = %p %p %p\n",at[0],at[1],at[2]); + printf(" TriaVertices pointer towards three vertices\n"); + printf(" TriaVertices[0] TriaVertices[1] TriaVertices[2] = %p %p %p\n",TriaVertices[0],TriaVertices[1],TriaVertices[2]); + printf(" TriaAdjTriangles pointer towards three adjacent triangles\n"); + printf(" TriaAdjTriangles[0] TriaAdjTriangles[1] TriaAdjTriangles[2] = %p %p %p\n",TriaAdjTriangles[0],TriaAdjTriangles[1],TriaAdjTriangles[2]); printf(" det (integer triangle determinant) = %i\n",det); if (link){ @@ -85,6 +85,6 @@ printf("\nThree vertices:\n"); for(i=0;i<3;i++){ - if (ns[i]){ - ns[i]->Echo(); + if (TriaVertices[i]){ + TriaVertices[i]->Echo(); } else{ @@ -96,4 +96,150 @@ } /*}}}*/ + /*FUNCTION Triangle::Optim{{{1*/ + Int4 Triangle::Optim(Int2 i,int koption) { + // turne around in positif sens + Int4 NbSwap =0; + Triangle *t = this; + int k=0,j =OppositeEdge[i]; + int jp = PreviousEdge[j]; + // initialise tp, jp the previous triangle & edge + Triangle *tp= TriaAdjTriangles[jp]; + jp = TriaAdjSharedEdge[jp]&3; + do { + while (t->swap(j,koption)) + { + NbSwap++; + k++; + if (k>=20000){ + throw ErrorException(__FUNCT__,exprintf("k>=20000")); + } + t= tp->TriaAdjTriangles[jp]; // set unchange t qnd j for previous triangles + j= NextEdge[tp->TriaAdjSharedEdge[jp]&3]; + } + // end on this Triangle + tp = t; + jp = NextEdge[j]; + + t= tp->TriaAdjTriangles[jp]; // set unchange t qnd j for previous triangles + j= NextEdge[tp->TriaAdjSharedEdge[jp]&3]; + + } while( t != this); + return NbSwap; + } + /*}}}1*/ + /*FUNCTION Triangle::swap{{{1*/ + int Triangle::swap(Int2 a,int koption){ + if(a/4 !=0) return 0;// arete lock or MarkUnSwap + + register Triangle *t1=this,*t2=TriaAdjTriangles[a];// les 2 triangles adjacent + register Int1 a1=a,a2=TriaAdjSharedEdge[a];// les 2 numero de l arete dans les 2 triangles + if(a2/4 !=0) return 0; // arete lock or MarkUnSwap + + register Vertex *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]]; + register Vertex *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]]; + register Vertex *s1=t1->TriaVertices[OppositeVertex[a1]]; + register Vertex *s2=t2->TriaVertices[OppositeVertex[a2]]; + + Icoor2 det1=t1->det , det2=t2->det ; + Icoor2 detT = det1+det2; + Icoor2 detA = Abs(det1) + Abs(det2); + Icoor2 detMin = Min(det1,det2); + + int OnSwap = 0; + // si 2 triangle infini (bord) => detT = -2; + if (sa == 0) {// les deux triangles sont frontieres + det2=bamg::det(s2->i,sb->i,s1->i); + OnSwap = det2 >0;} + else if (sb == 0) { // les deux triangles sont frontieres + det1=bamg::det(s1->i,sa->i,s2->i); + OnSwap = det1 >0;} + else if(( s1 != 0) && (s2 != 0) ) { + det1 = bamg::det(s1->i,sa->i,s2->i); + det2 = detT - det1; + OnSwap = (Abs(det1) + Abs(det2)) < detA; + + Icoor2 detMinNew=Min(det1,det2); + // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test + if (! OnSwap &&(detMinNew>0)) { + OnSwap = detMin ==0; + if (! OnSwap) { + int kopt = koption; + while (1) + if(kopt) { + // critere de Delaunay pure isotrope + register Icoor2 xb1 = sb->i.x - s1->i.x, + x21 = s2->i.x - s1->i.x, + yb1 = sb->i.y - s1->i.y, + y21 = s2->i.y - s1->i.y, + xba = sb->i.x - sa->i.x, + x2a = s2->i.x - sa->i.x, + yba = sb->i.y - sa->i.y, + y2a = s2->i.y - sa->i.y; + register double + cosb12 = double(xb1*x21 + yb1*y21), + cosba2 = double(xba*x2a + yba*y2a) , + sinb12 = double(det2), + sinba2 = double(t2->det); + + + // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) + OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); + break; + } + else + { + // critere de Delaunay anisotrope + Real8 som; + I2 AB=(I2) *sb - (I2) *sa; + I2 MAB2=((I2) *sb + (I2) *sa); + R2 MAB(MAB2.x*0.5,MAB2.y*0.5); + I2 A1=(I2) *s1 - (I2) *sa; + I2 D = (I2) * s1 - (I2) * sb ; + R2 S2(s2->i.x,s2->i.y); + R2 S1(s1->i.x,s1->i.y); + { + Metric M=s1->m; + R2 ABo = M.Orthogonal(AB); + R2 A1o = M.Orthogonal(A1); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + if (Abs(d) > dd*1.e-3) { + R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + som = M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + + } + { + Metric M=s2->m; + R2 ABo = M.Orthogonal(AB); + R2 A1o = M.Orthogonal(A1); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + if(Abs(d) > dd*1.e-3) { + R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + som += M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + } + OnSwap = som < 2; + break; + } + + } // OnSwap + } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) + } + if( OnSwap ) + bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); + else { + t1->SetMarkUnSwap(a1); + } + return OnSwap; + } + /*}}}1*/ } Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 2969) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 2970) @@ -3554,8 +3554,5 @@ printf(" NbSwap of insertion: %i\n",NbSwap); printf(" NbSwap/nbv: %i\n",NbSwap/nbv); - printf(" NbUnSwap: %i\n",NbUnSwap); - printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); - } - NbUnSwap = 0; + } #ifdef NBLOOPOPTIM @@ -3572,8 +3569,5 @@ printf(" Optim Loop: %i\n",Nbloop); printf(" NbSwap/nbv: %i\n",NbSwap/nbv); - printf(" NbUnSwap: %i\n",NbUnSwap); - printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); - } - NbUnSwap = 0; + } if(!NbSwap) break; } @@ -3963,36 +3957,4 @@ } /*}}}1*/ -/*FUNCTION Triangles::Optim{{{1*/ -Int4 Triangle::Optim(Int2 i,int koption) { - // turne around in positif sens - Int4 NbSwap =0; - Triangle *t = this; - int k=0,j =OppositeEdge[i]; - int jp = PreviousEdge[j]; - // initialise tp, jp the previous triangle & edge - Triangle *tp= at[jp]; - jp = aa[jp]&3; - do { - while (t->swap(j,koption)) - { - NbSwap++; - k++; - if (k>=20000){ - throw ErrorException(__FUNCT__,exprintf("k>=20000")); - } - t= tp->at[jp]; // set unchange t qnd j for previous triangles - j= NextEdge[tp->aa[jp]&3]; - } - // end on this Triangle - tp = t; - jp = NextEdge[j]; - - t= tp->at[jp]; // set unchange t qnd j for previous triangles - j= NextEdge[tp->aa[jp]&3]; - - } while( t != this); - return NbSwap; -} -/*}}}1*/ /*FUNCTION Triangles::PreInit{{{1*/ void Triangles::PreInit(Int4 inbvx) { @@ -5181,119 +5143,4 @@ } /*}}}1*/ - /*FUNCTION Triangles::swap{{{1*/ - int Triangle::swap(Int2 a,int koption){ - if(a/4 !=0) return 0;// arete lock or MarkUnSwap - - register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent - register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles - if(a2/4 !=0) return 0; // arete lock or MarkUnSwap - - register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]]; - register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]]; - register Vertex *s1=t1->ns[OppositeVertex[a1]]; - register Vertex *s2=t2->ns[OppositeVertex[a2]]; - - Icoor2 det1=t1->det , det2=t2->det ; - Icoor2 detT = det1+det2; - Icoor2 detA = Abs(det1) + Abs(det2); - Icoor2 detMin = Min(det1,det2); - - int OnSwap = 0; - // si 2 triangle infini (bord) => detT = -2; - if (sa == 0) {// les deux triangles sont frontieres - det2=bamg::det(s2->i,sb->i,s1->i); - OnSwap = det2 >0;} - else if (sb == 0) { // les deux triangles sont frontieres - det1=bamg::det(s1->i,sa->i,s2->i); - OnSwap = det1 >0;} - else if(( s1 != 0) && (s2 != 0) ) { - det1 = bamg::det(s1->i,sa->i,s2->i); - det2 = detT - det1; - OnSwap = (Abs(det1) + Abs(det2)) < detA; - - Icoor2 detMinNew=Min(det1,det2); - // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test - if (! OnSwap &&(detMinNew>0)) { - OnSwap = detMin ==0; - if (! OnSwap) { - int kopt = koption; - while (1) - if(kopt) { - // critere de Delaunay pure isotrope - register Icoor2 xb1 = sb->i.x - s1->i.x, - x21 = s2->i.x - s1->i.x, - yb1 = sb->i.y - s1->i.y, - y21 = s2->i.y - s1->i.y, - xba = sb->i.x - sa->i.x, - x2a = s2->i.x - sa->i.x, - yba = sb->i.y - sa->i.y, - y2a = s2->i.y - sa->i.y; - register double - cosb12 = double(xb1*x21 + yb1*y21), - cosba2 = double(xba*x2a + yba*y2a) , - sinb12 = double(det2), - sinba2 = double(t2->det); - - - // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) - OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); - break; - } - else - { - // critere de Delaunay anisotrope - Real8 som; - I2 AB=(I2) *sb - (I2) *sa; - I2 MAB2=((I2) *sb + (I2) *sa); - R2 MAB(MAB2.x*0.5,MAB2.y*0.5); - I2 A1=(I2) *s1 - (I2) *sa; - I2 D = (I2) * s1 - (I2) * sb ; - R2 S2(s2->i.x,s2->i.y); - R2 S1(s1->i.x,s1->i.y); - { - Metric M=s1->m; - R2 ABo = M.Orthogonal(AB); - R2 A1o = M.Orthogonal(A1); - // (A+B)+ x ABo = (S1+B)/2+ y A1 - // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 - double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); - double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 - if (Abs(d) > dd*1.e-3) { - R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); - som = M(C - S2)/M(C - S1); - } else - {kopt=1;continue;} - - } - { - Metric M=s2->m; - R2 ABo = M.Orthogonal(AB); - R2 A1o = M.Orthogonal(A1); - // (A+B)+ x ABo = (S1+B)/2+ y A1 - // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 - double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); - double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 - if(Abs(d) > dd*1.e-3) { - R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); - som += M(C - S2)/M(C - S1); - } else - {kopt=1;continue;} - } - OnSwap = som < 2; - break; - } - - } // OnSwap - } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) - } - if( OnSwap ) - bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); - else { - NbUnSwap ++; - t1->SetMarkUnSwap(a1); - } - return OnSwap; - } - /*}}}1*/ /*FUNCTION Triangles::UnCrack{{{1*/ int Triangles::UnCrack() {