source: issm/trunk-jpl/src/c/classes/bamg/Triangle.cpp@ 15012

Last change on this file since 15012 was 15012, checked in by Mathieu Morlighem, 12 years ago

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