![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2970
- Timestamp:
- 02/05/10 14:11:56 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 3 edited
-
Mesh2.h (modified) (14 diffs)
-
objects/Triangle.cpp (modified) (5 diffs)
-
objects/Triangles.cpp (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2969 r2970 296 296 297 297 private: // les arete sont opposes a un sommet 298 Vertex* ns [3];// 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer299 Triangle* at [3]; // 3 adjacent triangles (nu)300 Int1 aa[3]; // les nu des arete dans le triangles (mod 3)298 Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer 299 Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles 300 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 301 301 public: 302 302 Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) … … 307 307 void Echo(); 308 308 void SetDet() { 309 if( ns[0] && ns[1] && ns[2]) det = bamg::det(*ns[0],*ns[1],*ns[2]);309 if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); 310 310 else det = -1; } 311 311 Triangle() {} … … 313 313 Triangle(Vertex *v0,Vertex *v1,Vertex *v2); 314 314 inline void Set(const Triangle &,const Triangles &,Triangles &); 315 inline int In(Vertex *v) const { return ns[0]==v || ns[1]==v || ns[2]==v ;}315 inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;} 316 316 TriangleAdjacent FindBoundaryEdge(int i) const; 317 317 318 318 void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){ 319 319 if (link >=tb && link <te) link = tb + renu[link -tb]; 320 if ( at[0] >=tb && at[0] <te) at[0] = tb + renu[at[0]-tb];321 if ( at[1] >=tb && at[1] <te) at[1] = tb + renu[at[1]-tb];322 if ( at[2] >=tb && at[2] <te) at[2] = tb + renu[at[2]-tb];320 if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb]; 321 if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb]; 322 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb]; 323 323 } 324 324 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){ 325 if ( ns[0] >=vb && ns[0] <ve) ns[0] = vb + renu[ns[0]-vb];326 if ( ns[1] >=vb && ns[1] <ve) ns[1] = vb + renu[ns[1]-vb];327 if ( ns[2] >=vb && ns[2] <ve) ns[2] = vb + renu[ns[2]-vb];328 } 329 330 const Vertex & operator[](int i) const {return * ns[i];};331 Vertex & operator[](int i) {return * ns[i];};332 333 const Vertex * operator()(int i) const {return ns[i];};334 Vertex * & operator()(int i) {return ns[i];};325 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb]; 326 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb]; 327 if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb]; 328 } 329 330 const Vertex & operator[](int i) const {return *TriaVertices[i];}; 331 Vertex & operator[](int i) {return *TriaVertices[i];}; 332 333 const Vertex * operator()(int i) const {return TriaVertices[i];}; 334 Vertex * & operator()(int i) {return TriaVertices[i];}; 335 335 336 336 TriangleAdjacent Adj(int i) const // triangle adjacent + arete 337 { return TriangleAdjacent( at[i],aa[i]&3);};337 { return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);}; 338 338 339 339 Triangle * TriangleAdj(int i) const 340 {return at[i&3];} // triangle adjacent + arete340 {return TriaAdjTriangles[i&3];} // triangle adjacent + arete 341 341 Int1 NuEdgeTriangleAdj(int i) const 342 {return aa[i&3]&3;} // Number of the adjacent edge in adj tria342 {return TriaAdjSharedEdge[i&3]&3;} // Number of the adjacent edge in adj tria 343 343 344 344 inline Real4 qualite() ; … … 346 346 void SetAdjAdj(Int1 a) 347 347 { a &= 3; 348 register Triangle *tt= at[a];349 aa[a] &= 55; // remove MarkUnSwap350 register Int1 aatt = aa[a] & 3;348 register Triangle *tt=TriaAdjTriangles[a]; 349 TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap 350 register Int1 aatt = TriaAdjSharedEdge[a] & 3; 351 351 if(tt){ 352 tt-> at[aatt]=this;353 tt-> aa[aatt]=a + (aa[a] & 60 ) ;}// Copy all the mark352 tt->TriaAdjTriangles[aatt]=this; 353 tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark 354 354 } 355 355 356 356 void SetAdj2(Int1 a,Triangle *t,Int1 aat){ 357 // at= pointer toward the adjacent triangles of this358 // aa= position of the edge in the triangle (mod 3)359 at[a]=t; //the adjacent triangle to the edge a is t360 aa[a]=aat; //position of the edge in the adjacent triangle357 //TriaAdjTriangles = pointer toward the adjacent triangles of this 358 //TriaAdjSharedEdge = position of the edge in the triangle (mod 3) 359 TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t 360 TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle 361 361 //if t!=NULL add adjacent triangle to t (this) 362 362 if(t) { 363 t-> at[aat]=this;364 t-> aa[aat]=a;363 t->TriaAdjTriangles[aat]=this; 364 t->TriaAdjSharedEdge[aat]=a; 365 365 } 366 366 } 367 367 368 368 void SetTriangleContainingTheVertex() { 369 if ( ns[0]) (ns[0]->t=this,ns[0]->vint=0);370 if ( ns[1]) (ns[1]->t=this,ns[1]->vint=1);371 if ( ns[2]) (ns[2]->t=this,ns[2]->vint=2);369 if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0); 370 if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1); 371 if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2); 372 372 } 373 373 … … 375 375 Int4 Optim(Int2 a,int =0); 376 376 377 int Locked(int a)const { return aa[a]&4;}378 int Hidden(int a)const { return aa[a]&16;}379 int Cracked(int a) const { return aa[a] & 32;}377 int Locked(int a)const { return TriaAdjSharedEdge[a]&4;} 378 int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} 379 int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;} 380 380 // for optimisation 381 int GetAllflag(int a){return aa[a] & 1020;}382 void SetAllFlag(int a,int f){ aa[a] = (aa[a] &3) + (1020 & f);}381 int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;} 382 void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);} 383 383 384 384 void SetHidden(int a){ 385 385 386 386 //Get Adjacent Triangle number a 387 register Triangle* t = at[a];387 register Triangle* t = TriaAdjTriangles[a]; 388 388 389 389 //if it exist 390 390 //C|=D -> C=(C|D) bitwise inclusive OR 391 if(t) t-> aa[aa[a] & 3] |=16;392 aa[a] |= 16;391 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16; 392 TriaAdjSharedEdge[a] |= 16; 393 393 } 394 394 void SetCracked(int a){ 395 register Triangle * t = at[a];396 if(t) t-> aa[aa[a] & 3] |=32;397 aa[a] |= 32;395 register Triangle * t = TriaAdjTriangles[a]; 396 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32; 397 TriaAdjSharedEdge[a] |= 32; 398 398 } 399 399 … … 403 403 void SetLocked(int a){ 404 404 //mark the edge as on Boundary 405 register Triangle * t = at[a];406 t-> aa[aa[a] & 3] |=4;407 aa[a] |= 4;405 register Triangle * t = TriaAdjTriangles[a]; 406 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4; 407 TriaAdjSharedEdge[a] |= 4; 408 408 } 409 409 410 410 void SetMarkUnSwap(int a){ 411 register Triangle * t = at[a];412 t-> aa[aa[a] & 3] |=8;413 aa[a] |=8 ;411 register Triangle * t = TriaAdjTriangles[a]; 412 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8; 413 TriaAdjSharedEdge[a] |=8 ; 414 414 } 415 415 416 416 void SetUnMarkUnSwap(int a){ 417 register Triangle * t = at[a];418 t-> aa[aa[a] & 3] &=55; // 23 + 32419 aa[a] &=55 ;417 register Triangle * t = TriaAdjTriangles[a]; 418 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 419 TriaAdjSharedEdge[a] &=55 ; 420 420 } 421 421 … … 913 913 if (link) { 914 914 int a=-1; 915 if ( aa[0] & 16 ) a=0;916 if ( aa[1] & 16 ) a=1;917 if ( aa[2] & 16 ) a=2;915 if (TriaAdjSharedEdge[0] & 16 ) a=0; 916 if (TriaAdjSharedEdge[1] & 16 ) a=1; 917 if (TriaAdjSharedEdge[2] & 16 ) a=2; 918 918 if (a>=0) { 919 t = at[a];919 t = TriaAdjTriangles[a]; 920 920 // if (t-this<0) return 0; 921 v2 = ns[VerticesOfTriangularEdge[a][0]];922 v0 = ns[VerticesOfTriangularEdge[a][1]];923 v1 = ns[OppositeEdge[a]];924 v3 = t-> ns[OppositeEdge[aa[a]&3]];921 v2 = TriaVertices[VerticesOfTriangularEdge[a][0]]; 922 v0 = TriaVertices[VerticesOfTriangularEdge[a][1]]; 923 v1 = TriaVertices[OppositeEdge[a]]; 924 v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]]; 925 925 } 926 926 } … … 931 931 { // first do the logique part 932 932 double q; 933 if (!link || aa[a] &4)933 if (!link || TriaAdjSharedEdge[a] &4) 934 934 q= -1; 935 935 else { 936 Triangle * t = at[a];936 Triangle * t = TriaAdjTriangles[a]; 937 937 if (t-this<0) q= -1;// because we do 2 times 938 938 else if (!t->link ) q= -1; 939 else if ( aa[0] & 16 || aa[1] & 16 || aa[2] & 16 || t->aa[0] & 16 || t->aa[1] & 16 || t->aa[2] & 16 )939 else if (TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16 || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 ) 940 940 q= -1; 941 941 else if(option) 942 942 { 943 const Vertex & v2 = * ns[VerticesOfTriangularEdge[a][0]];944 const Vertex & v0 = * ns[VerticesOfTriangularEdge[a][1]];945 const Vertex & v1 = * ns[OppositeEdge[a]];946 const Vertex & v3 = * t-> ns[OppositeEdge[aa[a]&3]];943 const Vertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]]; 944 const Vertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]]; 945 const Vertex & v1 = *TriaVertices[OppositeEdge[a]]; 946 const Vertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]]; 947 947 q = QuadQuality(v0,v1,v2,v3); // do the float part 948 948 } … … 991 991 { 992 992 *this = rec; 993 if ( ns[0] ) ns[0] = ThNew.vertices + Th.Number(ns[0]);994 if ( ns[1] ) ns[1] = ThNew.vertices + Th.Number(ns[1]);995 if ( ns[2] ) ns[2] = ThNew.vertices + Th.Number(ns[2]);996 if( at[0]) at[0] = ThNew.triangles + Th.Number(at[0]);997 if( at[1]) at[1] = ThNew.triangles + Th.Number(at[1]);998 if( at[2]) at[2] = ThNew.triangles + Th.Number(at[2]);993 if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices + Th.Number(TriaVertices[0]); 994 if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices + Th.Number(TriaVertices[1]); 995 if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices + Th.Number(TriaVertices[2]); 996 if(TriaAdjTriangles[0]) TriaAdjTriangles[0] = ThNew.triangles + Th.Number(TriaAdjTriangles[0]); 997 if(TriaAdjTriangles[1]) TriaAdjTriangles[1] = ThNew.triangles + Th.Number(TriaAdjTriangles[1]); 998 if(TriaAdjTriangles[2]) TriaAdjTriangles[2] = ThNew.triangles + Th.Number(TriaAdjTriangles[2]); 999 999 if (link >= Th.triangles && link < Th.triangles + Th.nbt) 1000 1000 link = ThNew.triangles + Th.Number(link); … … 1086 1086 //set Adjacent Triangle of a triangle 1087 1087 if(t) { 1088 t-> at[a]=ta.t;1089 t-> aa[a]=ta.a|l;1088 t->TriaAdjTriangles[a]=ta.t; 1089 t->TriaAdjSharedEdge[a]=ta.a|l; 1090 1090 } 1091 1091 if(ta.t) { 1092 ta.t-> at[ta.a] = t ;1093 ta.t-> aa[ta.a] = a| l ;1092 ta.t->TriaAdjTriangles[ta.a] = t ; 1093 ta.t->TriaAdjSharedEdge[ta.a] = a| l ; 1094 1094 } 1095 1095 } … … 1097 1097 1098 1098 inline int TriangleAdjacent::Locked() const 1099 { return t-> aa[a] &4;}1099 { return t->TriaAdjSharedEdge[a] &4;} 1100 1100 inline int TriangleAdjacent::Cracked() const 1101 { return t-> aa[a] &32;}1101 { return t->TriaAdjSharedEdge[a] &32;} 1102 1102 inline int TriangleAdjacent::GetAllFlag_UnSwap() const 1103 { return t-> aa[a] & 1012;} // take all flag except MarkUnSwap1103 { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap 1104 1104 1105 1105 inline int TriangleAdjacent::MarkUnSwap() const 1106 { return t-> aa[a] &8;}1106 { return t->TriaAdjSharedEdge[a] &8;} 1107 1107 1108 1108 inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);} … … 1114 1114 1115 1115 inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const 1116 {return t-> ns[VerticesOfTriangularEdge[a][i]]; }1116 {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; } 1117 1117 inline Vertex * TriangleAdjacent::OppositeVertex() const 1118 {return t-> ns[bamg::OppositeVertex[a]]; }1118 {return t->TriaVertices[bamg::OppositeVertex[a]]; } 1119 1119 inline Icoor2 & TriangleAdjacent::det() const 1120 1120 { return t->det;} … … 1158 1158 throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv")); 1159 1159 } 1160 ns[0]=v+i;1161 ns[1]=v+j;1162 ns[2]=v+k;1163 at[0]=at[1]=at[2]=0;1164 aa[0]=aa[1]=aa[2]=0;1160 TriaVertices[0]=v+i; 1161 TriaVertices[1]=v+j; 1162 TriaVertices[2]=v+k; 1163 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0; 1164 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0; 1165 1165 det=0; 1166 1166 } 1167 1167 1168 1168 inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){ 1169 ns[0]=v0;1170 ns[1]=v1;1171 ns[2]=v2;1172 at[0]=at[1]=at[2]=0;1173 aa[0]=aa[1]=aa[2]=0;1169 TriaVertices[0]=v0; 1170 TriaVertices[1]=v1; 1171 TriaVertices[2]=v2; 1172 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0; 1173 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0; 1174 1174 if (v0) det=0; 1175 1175 else { … … 1180 1180 inline Real4 Triangle::qualite() 1181 1181 { 1182 return det < 0 ? -1 : bamg::qualite(* ns[0],*ns[1],*ns[2]);1182 return det < 0 ? -1 : bamg::qualite(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); 1183 1183 } 1184 1184 -
issm/trunk/src/c/Bamgx/objects/Triangle.cpp
r2969 r2970 43 43 k++; 44 44 //Get ttc, adjacent triangle of t with respect to vertex j 45 ttc = t-> at[j];45 ttc = t->TriaAdjTriangles[j]; 46 46 //is the current triangle inside or outside? 47 47 outside = !ttc->link; … … 52 52 t = ttc; 53 53 //NextEdge[3] = {1,2,0}; 54 jc = NextEdge[t-> aa[j]&3];54 jc = NextEdge[t->TriaAdjSharedEdge[j]&3]; 55 55 j = NextEdge[jc]; 56 56 … … 71 71 72 72 printf("Triangle:\n"); 73 printf(" ns pointer towards three vertices\n");74 printf(" ns[0] ns[1] ns[2] = %p %p %p\n",ns[0],ns[1],ns[2]);75 printf(" atpointer towards three adjacent triangles\n");76 printf(" at[0] at[1] at[2] = %p %p %p\n",at[0],at[1],at[2]);73 printf(" TriaVertices pointer towards three vertices\n"); 74 printf(" TriaVertices[0] TriaVertices[1] TriaVertices[2] = %p %p %p\n",TriaVertices[0],TriaVertices[1],TriaVertices[2]); 75 printf(" TriaAdjTriangles pointer towards three adjacent triangles\n"); 76 printf(" TriaAdjTriangles[0] TriaAdjTriangles[1] TriaAdjTriangles[2] = %p %p %p\n",TriaAdjTriangles[0],TriaAdjTriangles[1],TriaAdjTriangles[2]); 77 77 printf(" det (integer triangle determinant) = %i\n",det); 78 78 if (link){ … … 85 85 printf("\nThree vertices:\n"); 86 86 for(i=0;i<3;i++){ 87 if ( ns[i]){88 ns[i]->Echo();87 if (TriaVertices[i]){ 88 TriaVertices[i]->Echo(); 89 89 } 90 90 else{ … … 96 96 } 97 97 /*}}}*/ 98 /*FUNCTION Triangle::Optim{{{1*/ 99 Int4 Triangle::Optim(Int2 i,int koption) { 100 // turne around in positif sens 101 Int4 NbSwap =0; 102 Triangle *t = this; 103 int k=0,j =OppositeEdge[i]; 104 int jp = PreviousEdge[j]; 105 // initialise tp, jp the previous triangle & edge 106 Triangle *tp= TriaAdjTriangles[jp]; 107 jp = TriaAdjSharedEdge[jp]&3; 108 do { 109 while (t->swap(j,koption)) 110 { 111 NbSwap++; 112 k++; 113 if (k>=20000){ 114 throw ErrorException(__FUNCT__,exprintf("k>=20000")); 115 } 116 t= tp->TriaAdjTriangles[jp]; // set unchange t qnd j for previous triangles 117 j= NextEdge[tp->TriaAdjSharedEdge[jp]&3]; 118 } 119 // end on this Triangle 120 tp = t; 121 jp = NextEdge[j]; 122 123 t= tp->TriaAdjTriangles[jp]; // set unchange t qnd j for previous triangles 124 j= NextEdge[tp->TriaAdjSharedEdge[jp]&3]; 125 126 } while( t != this); 127 return NbSwap; 128 } 129 /*}}}1*/ 130 /*FUNCTION Triangle::swap{{{1*/ 131 int Triangle::swap(Int2 a,int koption){ 132 if(a/4 !=0) return 0;// arete lock or MarkUnSwap 133 134 register Triangle *t1=this,*t2=TriaAdjTriangles[a];// les 2 triangles adjacent 135 register Int1 a1=a,a2=TriaAdjSharedEdge[a];// les 2 numero de l arete dans les 2 triangles 136 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap 137 138 register Vertex *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]]; 139 register Vertex *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]]; 140 register Vertex *s1=t1->TriaVertices[OppositeVertex[a1]]; 141 register Vertex *s2=t2->TriaVertices[OppositeVertex[a2]]; 142 143 Icoor2 det1=t1->det , det2=t2->det ; 144 Icoor2 detT = det1+det2; 145 Icoor2 detA = Abs(det1) + Abs(det2); 146 Icoor2 detMin = Min(det1,det2); 147 148 int OnSwap = 0; 149 // si 2 triangle infini (bord) => detT = -2; 150 if (sa == 0) {// les deux triangles sont frontieres 151 det2=bamg::det(s2->i,sb->i,s1->i); 152 OnSwap = det2 >0;} 153 else if (sb == 0) { // les deux triangles sont frontieres 154 det1=bamg::det(s1->i,sa->i,s2->i); 155 OnSwap = det1 >0;} 156 else if(( s1 != 0) && (s2 != 0) ) { 157 det1 = bamg::det(s1->i,sa->i,s2->i); 158 det2 = detT - det1; 159 OnSwap = (Abs(det1) + Abs(det2)) < detA; 160 161 Icoor2 detMinNew=Min(det1,det2); 162 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test 163 if (! OnSwap &&(detMinNew>0)) { 164 OnSwap = detMin ==0; 165 if (! OnSwap) { 166 int kopt = koption; 167 while (1) 168 if(kopt) { 169 // critere de Delaunay pure isotrope 170 register Icoor2 xb1 = sb->i.x - s1->i.x, 171 x21 = s2->i.x - s1->i.x, 172 yb1 = sb->i.y - s1->i.y, 173 y21 = s2->i.y - s1->i.y, 174 xba = sb->i.x - sa->i.x, 175 x2a = s2->i.x - sa->i.x, 176 yba = sb->i.y - sa->i.y, 177 y2a = s2->i.y - sa->i.y; 178 register double 179 cosb12 = double(xb1*x21 + yb1*y21), 180 cosba2 = double(xba*x2a + yba*y2a) , 181 sinb12 = double(det2), 182 sinba2 = double(t2->det); 183 184 185 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) 186 OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); 187 break; 188 } 189 else 190 { 191 // critere de Delaunay anisotrope 192 Real8 som; 193 I2 AB=(I2) *sb - (I2) *sa; 194 I2 MAB2=((I2) *sb + (I2) *sa); 195 R2 MAB(MAB2.x*0.5,MAB2.y*0.5); 196 I2 A1=(I2) *s1 - (I2) *sa; 197 I2 D = (I2) * s1 - (I2) * sb ; 198 R2 S2(s2->i.x,s2->i.y); 199 R2 S1(s1->i.x,s1->i.y); 200 { 201 Metric M=s1->m; 202 R2 ABo = M.Orthogonal(AB); 203 R2 A1o = M.Orthogonal(A1); 204 // (A+B)+ x ABo = (S1+B)/2+ y A1 205 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 206 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); 207 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 208 if (Abs(d) > dd*1.e-3) { 209 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 210 som = M(C - S2)/M(C - S1); 211 } else 212 {kopt=1;continue;} 213 214 } 215 { 216 Metric M=s2->m; 217 R2 ABo = M.Orthogonal(AB); 218 R2 A1o = M.Orthogonal(A1); 219 // (A+B)+ x ABo = (S1+B)/2+ y A1 220 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 221 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); 222 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 223 if(Abs(d) > dd*1.e-3) { 224 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 225 som += M(C - S2)/M(C - S1); 226 } else 227 {kopt=1;continue;} 228 } 229 OnSwap = som < 2; 230 break; 231 } 232 233 } // OnSwap 234 } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) 235 } 236 if( OnSwap ) 237 bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); 238 else { 239 t1->SetMarkUnSwap(a1); 240 } 241 return OnSwap; 242 } 243 /*}}}1*/ 98 244 99 245 } -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2969 r2970 3554 3554 printf(" NbSwap of insertion: %i\n",NbSwap); 3555 3555 printf(" NbSwap/nbv: %i\n",NbSwap/nbv); 3556 printf(" NbUnSwap: %i\n",NbUnSwap); 3557 printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); 3558 } 3559 NbUnSwap = 0; 3556 } 3560 3557 3561 3558 #ifdef NBLOOPOPTIM … … 3572 3569 printf(" Optim Loop: %i\n",Nbloop); 3573 3570 printf(" NbSwap/nbv: %i\n",NbSwap/nbv); 3574 printf(" NbUnSwap: %i\n",NbUnSwap); 3575 printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); 3576 } 3577 NbUnSwap = 0; 3571 } 3578 3572 if(!NbSwap) break; 3579 3573 } … … 3963 3957 } 3964 3958 /*}}}1*/ 3965 /*FUNCTION Triangles::Optim{{{1*/3966 Int4 Triangle::Optim(Int2 i,int koption) {3967 // turne around in positif sens3968 Int4 NbSwap =0;3969 Triangle *t = this;3970 int k=0,j =OppositeEdge[i];3971 int jp = PreviousEdge[j];3972 // initialise tp, jp the previous triangle & edge3973 Triangle *tp= at[jp];3974 jp = aa[jp]&3;3975 do {3976 while (t->swap(j,koption))3977 {3978 NbSwap++;3979 k++;3980 if (k>=20000){3981 throw ErrorException(__FUNCT__,exprintf("k>=20000"));3982 }3983 t= tp->at[jp]; // set unchange t qnd j for previous triangles3984 j= NextEdge[tp->aa[jp]&3];3985 }3986 // end on this Triangle3987 tp = t;3988 jp = NextEdge[j];3989 3990 t= tp->at[jp]; // set unchange t qnd j for previous triangles3991 j= NextEdge[tp->aa[jp]&3];3992 3993 } while( t != this);3994 return NbSwap;3995 }3996 /*}}}1*/3997 3959 /*FUNCTION Triangles::PreInit{{{1*/ 3998 3960 void Triangles::PreInit(Int4 inbvx) { … … 5181 5143 } 5182 5144 /*}}}1*/ 5183 /*FUNCTION Triangles::swap{{{1*/5184 int Triangle::swap(Int2 a,int koption){5185 if(a/4 !=0) return 0;// arete lock or MarkUnSwap5186 5187 register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent5188 register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles5189 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap5190 5191 register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]];5192 register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]];5193 register Vertex *s1=t1->ns[OppositeVertex[a1]];5194 register Vertex *s2=t2->ns[OppositeVertex[a2]];5195 5196 Icoor2 det1=t1->det , det2=t2->det ;5197 Icoor2 detT = det1+det2;5198 Icoor2 detA = Abs(det1) + Abs(det2);5199 Icoor2 detMin = Min(det1,det2);5200 5201 int OnSwap = 0;5202 // si 2 triangle infini (bord) => detT = -2;5203 if (sa == 0) {// les deux triangles sont frontieres5204 det2=bamg::det(s2->i,sb->i,s1->i);5205 OnSwap = det2 >0;}5206 else if (sb == 0) { // les deux triangles sont frontieres5207 det1=bamg::det(s1->i,sa->i,s2->i);5208 OnSwap = det1 >0;}5209 else if(( s1 != 0) && (s2 != 0) ) {5210 det1 = bamg::det(s1->i,sa->i,s2->i);5211 det2 = detT - det1;5212 OnSwap = (Abs(det1) + Abs(det2)) < detA;5213 5214 Icoor2 detMinNew=Min(det1,det2);5215 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test5216 if (! OnSwap &&(detMinNew>0)) {5217 OnSwap = detMin ==0;5218 if (! OnSwap) {5219 int kopt = koption;5220 while (1)5221 if(kopt) {5222 // critere de Delaunay pure isotrope5223 register Icoor2 xb1 = sb->i.x - s1->i.x,5224 x21 = s2->i.x - s1->i.x,5225 yb1 = sb->i.y - s1->i.y,5226 y21 = s2->i.y - s1->i.y,5227 xba = sb->i.x - sa->i.x,5228 x2a = s2->i.x - sa->i.x,5229 yba = sb->i.y - sa->i.y,5230 y2a = s2->i.y - sa->i.y;5231 register double5232 cosb12 = double(xb1*x21 + yb1*y21),5233 cosba2 = double(xba*x2a + yba*y2a) ,5234 sinb12 = double(det2),5235 sinba2 = double(t2->det);5236 5237 5238 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2)5239 OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12);5240 break;5241 }5242 else5243 {5244 // critere de Delaunay anisotrope5245 Real8 som;5246 I2 AB=(I2) *sb - (I2) *sa;5247 I2 MAB2=((I2) *sb + (I2) *sa);5248 R2 MAB(MAB2.x*0.5,MAB2.y*0.5);5249 I2 A1=(I2) *s1 - (I2) *sa;5250 I2 D = (I2) * s1 - (I2) * sb ;5251 R2 S2(s2->i.x,s2->i.y);5252 R2 S1(s1->i.x,s1->i.y);5253 {5254 Metric M=s1->m;5255 R2 ABo = M.Orthogonal(AB);5256 R2 A1o = M.Orthogonal(A1);5257 // (A+B)+ x ABo = (S1+B)/2+ y A15258 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/25259 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);5260 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/25261 if (Abs(d) > dd*1.e-3) {5262 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));5263 som = M(C - S2)/M(C - S1);5264 } else5265 {kopt=1;continue;}5266 5267 }5268 {5269 Metric M=s2->m;5270 R2 ABo = M.Orthogonal(AB);5271 R2 A1o = M.Orthogonal(A1);5272 // (A+B)+ x ABo = (S1+B)/2+ y A15273 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/25274 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);5275 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/25276 if(Abs(d) > dd*1.e-3) {5277 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));5278 som += M(C - S2)/M(C - S1);5279 } else5280 {kopt=1;continue;}5281 }5282 OnSwap = som < 2;5283 break;5284 }5285 5286 } // OnSwap5287 } // (! OnSwap &&(det1 > 0) && (det2 > 0) )5288 }5289 if( OnSwap )5290 bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2);5291 else {5292 NbUnSwap ++;5293 t1->SetMarkUnSwap(a1);5294 }5295 return OnSwap;5296 }5297 /*}}}1*/5298 5145 /*FUNCTION Triangles::UnCrack{{{1*/ 5299 5146 int Triangles::UnCrack() {
Note:
See TracChangeset
for help on using the changeset viewer.
