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

Last change on this file since 10937 was 9371, checked in by Mathieu Morlighem, 14 years ago

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