source: issm/trunk-jpl/src/c/objects/Bamg/Triangle.cpp@ 12493

Last change on this file since 12493 was 12493, checked in by Mathieu Morlighem, 13 years ago

replaced all _error_ to _error2_
File size: 14.5 KB
RevLine 
[3913]1#include <cstdio>
2#include <cstring>
3#include <cmath>
4#include <ctime>
5
6#include "../objects.h"
7
8namespace bamg {
9
10 /*Constructors/Destructors*/
[12365]11 /*FUNCTION Triangle(){{{*/
[9371]12 Triangle::Triangle(void){
13
14 }
15 /*}}}*/
[12365]16 /*FUNCTION Triangle(Mesh *Th,long i,long j,long k) {{{*/
[5095]17 Triangle::Triangle(Mesh *Th,long i,long j,long k) {
[5120]18 BamgVertex *v=Th->vertices;
[3913]19 long nbv = Th->nbv;
20 if (i<0 || j<0 || k<0){
[12493]21 _error2_("i<0 || j<0 || k<0");
[3913]22 }
23 if (i>=nbv || j>=nbv || k>=nbv){
[12493]24 _error2_("i>=nbv || j>=nbv || k>=nbv");
[3913]25 }
[5148]26 vertices[0]=v+i;
27 vertices[1]=v+j;
28 vertices[2]=v+k;
29 adj[0]=adj[1]=adj[2]=0;
[5401]30 AdjEdgeIndex[0]=AdjEdgeIndex[1]=AdjEdgeIndex[2]=0;
[3913]31 det=0;
32 }
33 /*}}}*/
[12365]34 /*FUNCTION Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2) {{{*/
[5120]35 Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){
[5148]36 vertices[0]=v0;
37 vertices[1]=v1;
38 vertices[2]=v2;
39 adj[0]=adj[1]=adj[2]=0;
[5401]40 AdjEdgeIndex[0]=AdjEdgeIndex[1]=AdjEdgeIndex[2]=0;
[3913]41 if (v0) det=0;
42 else {
43 det=-1;
44 link=NULL;};
45 }
46 /*}}}*/
47
48 /*Methods*/
[12365]49 /*FUNCTION Triangle::Adj{{{*/
[5151]50 AdjacentTriangle Triangle::Adj(int i) const {
[5401]51 return AdjacentTriangle(adj[i],AdjEdgeIndex[i]&3);
[5151]52 };/*}}}*/
[12365]53 /*FUNCTION Triangle::Anisotropy{{{*/
[11584]54 double Triangle::Anisotropy() const{
55
56 double lmin,lmax;
57
58 /*Get three vertices A,B and C*/
59 R2 A=*this->vertices[0];
60 R2 B=*this->vertices[1];
61 R2 C=*this->vertices[2];
62
63 /*Compute edges*/
64 R2 e1=B-A;
65 R2 e2=C-A;
66 R2 e3=B-C;
67
68 /*Compute edge length*/
69 double l1=Norme2(e1);
70 double l2=Norme2(e2);
71 double l3=Norme2(e3);
72
73 lmin=l1;
74 lmin=min(lmin,l2);
75 lmin=min(lmin,l3);
76 lmax=l1;
77 lmax=max(lmax,l2);
78 lmax=max(lmax,l3);
79
80 return lmax/lmin;
81 };/*}}}*/
[12365]82 /*FUNCTION Triangle::Length{{{*/
[11584]83 double Triangle::Length() const{
84
85 double l;
86
87 /*Get three vertices A,B and C*/
88 R2 A=*this->vertices[0];
89 R2 B=*this->vertices[1];
90 R2 C=*this->vertices[2];
91
92 /*Compute edges*/
93 R2 e1=B-A;
94 R2 e2=C-A;
95 R2 e3=B-C;
96
97 /*Compute edge length*/
98 l=Norme2(e1);
99 l=max(l,Norme2(e2));
100 l=max(l,Norme2(e3));
101
102 return l;
103 };/*}}}*/
[12365]104 /*FUNCTION Triangle::Echo {{{*/
[3913]105 void Triangle::Echo(void){
106
107 int i;
108
109 printf("Triangle:\n");
[5148]110 printf(" vertices pointer towards three vertices\n");
111 printf(" vertices[0] vertices[1] vertices[2] = %p %p %p\n",vertices[0],vertices[1],vertices[2]);
112 printf(" adj pointer towards three adjacent triangles\n");
113 printf(" adj[0] adj[1] adj[2] = %p %p %p\n",adj[0],adj[1],adj[2]);
[3913]114 printf(" det (integer triangle determinant) = %i\n",det);
115 if (link){
116 printf(" link (pointer toward duplicate triangle)= %p\n",link);
117 }
118 else{
119 printf(" color = %i\n",color);
120 }
121
122 printf("\nThree vertices:\n");
123 for(i=0;i<3;i++){
[5148]124 if (vertices[i]){
125 vertices[i]->Echo();
[3913]126 }
127 else{
128 printf(" vertex %i does not exist\n",i+1);
129 }
130 }
131
132 return;
133 }
134 /*}}}*/
[12365]135 /*FUNCTION Triangle::FindBoundaryEdge{{{*/
[5143]136 AdjacentTriangle Triangle::FindBoundaryEdge(int i) const{
[3913]137 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindBoundaryEdge)*/
138
139 /*Intermediary*/
140 Triangle* ttc=NULL;
141 int k,j,jc;
142
143 // call current triangle t
144 Triangle* t = (Triangle*)this;
145
146 //is the current triangle inside or outside?
147 int outside=!link ;
148
149 // EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}};
150 // initialize j as the first vertex of the ith edge
151 j=EdgesVertexTriangle[i][0];
152
153 //Loop over the adjacent triangle of t
154 k=0;
155 do{
156 //keep track of outside
157 int outsidep = outside;
158 //increment k
159 k++;
160 //Get ttc, adjacent triangle of t with respect to vertex j
[5148]161 ttc = t->adj[j];
[3913]162 //is the current triangle inside or outside?
163 outside = !ttc->link;
164 //if both previous triangle are outside, return
[5143]165 if (outside+outsidep == 1) return AdjacentTriangle(t,j);
[3913]166
167 //update t and j
168 t = ttc;
169 //NextEdge[3] = {1,2,0};
[5401]170 jc = NextEdge[t->AdjEdgeIndex[j]&3];
[3913]171 j = NextEdge[jc];
172
173 //check number of iterations
174 if (k>=2000){
[12493]175 _error2_("too many iteration in Triangle::FindBoundaryEdge (k>=2000)");
[3913]176 }
177 } while (this!= t);
178 //not found, return empty triangle
[5143]179 return AdjacentTriangle(NULL,0);
[3913]180 }
[12365]181 /*}}}*/
182 /*FUNCTION Triangle::GetAllflag{{{*/
[5151]183 int Triangle::GetAllflag(int a){
[5401]184 return AdjEdgeIndex[a] & 1020;
[5151]185 }/*}}}*/
[12365]186 /*FUNCTION Triangle::Hidden{{{*/
[5151]187 int Triangle::Hidden(int a)const {
[5401]188 return AdjEdgeIndex[a]&16;
[5151]189 } /*}}}*/
[12365]190 /*FUNCTION Triangle::Locked{{{*/
[5151]191 int Triangle::Locked(int a)const {
[5401]192 return AdjEdgeIndex[a]&4;
[5151]193 } /*}}}*/
[12365]194 /*FUNCTION Triangle::NuEdgeTriangleAdj{{{*/
[5151]195 short Triangle::NuEdgeTriangleAdj(int i) const {
[5187]196 /*Number of the adjacent edge in adj tria (make sure it is between 0 and 2*/
[5401]197 return AdjEdgeIndex[i&3]&3;
[5151]198 }/*}}}*/
[12365]199 /*FUNCTION Triangle::Optim{{{*/
[3913]200 long Triangle::Optim(short i,int koption) {
201 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Optim)*/
202
203 // turn around (positive direction)
204 Triangle *t=this;
205 long NbSwap =0;
206 int k = 0;
207 int j = OppositeEdge[i];
208 int jp= PreviousEdge[j];
209
210 // initialize tp, jp the previous triangle & edge
[5148]211 Triangle *tp=adj[jp];
[5401]212 jp = AdjEdgeIndex[jp]&3;
[3913]213 do {
214 while (t->swap(j,koption)){
[12493]215 if (k>=20000) _error2_("k>=20000");
[3913]216 NbSwap++;
217 k++;
[5148]218 t= tp->adj[jp]; // set unchange t qnd j for previous triangles
[5401]219 j= NextEdge[tp->AdjEdgeIndex[jp]&3];
[3913]220 }
221 // end on this Triangle
222 tp = t;
223 jp = NextEdge[j];
224
[5148]225 t= tp->adj[jp]; // set unchange t qnd j for previous triangles
[5401]226 j= NextEdge[tp->AdjEdgeIndex[jp]&3];
[3913]227
228 } while( t != this);
229 return NbSwap;
230 }
[12365]231 /*}}}*/
232 /*FUNCTION Triangle::Quadrangle {{{*/
[5120]233 Triangle* Triangle::Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const{
[3913]234 // return the other triangle of the quad if a quad or 0 if not a quat
235 Triangle * t =0;
236 if (link) {
237 int a=-1;
[5401]238 if (AdjEdgeIndex[0] & 16 ) a=0;
239 if (AdjEdgeIndex[1] & 16 ) a=1;
240 if (AdjEdgeIndex[2] & 16 ) a=2;
[3913]241 if (a>=0) {
[5148]242 t = adj[a];
[3913]243 // if (t-this<0) return 0;
[5148]244 v2 = vertices[VerticesOfTriangularEdge[a][0]];
245 v0 = vertices[VerticesOfTriangularEdge[a][1]];
246 v1 = vertices[OppositeEdge[a]];
[5401]247 v3 = t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]];
[3913]248 }
249 }
250 return t;
251 }
252 /*}}}*/
[12365]253 /*FUNCTION Triangle::QualityQuad {{{*/
[3913]254 double Triangle::QualityQuad(int a,int option) const{
255 double q;
[5401]256 if (!link || AdjEdgeIndex[a] &4)
[3913]257 q= -1;
258 else {
[5148]259 Triangle * t = adj[a];
[3913]260 if (t-this<0) q= -1;// because we do 2 times
261 else if (!t->link ) q= -1;
[5401]262 else if (AdjEdgeIndex[0] & 16 || AdjEdgeIndex[1] & 16 || AdjEdgeIndex[2] & 16 || t->AdjEdgeIndex[0] & 16 || t->AdjEdgeIndex[1] & 16 || t->AdjEdgeIndex[2] & 16 )
[3913]263 q= -1;
264 else if(option){
[5148]265 const BamgVertex & v2 = *vertices[VerticesOfTriangularEdge[a][0]];
266 const BamgVertex & v0 = *vertices[VerticesOfTriangularEdge[a][1]];
267 const BamgVertex & v1 = *vertices[OppositeEdge[a]];
[5401]268 const BamgVertex & v3 = * t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]];
[3913]269 q = QuadQuality(v0,v1,v2,v3); // do the float part
270 }
271 else q= 1;
272 }
273 return q;
274 }
275 /*}}}*/
[12365]276 /*FUNCTION Triangle::Renumbering(Triangle *tb,Triangle *te, long *renu){{{*/
[5151]277 void Triangle::Renumbering(Triangle *tb,Triangle *te, long *renu){
278
279 if (link >=tb && link <te) link = tb + renu[link -tb];
280 if (adj[0] >=tb && adj[0] <te) adj[0] = tb + renu[adj[0]-tb];
281 if (adj[1] >=tb && adj[1] <te) adj[1] = tb + renu[adj[1]-tb];
282 if (adj[2] >=tb && adj[2] <te) adj[2] = tb + renu[adj[2]-tb];
283 }/*}}}*/
[12365]284 /*FUNCTION Triangle::Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){{{*/
[5151]285 void Triangle::Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
286 if (vertices[0] >=vb && vertices[0] <ve) vertices[0] = vb + renu[vertices[0]-vb];
287 if (vertices[1] >=vb && vertices[1] <ve) vertices[1] = vb + renu[vertices[1]-vb];
288 if (vertices[2] >=vb && vertices[2] <ve) vertices[2] = vb + renu[vertices[2]-vb];
289 }/*}}}*/
[12365]290 /*FUNCTION Triangle::Set {{{*/
[5095]291 void Triangle::Set(const Triangle & rec,const Mesh & Th ,Mesh & ThNew){
[3913]292 *this = rec;
[5149]293 if ( vertices[0] ) vertices[0] = ThNew.vertices + Th.GetId(vertices[0]);
294 if ( vertices[1] ) vertices[1] = ThNew.vertices + Th.GetId(vertices[1]);
295 if ( vertices[2] ) vertices[2] = ThNew.vertices + Th.GetId(vertices[2]);
296 if(adj[0]) adj[0] = ThNew.triangles + Th.GetId(adj[0]);
297 if(adj[1]) adj[1] = ThNew.triangles + Th.GetId(adj[1]);
298 if(adj[2]) adj[2] = ThNew.triangles + Th.GetId(adj[2]);
[3913]299 if (link >= Th.triangles && link < Th.triangles + Th.nbt)
[5149]300 link = ThNew.triangles + Th.GetId(link);
[3913]301 }
302 /*}}}*/
[12365]303 /*FUNCTION Triangle::SetAdjAdj{{{*/
[5151]304 void Triangle::SetAdjAdj(short a){
305 // Copy all the mark
306 a &= 3;
307 register Triangle *tt=adj[a];
[5401]308 AdjEdgeIndex [a] &= 55; // remove MarkUnSwap
309 register short aatt = AdjEdgeIndex[a] & 3;
[5151]310 if(tt){
311 tt->adj[aatt]=this;
[5401]312 tt->AdjEdgeIndex[aatt]=a + (AdjEdgeIndex[a] & 60 ) ;
[5151]313 }
314 }/*}}}*/
[12365]315 /*FUNCTION Triangle::SetAdj2{{{*/
[5151]316 void Triangle::SetAdj2(short a,Triangle *t,short aat){
[5605]317 /*For current triangle:
318 * - a is the index of the edge were the adjency is set (in [0 2])
319 * - t is the adjacent triangle
320 * - aat is the index of the same edge in the adjacent triangle*/
321 adj[a]=t;
322 AdjEdgeIndex[a]=aat;
323 if(t){ //if t!=NULL add adjacent triangle to t (this)
[5151]324 t->adj[aat]=this;
[5401]325 t->AdjEdgeIndex[aat]=a;
[5151]326 }
327 }/*}}}*/
[12365]328 /*FUNCTION Triangle::SetAllFlag{{{*/
[5151]329 void Triangle::SetAllFlag(int a,int f){
[5401]330 AdjEdgeIndex[a] = (AdjEdgeIndex[a] &3) + (1020 & f);
[5151]331 }/*}}}*/
[12365]332 /*FUNCTION Triangle::SetDet{{{*/
[5151]333 void Triangle::SetDet() {
334 if(vertices[0] && vertices[1] && vertices[2]) det = bamg::det(*vertices[0],*vertices[1],*vertices[2]);
335 else det = -1;
336 }/*}}}*/
[12365]337 /*FUNCTION Triangle::SetHidden{{{*/
[5151]338 void Triangle::SetHidden(int a){
339 //Get Adjacent Triangle number a
340 register Triangle* t = adj[a];
341 //if it exist
342 //C|=D -> C=(C|D) bitwise inclusive OR
[5401]343 if(t) t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=16;
344 AdjEdgeIndex[a] |= 16;
[5151]345 }/*}}}*/
[12365]346 /*FUNCTION Triangle::SetLocked{{{*/
[5151]347 void Triangle::SetLocked(int a){
348 //mark the edge as on Boundary
349 register Triangle * t = adj[a];
[5401]350 t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=4;
351 AdjEdgeIndex[a] |= 4;
[5151]352 }/*}}}*/
[12365]353 /*FUNCTION Triangle::SetMarkUnSwap{{{*/
[5151]354 void Triangle::SetMarkUnSwap(int a){
355 register Triangle * t = adj[a];
[5401]356 t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=8;
357 AdjEdgeIndex[a] |=8 ;
[5151]358 }/*}}}*/
[12365]359 /*FUNCTION Triangle::SetSingleVertexToTriangleConnectivity{{{*/
[5151]360 void Triangle::SetSingleVertexToTriangleConnectivity() {
[5460]361 if (vertices[0]) (vertices[0]->t=this,vertices[0]->IndexInTriangle=0);
362 if (vertices[1]) (vertices[1]->t=this,vertices[1]->IndexInTriangle=1);
363 if (vertices[2]) (vertices[2]->t=this,vertices[2]->IndexInTriangle=2);
[5151]364 }/*}}}*/
[12365]365 /*FUNCTION Triangle::SetUnMarkUnSwap{{{*/
[5151]366 void Triangle::SetUnMarkUnSwap(int a){
367 register Triangle * t = adj[a];
[5401]368 t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] &=55; // 23 + 32
369 AdjEdgeIndex[a] &=55 ;
[5151]370 }/*}}}*/
[12365]371 /*FUNCTION Triangle::swap{{{*/
[3913]372 int Triangle::swap(short a,int koption){
373 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/
374
375 if(a/4 !=0) return 0;// arete lock or MarkUnSwap
376
[5148]377 register Triangle *t1=this,*t2=adj[a];// les 2 triangles adjacent
[5401]378 register short a1=a,a2=AdjEdgeIndex[a];// les 2 numero de l arete dans les 2 triangles
[3913]379 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
380
[5148]381 register BamgVertex *sa=t1->vertices[VerticesOfTriangularEdge[a1][0]];
382 register BamgVertex *sb=t1->vertices[VerticesOfTriangularEdge[a1][1]];
383 register BamgVertex *s1=t1->vertices[OppositeVertex[a1]];
384 register BamgVertex *s2=t2->vertices[OppositeVertex[a2]];
[3913]385
386 Icoor2 det1=t1->det , det2=t2->det ;
387 Icoor2 detT = det1+det2;
388 Icoor2 detA = Abs(det1) + Abs(det2);
389 Icoor2 detMin = Min(det1,det2);
390
391 int OnSwap = 0;
392 // si 2 triangle infini (bord) => detT = -2;
393 if (sa == 0) {// les deux triangles sont frontieres
394 det2=bamg::det(s2->i,sb->i,s1->i);
395 OnSwap = det2 >0;}
396 else if (sb == 0) { // les deux triangles sont frontieres
397 det1=bamg::det(s1->i,sa->i,s2->i);
398 OnSwap = det1 >0;}
399 else if(( s1 != 0) && (s2 != 0) ) {
400 det1 = bamg::det(s1->i,sa->i,s2->i);
401 det2 = detT - det1;
402 OnSwap = (Abs(det1) + Abs(det2)) < detA;
403
404 Icoor2 detMinNew=Min(det1,det2);
405 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test
406 if (! OnSwap &&(detMinNew>0)) {
407 OnSwap = detMin ==0;
408 if (! OnSwap) {
409 int kopt = koption;
410 while (1)
411 if(kopt) {
412 // critere de Delaunay pure isotrope
413 register Icoor2 xb1 = sb->i.x - s1->i.x,
414 x21 = s2->i.x - s1->i.x,
415 yb1 = sb->i.y - s1->i.y,
416 y21 = s2->i.y - s1->i.y,
417 xba = sb->i.x - sa->i.x,
418 x2a = s2->i.x - sa->i.x,
419 yba = sb->i.y - sa->i.y,
420 y2a = s2->i.y - sa->i.y;
421 register double
422 cosb12 = double(xb1*x21 + yb1*y21),
423 cosba2 = double(xba*x2a + yba*y2a) ,
424 sinb12 = double(det2),
425 sinba2 = double(t2->det);
426
427
428 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2)
429 OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12);
430 break;
431 }
432 else {
433 // critere de Delaunay anisotrope
434 double som;
435 I2 AB=(I2) *sb - (I2) *sa;
436 I2 MAB2=((I2) *sb + (I2) *sa);
437 R2 MAB(MAB2.x*0.5,MAB2.y*0.5);
438 I2 A1=(I2) *s1 - (I2) *sa;
439 I2 D = (I2) * s1 - (I2) * sb ;
440 R2 S2(s2->i.x,s2->i.y);
441 R2 S1(s1->i.x,s1->i.y);
442 {
443 Metric M=s1->m;
444 R2 ABo = M.Orthogonal(AB);
445 R2 A1o = M.Orthogonal(A1);
446 // (A+B)+ x ABo = (S1+B)/2+ y A1
447 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
448 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
449 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
450 if (Abs(d) > dd*1.e-3) {
451 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
452 som = M(C - S2)/M(C - S1);
453 } else
454 {kopt=1;continue;}
455
456 }
457 {
458 Metric M=s2->m;
459 R2 ABo = M.Orthogonal(AB);
460 R2 A1o = M.Orthogonal(A1);
461 // (A+B)+ x ABo = (S1+B)/2+ y A1
462 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
463 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
464 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
465 if(Abs(d) > dd*1.e-3) {
466 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
467 som += M(C - S2)/M(C - S1);
468 } else
469 {kopt=1;continue;}
470 }
471 OnSwap = som < 2;
472 break;
473 }
474
475 } // OnSwap
476 } // (! OnSwap &&(det1 > 0) && (det2 > 0) )
477 }
478 if( OnSwap )
479 bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2);
480 else {
481 t1->SetMarkUnSwap(a1);
482 }
483 return OnSwap;
484 }
[12365]485 /*}}}*/
486 /*FUNCTION Triangle::TriangleAdj{{{*/
[5151]487 Triangle* Triangle::TriangleAdj(int i) const {
488 return adj[i&3];
489 }/*}}}*/
[3913]490
491}
Note: See TracBrowser for help on using the repository browser.