1 // -*- Mode : c++ -*-
2 //
3 // SUMMARY :
4 // USAGE :
5 // ORG :
6 // AUTHOR : Frederic Hecht
7 // E-MAIL : hecht@ann.jussieu.fr
8 //
9
10 /*
11
12 This file is part of Freefem++
13
14 Freefem++ is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 2.1 of the License, or
17 (at your option) any later version.
18
19 Freefem++ is distributed in the hope that it will be useful,
20 but WITHOUT ANY WARRANTY; without even the implied warranty of
21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 GNU Lesser General Public License for more details.
23
24 You should have received a copy of the GNU Lesser General Public License
25 along with Freefem++; if not, write to the Free Software
26 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 */
28 #ifdef __MWERKS__
29 #ifdef __INTEL__
30 //#pragma global_optimizer off
31 //#pragma inline_depth(0)
32 //#pragma optimization_level 2
33 #endif
34 //#pragma inline_depth 0
35 #endif
36 extern bool withrgraphique;
37 #include <stdio.h>
38 #include <string.h>
39 #include <math.h>
40 #include <time.h>
41 #include <iostream>
42 using namespace std;
43
44 #include "Mesh2.h"
45 #include "QuadTree.h"
46 #include "SetOfE4.h"
47
48 namespace bamg {
49
50
51 #ifdef DEBUG1
52 extern int SHOW ; // for debugging
53 int SHOW = 0; // for debugging
54
55 #endif
56
57 int Triangles::counter = 0;
58
59 Triangles * CurrentTh =0;
60
61 int hinterpole=1;
62
63
64 long NbUnSwap =0;
65 int ForDebugging = 0;
66 const Direction NoDirOfSearch = Direction();
67 #ifndef NDEBUG
MyAssert(int i,char * ex,char * file,long line)68 inline void MyAssert(int i,char*ex,char * file,long line)
69 {
70 if( i) {
71 cerr << "Error Assert:" << ex << " in " << file << " line: " << line << endl;
72 #ifdef NOTFREEFEM
73 exit(1);
74 #else
75 throw(ErrorExec("exit",1000));
76 #endif
77 }
78 }
79 #endif
80
AGoodNumberPrimeWith(Int4 n)81 Int4 AGoodNumberPrimeWith(Int4 n)
82 {
83 const Int4 BigPrimeNumber[] ={ 567890359L,
84 567890431L, 567890437L, 567890461L, 567890471L,
85 567890483L, 567890489L, 567890497L, 567890507L,
86 567890591L, 567890599L, 567890621L, 567890629L , 0};
87
88 Int4 o = 0;
89 Int4 pi = BigPrimeNumber[1];
90 for (int i=0; BigPrimeNumber[i]; i++) {
91 Int4 r = BigPrimeNumber[i] % n;
92 Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r)));
93 if ( o < oo)
94 o=oo,pi=BigPrimeNumber[i];}
95 // cout << " AGoodNumberPrimeWith " << n << " " <<pi << " "<< o << endl;
96 return pi;
97 }
98
99 class Triangles;
MeshError(int Err,Triangles * Th)100 void MeshError(int Err,Triangles *Th){
101 cerr << " Fatal error in the meshgenerator " << Err << endl ;
102 #ifdef NOTFREEFEM
103 exit(1);
104 #else
105 throw(ErrorMesh("Bamg",Err,Th));
106 #endif
107 }
108
operator <<(ostream & f,const Triangle & ta)109 ostream& operator <<(ostream& f, const Triangle & ta)
110 {
111 if(CurrentTh)
112 f << "[" << CurrentTh->Number(ta) << "::"
113 << CurrentTh->Number(ta.ns[0]) << ","
114 << CurrentTh->Number(ta.ns[1]) << ","
115 << CurrentTh->Number(ta.ns[2]) << ","
116 << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} "
117 << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} "
118 << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} "
119 << "]" ;
120 else
121 f << "["
122 << ta.ns[0] << ","
123 << ta.ns[1] << ","
124 << ta.ns[2] << ","
125 << "{" << ta.at[0] << " " << ta.aa[0] << "} "
126 << "{" << ta.at[1] << " " << ta.aa[1] << "} "
127 << "{" << ta.at[2] << " " << ta.aa[2] << "} "
128 << "]" ;
129 return f;}
130
swap(Triangle * t1,Int1 a1,Triangle * t2,Int1 a2,Vertex * s1,Vertex * s2,Icoor2 det1,Icoor2 det2)131 void swap(Triangle *t1,Int1 a1,
132 Triangle *t2,Int1 a2,
133 Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2)
134 { // swap
135 // --------------------------------------------------------------
136 // Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles
137 //
138 // sb sb
139 // / | \ / \ !
140 // as1/ | \ /a2 \ !
141 // / | \ / t2 \ !
142 // s1 /t1 | t2 \s2 --> s1 /___as2___\s2 !
143 // \ a1|a2 / \ as1 /
144 // \ | / \ t1 /
145 // \ | / as2 \ a1/
146 // \ | / \ /
147 // sa sa
148 // -------------------------------------------------------------
149 int as1 = NextEdge[a1];
150 int as2 = NextEdge[a2];
151 int ap1 = PreviousEdge[a1];
152 int ap2 = PreviousEdge[a2];
153 #ifdef DRAWING1
154 couleur(0);
155 t1->Draw();
156 t2->Draw();
157 #endif
158 #ifdef DEBUG1
159 t1->check();
160 t2->check();
161 #endif
162 (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb
163 (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa
164 // mise a jour des 2 adjacences externes
165 TriangleAdjacent taas1 = t1->Adj(as1),
166 taas2 = t2->Adj(as2),
167 tas1(t1,as1), tas2(t2,as2),
168 ta1(t1,a1),ta2(t2,a2);
169 #ifdef DEBUG
170 assert( ! ta1.Locked());
171 assert( ! ta2.Locked());
172 #endif
173 // externe haut gauche
174 taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap());
175 // externe bas droite
176 taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap());
177 // remove the Mark UnMarkSwap
178 t1->SetUnMarkUnSwap(ap1);
179 t2->SetUnMarkUnSwap(ap2);
180 // interne
181 tas1.SetAdj2(tas2);
182
183 t1->det = det1;
184 t2->det = det2;
185
186 t1->SetTriangleContainingTheVertex();
187 t2->SetTriangleContainingTheVertex();
188 #ifdef DEBUG1
189 t1->check();
190 t2->check();
191 #endif
192 #ifdef DRAWING1
193 couleur(1);
194 t1->Draw();
195 t2->Draw();
196 #endif
197 #ifdef DRAWING1
198 if( CurrentTh)
199 CurrentTh->inquire();
200 #endif
201
202 } // end swap
203
204
205
206
207
FindTriangle(Triangles & Th,Real8 x,Real8 y,double * a,int & inside)208 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside)
209 {
210 CurrentTh=&Th;
211 assert(&Th);
212 I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y)));
213 Icoor2 dete[3];
214 Triangle & tb = *Th.FindTriangleContening(I,dete);
215
216 if (tb.link)
217 { // internal point in a true triangles
218 a[0]= (Real8) dete[0]/ tb.det;
219 a[1]= (Real8) dete[1] / tb.det;
220 a[2] = (Real8) dete[2] / tb.det;
221 inside = 1;
222 return Th.Number(tb);
223 }
224 else
225 {
226 inside = 0;
227 double aa,bb;
228 TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb);
229 int k = ta;
230 Triangle * tc = ta;
231 if (!tc->link)
232 { ta = ta.Adj();
233 tc=ta;
234 k = ta;
235 Exchange(aa,bb);
236 assert(tc->link);
237 }
238 a[VerticesOfTriangularEdge[k][0]] = aa;
239 a[VerticesOfTriangularEdge[k][1]] = bb;
240 a[OppositeVertex[k]] = 1- aa -bb;
241 return Th.Number(tc);
242 }
243 }
244
245
CloseBoundaryEdge(I2 A,Triangle * t,double & a,double & b)246 TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) {
247 //
248 // cout << " - ";
249 int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
250 int dir=0;
251 assert(k>=0);
252 int kkk=0;
253 Icoor2 IJ_IA,IJ_AJ;
254 TriangleAdjacent edge(t,OppositeEdge[k]);
255 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge))))
256 {
257
258 assert(kkk++<1000);
259 Vertex &vI = *edge.EdgeVertex(0);
260 Vertex &vJ = *edge.EdgeVertex(1);
261 I2 I=vI, J=vJ, IJ= J-I;
262 IJ_IA = (IJ ,(A-I));
263 // cout << A << vI.i << vJ.i << edge << " " << IJ_IA << " dir " << dir <<endl;
264 if (IJ_IA<0) {
265 if (dir>0) {a=1;b=0;return edge;}// change of signe => I
266 else {dir=-1;
267 continue;}};// go in direction i
268 IJ_AJ = (IJ ,(J-A));
269 if (IJ_AJ<0) {
270 if(dir<0) {a=0;b=1;return edge;}
271 else {dir = 1;
272 continue;}}// go in direction j
273 double IJ2 = IJ_IA + IJ_AJ;
274 assert(IJ2);
275 a= IJ_AJ/IJ2;
276 b= IJ_IA/IJ2;
277 // cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl;
278 return edge;
279 }
280 }
281
FindBoundaryEdge(int i) const282 TriangleAdjacent Triangle::FindBoundaryEdge(int i) const
283 {
284 // turn around the vertex ns[i] also call s
285 #ifdef DEBUG
286 register Vertex * s = ns[i];
287 #endif
288 Triangle *t = (Triangle *) this , *ttc;
289 int k=0,j = EdgesVertexTriangle[i][0],jc;
290 int exterieur = !link ;
291
292 do
293 {
294 int exterieurp = exterieur;
295 k++;
296 #ifdef DEBUG
297 assert( s == & (*t)[VerticesOfTriangularEdge[j][1]] );
298 #endif
299 ttc = t->at[j];
300 exterieur = !ttc->link;
301 if (exterieur+exterieurp == 1)
302 return TriangleAdjacent(t,j);
303 jc = NextEdge[t->aa[j]&3];
304 t = ttc;
305 j = NextEdge[jc];
306 assert(k<2000);
307 } while ( (this!= t));
308 return TriangleAdjacent(0,0);
309
310 }
311
312
CloseBoundaryEdgeV2(I2 C,Triangle * t,double & a,double & b)313 TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
314 {
315 // walk around the vertex
316 // version 2 for remove the probleme if we fill the hole
317 //int bug=1;
318 // Triangle *torigine = t;
319 // restart:
320 // int dir=0;
321 assert(t->link == 0);
322 // to have a starting edges
323 // try the 3 edge bourna-- in case of internal hole
324 // and choice the best
325 //
326 //
327 // the probleme is in case of the fine and long internal hole
328 // for exemple neart the training edge of a wing
329 //
330 Vertex * s=0,*s1=0, *s0=0;
331 Icoor2 imax = MaxICoor22;
332 Icoor2 l0 = imax,l1 = imax;
333 double dd2 = imax;// infinity
334 TriangleAdjacent er;
335 int cas=-2;
336 for (int j=0;j<3;j++)
337 {
338 TriangleAdjacent ta=t->FindBoundaryEdge(j);
339 if (! (Triangle *) ta) continue;
340 s0 = ta.EdgeVertex(0);
341 s1 = ta.EdgeVertex(1);
342 I2 A = * s0;
343 I2 B = *ta.EdgeVertex(1);
344 I2 AB = B-A,AC=C-A,BC=B-C;
345 Icoor2 ACAC = (AC,AC), BCBC = (BC,BC);
346 Icoor2 AB2 = Norme2_2(AB); // ||AB||^2
347 Icoor2 ABAC = (AB,AC); // AB.AC|
348
349 double d2;
350 if ( ABAC < 0 ) // DIST A
351 {
352 if ( (d2=(double) ACAC) < dd2)
353 {
354 // cout << " A " << d2 << " " << dd2;
355 er = ta;
356 l0 = ACAC;
357 l1 = BCBC;
358 cas = 0;
359 s = s0;
360 }
361 }
362 else if (ABAC > AB2) // DIST B
363 {
364 if ( (d2=(double) BCBC) < dd2)
365 {
366 // cout << " B " << d2 << " " << dd2;
367 dd2 = d2;
368 er = Adj(ta); // other direction
369 l0 = BCBC;
370 l1 = ACAC;
371 cas = 1;
372 s = s1;
373 }
374 }
375 else // DIST AB
376 {
377
378 double det_2 = (double) Det(AB,AC);
379 det_2 *= det_2; // square of area*2 of triangle ABC
380 d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC
381 // cout << " AB " << d2 << " " << dd2
382 // << " " << CurrentTh->Number(ta.EdgeVertex(0))
383 // << " " << CurrentTh->Number(ta.EdgeVertex(1)) << " " ;
384
385 if (d2 < dd2)
386 {
387 dd2 = d2;
388 er = ta;
389 l0 = (AC,AC);
390 l1 = (BC,BC);
391 s = 0;
392 cas = -1;
393 // cout << " ABAC " << ABAC << " ABAC " << ABAC
394 // << " AB2 " << AB2 << endl;
395 b = ((double) ABAC/(double) AB2);
396 a = 1 - b;
397 }
398 }
399 }
400 assert(cas !=-2);
401 // l1 = ||C s1|| , l0 = ||C s0||
402 // where s0,s1 are the vertex of the edge er
403
404 if ( s)
405 {
406 t=er;
407 TriangleAdjacent edge(er);
408
409 int kkk=0;
410 int linkp = t->link == 0;
411
412 Triangle * tt=t=edge=Adj(Previous(edge));
413 // cout << CurrentTh->Number(t) << " " << linkp << endl;
414 do { // loop around vertex s
415
416 assert(edge.EdgeVertex(0)==s && kkk++<10000);
417
418 int link = tt->link == 0;
419 // cout << CurrentTh->Number(tt) << " " << link << " " << CurrentTh->Number(s)
420 // << " " << CurrentTh->Number(er.EdgeVertex(0))
421 // << " " << CurrentTh->Number(er.EdgeVertex(1))
422 // << " " << CurrentTh->Number(edge.EdgeVertex(0))
423 // << " " << CurrentTh->Number(edge.EdgeVertex(1))
424 // << endl;
425 if ((link + linkp) == 1)
426 { // a boundary edge
427 Vertex * st = edge.EdgeVertex(1);
428 I2 I=*st;
429 Icoor2 ll = Norme2_2 (C-I);
430 if (ll < l1) { // the other vertex is neart
431 s1=st;
432 l1=ll;
433 er = edge;
434 if(ll<l0) { // change of direction --
435 s1=s;
436 l1=l0;
437 s=st;
438 l0=ll;
439 t=tt;
440 edge=Adj(edge);
441 link=linkp;
442 er = edge;
443 }
444 }
445 }
446
447 linkp=link;
448 edge=Adj(Previous(edge));
449 tt = edge;
450 } while (t!=tt);
451
452 assert((Triangle *) er);
453 I2 A((I2)*er.EdgeVertex(0));
454 I2 B((I2)*er.EdgeVertex(1));
455 I2 AB=B-A,AC=C-A,CB=B-C;
456 double aa = (double) (AB,AC);
457 double bb = (double) (AB,CB);
458 // cout << " " << aa << " " << bb
459 // << " " << CurrentTh->Number(er.EdgeVertex(0))
460 // << " " << CurrentTh->Number(er.EdgeVertex(1)) ;
461 if (aa<0) a=1,b=0;
462 else if(bb<0) a=0,b=1;
463 else
464 {
465 a = bb/(aa+bb);
466 b = aa/(aa+bb);
467 }
468 }
469
470 // cout <<" return= " << CurrentTh->Number(er.EdgeVertex(0)) << " "
471 // << CurrentTh->Number(er.EdgeVertex(1)) << " " << a
472 // << " " << b <<" " << l0 << " " <<l1 <<endl;
473 return er;
474 }
475
476
477
MetricAt(const R2 & A) const478 Metric Triangles::MetricAt (const R2 & A) const
479 { //if ((vertices <= &v) && (vertices < v+nbv)) return v.m;
480 I2 a = toI2(A);
481 Icoor2 deta[3];
482 Triangle * t =FindTriangleContening(a,deta);
483 if (t->det <0) { // outside
484 double ba,bb;
485 TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ;
486 return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}
487 else { // inside
488 Real8 aa[3];
489 Real8 s = deta[0]+deta[1]+deta[2];
490 aa[0]=deta[0]/s;
491 aa[1]=deta[1]/s;
492 aa[2]=deta[2]/s;
493 return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);
494 }
495 }
496
497
SplitEdge(const Triangles & Bh,const R2 & A,const R2 & B,int nbegin)498 void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
499 const R2 &A,const R2 &B,int nbegin)
500 { // SplitEdge
501 // if(SHOW) cout << " splitedge " << A << B << " " << nbegin << endl;
502 Triangle *tbegin, *t;
503
504 Icoor2 deta[3], deti,detj;
505 Real8 ba[3];
506 int nbt =0,ifirst=-1,ilast;
507 int i0,i1,i2;
508 int ocut,i,j,k=-1;
509 // int OnAVertices =0;
510 Icoor2 dt[3];
511 I2 a = Bh.toI2(A) ,b= Bh.toI2(B);// compute the Icoor a,b
512 I2 vi,vj;
513 int iedge =-1;// not a edge
514
515 if(nbegin) {// optimisation
516 // we suppose knowing the starting triangle
517 t=tbegin=lIntTria[ilast=(Size-1)].t;
518 if (tbegin->det>=0)
519 ifirst = ilast;}
520 else {// not optimisation
521 init();
522 t=tbegin = Bh.FindTriangleContening(a,deta);
523 // if(SHOW) cout <<t << " " << Real8(deta[0])/t->det<< " " << Real8(deta[1])/t->det
524 // << " " << Real8(deta[2])/t->det << endl;
525 if( t->det>=0)
526 ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det);
527 else
528 {// find the nearest boundary edge of the vertex A
529 // find a edge or such normal projection a the edge IJ is on the edge
530 // <=> IJ.IA >=0 && IJ.AJ >=0
531 ilast=ifirst;
532 double ba,bb;
533 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
534 Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
535 NewItem(A,Metric(ba,v0,bb,v1));
536 t=edge;
537 // test if the point b is in the same side
538 if (det(v0.i,v1.i,b)>=0) {
539 //cout << " All the edge " << A << B << endl;
540 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
541 Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
542 NewItem(A,Metric(ba,v0,bb,v1));
543 return;
544 }
545 } // find the nearest boundary edge of the vertex A
546 } // end not optimisation
547 if (t->det<0) { // outside departure
548 while (t->det <0) { // intersection boundary edge and a,b,
549 k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
550 assert(k>=0);
551 ocut = OppositeEdge[k];
552 i=VerticesOfTriangularEdge[ocut][0];
553 j=VerticesOfTriangularEdge[ocut][1];
554 vi=(*t)[i];
555 vj=(*t)[j];
556 deti = bamg::det(a,b,vi);
557 detj = bamg::det(a,b,vj);
558 // if(SHOW) { penthickness(3);
559 // Move(vi);Line(vj);CurrentTh->inquire();penthickness(1);
560 // cout << Bh.Number(tbegin) << " " << Bh.Number(t) << " i= " << i <<" j= " << j << " k=" << k
561 // << " deti= " << deti << " detj= " << detj
562 // << " v = " << Bh.Number((*t)[i]) << (*t)[i].r << " " << Bh.Number((*t)[j]) << (*t)[j].r << endl;}
563 if (deti>0) // go to i direction on gamma
564 ocut = PreviousEdge[ocut];
565 else if (detj<=0) // go to j direction on gamma
566 ocut = NextEdge[ocut];
567 TriangleAdjacent tadj =t->Adj(ocut);
568 t = tadj;
569 iedge= tadj;
570 if (t == tbegin) { //
571 double ba,bb;
572 if (verbosity>7)
573 cout << " SplitEdge: All the edge " << A << B << nbegin << det(vi,vj,b)
574 << " deti= " << deti << " detj=" <<detj << endl;
575 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
576 Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
577 NewItem(A,Metric(ba,v0,bb,v1));
578 return;
579 /*
580 cerr << nbegin << det(vi,vj,b) << " deti= " << deti << " detj=" <<detj << endl;
581 cerr << "SplitEdge on boucle A" << A << " B = " << B << endl;
582
583 #ifdef DRAWING
584 reffecran();
585 Bh.Draw();
586 penthickness(5);
587 Move(A);
588 Line(B);
589 penthickness(1);
590
591 Bh.inquire();
592 penthickness(5);
593 Move(A);
594 Line(B);
595 penthickness(1);
596 Bh.inquire();
597 #endif
598 MeshError(997);*/
599 }
600 } // end while (t->det <0)
601 // theoriticaly we have: deti =<0 and detj>0
602
603 // computation of barycentric coor
604 // test if the point b is on size on t
605 // we revert vi,vj because vi,vj is def in Adj triangle
606 if ( det(vi,vj,b)>=0) {
607 if (verbosity>7)
608 cout << " SplitEdge: all AB outside " << A << B << endl;
609 t=tbegin;
610 Real8 ba,bb;
611 TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);
612 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
613 return;
614 }
615 else
616 {
617 k = OppositeVertex[iedge];
618 i=VerticesOfTriangularEdge[iedge][0];
619 j=VerticesOfTriangularEdge[iedge][1];
620 Real8 dij = detj-deti;
621 assert(i+j+k == 0 + 1 +2);
622 ba[j] = detj/dij;
623 ba[i] = -deti/dij;
624 ba[k] = 0;
625 // if(SHOW) cout << i << " " << j << " " << k << " " << ba[i] << " " << ba[j] << endl;
626 ilast=NewItem(t,ba[0],ba[1],ba[2]); }
627 } // outside departure
628
629
630
631 // recherche the intersection of [a,b] with Bh Mesh.
632 // we know a triangle ta contening the vertex a
633 // we have 2 case for intersection [a,b] with a edge [A,B] of Bh
634 // 1) the intersection point is in ]A,B[
635 // 2) is A or B
636 // first version ---
637 for (;;) {
638 // t->Draw();
639 if (iedge < 0) {
640 i0 =0;i1=1;i2=2;
641 dt[0] =bamg::det(a,b,(*t)[0]);
642 dt[1] =bamg::det(a,b,(*t)[1]);
643 dt[2] =bamg::det(a,b,(*t)[2]);}
644 else {
645 i2 = iedge;
646 i0 = NextEdge[i2];
647 i1 = NextEdge[i0];
648 dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because
649 dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side
650 dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);}
651
652 // so we have just to see the transition from - to + of the det0..2 on edge of t
653 // because we are going from a to b
654 if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
655 ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))
656 ocut =i0;
657 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
658 (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
659 ocut =i1;
660 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
661 (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
662 ocut =i2;
663 else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) &&
664 ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))
665 ocut =i0;
666 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) &&
667 (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
668 ocut =i1;
669 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) &&
670 (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
671 ocut =i2;
672 else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
673 ( dt[j=VerticesOfTriangularEdge[i0][1]] == 0))
674 ocut =i0;
675 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
676 (dt[j=VerticesOfTriangularEdge[i1][1]] == 0))
677 ocut =i1;
678 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
679 (dt[j=VerticesOfTriangularEdge[i2][1]] == 0))
680 ocut =i2;
681 else { // On a edge (2 zero)
682 k =0;
683 if (dt[0]) ocut=0,k++;
684 if (dt[1]) ocut=1,k++;
685 if (dt[2]) ocut=2,k++;
686 if(k == 1) {
687 if (dt[ocut] >0) // triangle upper AB
688 ocut = NextEdge[ocut];
689 i= VerticesOfTriangularEdge[ocut][0];
690 j= VerticesOfTriangularEdge[ocut][1];
691 }
692 else {
693 cerr << " Bug Split Edge " << endl;
694 cerr << " dt[0]= " << dt[0]
695 << " dt[1]= " << dt[1]
696 << " dt[2]= "<< dt[2] << endl;
697 cerr << i0 << " " << i1 << " " << i2 << endl;
698 cerr << " A = " << A << " B= " << B << endl;
699 cerr << " Triangle t = " << *t << endl;
700 cerr << (*t)[0] << (*t)[1] << (*t)[0] << endl;
701 cerr << " nbt = " << nbt << endl;
702 MeshError(100);}}
703
704 k = OppositeVertex[ocut];
705
706 Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);
707
708
709 if (detbij >= 0) { //we find the triangle contening b
710 dt[0]=bamg::det((*t)[1],(*t)[2],b);
711 dt[1]=bamg::det((*t)[2],(*t)[0],b);
712 dt[2]=bamg::det((*t)[0],(*t)[1],b);
713 #ifdef DEBUG
714 assert(dt[0] >= 0);
715 assert(dt[1] >= 0);
716 assert(dt[2] >= 0);
717 #endif
718 Real8 dd = t->det;
719 NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);
720 return ;}
721 else { // next triangle by adjacent by edge ocut
722 deti = dt[i];
723 detj = dt[j];
724 Real4 dij = detj-deti;
725 ba[i] = detj/dij;
726 ba[j] = -deti/dij;
727 ba[3-i-j ] = 0;
728 ilast=NewItem(t, ba[0],ba[1],ba[2]);
729
730 TriangleAdjacent ta =t->Adj(ocut);
731 t = ta;
732 iedge= ta;
733 if (t->det <= 0) {
734 double ba,bb;
735 TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);
736 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
737 // cout << " return " << ba << " " << bb << endl;
738 // ajoute le 03 frev 1997 par F. hecht
739 return;
740 }
741 }// we go outside of omega
742 } // for(;;)
743
744
745 } // routine SplitEdge
746
747
NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2)748 int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2) {
749 register int n;
750 R2 x(0,0);
751 if ( d0) x = (*tt)[0].r * d0;
752 if ( d1) x = x + (*tt)[1].r * d1;
753 if ( d2) x = x + (*tt)[2].r * d2;
754 // newer add same point
755 if(!Size || Norme2_2(lIntTria[Size-1].x-x)) {
756 if (Size==MaxSize) ReShape();
757 lIntTria[Size].t=tt;
758 lIntTria[Size].bary[0]=d0;
759 lIntTria[Size].bary[1]=d1;
760 lIntTria[Size].bary[2]=d2;
761 lIntTria[Size].x = x;
762 Metric m0,m1,m2;
763 register Vertex * v;
764 if ((v=(*tt)(0))) m0 = v->m;
765 if ((v=(*tt)(1))) m1 = v->m;
766 if ((v=(*tt)(2))) m2 = v->m;
767 lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2);
768 #ifdef DEBUG1
769 if(SHOW) { cout << "SHOW ++ NewItem =" << Size << x ;
770 cout << " " << d0 << " " << d1 << " " << d2 <<endl;}
771 #endif
772 n=Size++;}
773 else n=Size-1;
774 return n;
775 }
NewItem(R2 A,const Metric & mm)776 int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
777 register int n;
778 if(!Size || Norme2_2(lIntTria[Size-1].x-A)) {
779 if (Size==MaxSize) ReShape();
780 lIntTria[Size].t=0;
781 lIntTria[Size].x=A;
782 lIntTria[Size].m=mm;
783 #ifdef DEBUG1
784 if (SHOW) cout << "SHOW ++ NewItem A" << Size << A << endl;
785 #endif
786 n=Size++;
787 }
788 else n=Size-1;
789 return n;
790 }
791
Length()792 Real8 ListofIntersectionTriangles::Length()
793 {
794 // cout << " n= " << Size << ":" ;
795 assert(Size>0);
796 // computation of the length
797 R2 C;
798 Metric Mx,My;
799 int ii,jj;
800 R2 x,y,xy;
801
802 SegInterpolation *SegI=lSegsI;
803 SegI=lSegsI;
804 lSegsI[NbSeg].last=Size;// improvement
805
806 int EndSeg=Size;
807
808 y = lIntTria[0].x;
809 Real8 sxy, s = 0;
810 lIntTria[0].s =0;
811 SegI->lBegin=s;
812
813 for (jj=0,ii=1;ii<Size;jj=ii++)
814 {
815 // seg jj,ii
816 x=y;
817 y = lIntTria[ii].x;
818 xy = y-x;
819 Mx = lIntTria[ii].m;
820 My = lIntTria[jj].m;
821 // Real8 &sx= lIntTria[ii].sp; // previous seg
822 // Real8 &sy= lIntTria[jj].sn; // next seg
823 // sx = Mx(xy);
824 // sy = My(xy);
825 // sxy = (Mx(xy)+ My(xy))/2.0;
826 sxy = LengthInterpole(Mx,My,xy);
827 s += sxy;
828 lIntTria[ii].s = s;
829 if (ii == EndSeg)
830 SegI->lEnd=s,
831 SegI++,
832 EndSeg=SegI->last,
833 SegI->lBegin=s;
834
835 // cout << ii << " " << jj << x<< y <<xy << s << lIntTria[ii].m ;
836 }
837 len = s;
838 SegI->lEnd=s;
839
840 // cout << " len= " << s << endl;
841 return s;
842 }
843
NewPoints(Vertex * vertices,Int4 & nbv,Int4 nbvx)844 Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4 nbvx)
845 {
846
847 const Int4 nbvold = nbv;
848 Real8 s = Length();
849 if (s < 1.5 ) return 0;
850 //////////////////////
851 int ii = 1 ;
852 R2 y,x;
853 Metric My,Mx ;
854 Real8 sx =0,sy;
855 int nbi = Max(2,(int) (s+0.5));
856 Real8 sint = s/nbi;
857 Real8 si = sint;
858
859 int EndSeg=Size;
860 SegInterpolation *SegI=0;
861 if (NbSeg)
862 SegI=lSegsI,EndSeg=SegI->last;
863
864 for (int k=1;k<nbi;k++)
865 {
866 while ((ii < Size) && ( lIntTria[ii].s <= si ))
867 if (ii++ == EndSeg)
868 SegI++,EndSeg=SegI->last;
869
870 int ii1=ii-1;
871 x =lIntTria[ii1].x;
872 sx =lIntTria[ii1].s;
873 Metric Mx=lIntTria[ii1].m;
874 #ifdef DEBUG
875 double lx = lIntTria[ii-1].sn;
876 #endif
877 y =lIntTria[ii].x;
878 sy =lIntTria[ii].s;
879 Metric My=lIntTria[ii].m;
880 #ifdef DEBUG
881 double ly =lIntTria[ii].sp;
882 assert( sx <= si);
883 assert( si <= sy);
884 assert( sy != sx);
885 #endif
886
887 Real8 lxy = sy-sx;
888 Real8 cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
889
890 R2 C;
891 Real8 cx = 1-cy;
892 C = SegI ? SegI->F(si): x * cx + y *cy;
893
894 si += sint;
895 if ( nbv<nbvx) {
896 vertices[nbv].r = C;
897 vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
898 if((verbosity/100%10)==2)
899 cout << " -- Add point " << nbv-1 << " " << vertices[nbv-1] << " " << vertices[nbv-1].m << endl;
900
901 #ifdef DEBUG
902 if(k>1) {
903 R2 AB = vertices[nbv-2].r - vertices[nbv-1].r ;
904 Real8 dp = LengthInterpole(vertices[nbv-2].m,vertices[nbv-1].m,AB);
905 if (dp > 1.6) {
906 cerr << "PB calcul new Int. points trop loin l=" << dp << " v=" << nbv-1 << " " << nbv-2 <<Mx<<My<<y-x << endl;
907 }
908 }
909 #endif
910 }
911 else return nbv-nbvold;
912 }
913 return nbv-nbvold;
914 }
915
SwapForForcingEdge(Vertex * & pva,Vertex * & pvb,TriangleAdjacent & tt1,Icoor2 & dets1,Icoor2 & detsa,Icoor2 & detsb,int & NbSwap)916 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
917 TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap)
918 { // l'arete ta coupe l'arete pva pvb
919 // de cas apres le swap sa coupe toujours
920 // on cherche l'arete suivante
921 // on suppose que detsa >0 et detsb <0
922 // attention la routine echange pva et pvb
923
924 if(tt1.Locked()) return 0; // frontiere croise
925
926 TriangleAdjacent tt2 = Adj(tt1);
927 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent
928 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
929 assert ( a1 >= 0 && a1 < 3 );
930
931 Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
932 Vertex & s1= (*t1)[OppositeVertex[a1]];
933 Vertex & s2= (*t2)[OppositeVertex[a2]];
934
935
936 Icoor2 dets2 = det(*pva,*pvb,s2);
937
938 #ifdef DEBUG
939 Vertex & sb= (*t1)[VerticesOfTriangularEdge[a1][1]];
940 Icoor2 wdets1 = det(*pva,*pvb,s1);
941 Icoor2 wdetsa = det(*pva,*pvb,sa);
942 Icoor2 wdetsb = det(*pva,*pvb,sb);
943 assert(wdets1 == dets1);
944 assert(wdetsa == detsa);
945 assert(wdetsb == detsb);
946 #endif
947
948 Icoor2 det1=t1->det , det2=t2->det ;
949 #ifdef DEBUG
950 assert(det1>0 && det2 >0);
951 Icoor2 ddet1 = det((*t1)[0],(*t1)[1],(*t1)[2]);
952 Icoor2 ddet2 = det((*t2)[0],(*t2)[1],(*t2)[2]);
953 if ((det1 != ddet1) || (det2 != ddet2) )
954 {
955 assert(det1 == ddet1);
956 assert(det2 == ddet2);
957 }
958 Icoor2 detvasasb = det(*pva,sa,sb);
959 Icoor2 detvbsasb = det(*pvb,sa,sb);
960 if ( CurrentTh && ! ( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0))))
961 {
962 cout << " detvasasb =" << detvasasb << "detvbsasb = " << detvbsasb
963 << " " << pva << " " << pvb << " " <<CurrentTh <<endl;
964 #ifdef DRAWING1
965 reffecran();
966 CurrentTh->Draw();
967 penthickness(10);
968 pva->MoveTo();pvb->LineTo();
969 penthickness(1);
970 CurrentTh->inquire();
971 #endif
972 }
973 assert( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0)));
974 #endif
975
976 Icoor2 detT = det1+det2;
977 assert((det1>0 ) && (det2 > 0));
978 assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb
979 Icoor2 ndet1 = bamg::det(s1,sa,s2);
980 Icoor2 ndet2 = detT - ndet1;
981
982 int ToSwap =0; //pas de swap
983 if ((ndet1 >0) && (ndet2 >0))
984 { // on peut swaper
985 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))
986 ToSwap =1;
987 else // swap alleatoire
988 if (BinaryRand())
989 ToSwap =2;
990 }
991 #ifdef DEBUG
992 if (ForDebugging) {
993 cerr << "swap = " << ToSwap << " ndet1 " << ndet1 << ", ndet2 " << ndet2 << "det1 " << det1 << " det2 " << det2
994 << " if1 = " << ((ndet1 >0) && (ndet2 >0))
995 << " if2 = " << ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) << endl;
996 #ifdef DRAWING
997 couleur(0);
998 t1->Draw();
999 t2->Draw();
1000 #endif
1001 }
1002 #endif
1003 if (ToSwap) NbSwap++,
1004 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);
1005
1006 #ifdef DEBUG
1007 if (ForDebugging) {
1008 #ifdef DRAWING
1009 couleur(4);
1010 t1->Draw();
1011 t2->Draw();
1012 rattente(1);
1013 #endif
1014 }
1015 #endif
1016 int ret=1;
1017
1018 if (dets2 < 0) {// haut
1019 dets1 = ToSwap ? dets1 : detsa ;
1020 detsa = dets2;
1021 tt1 = Previous(tt2) ;}
1022 else if (dets2 > 0){// bas
1023 dets1 = ToSwap ? dets1 : detsb ;
1024 detsb = dets2;
1025 //xxxx tt1 = ToSwap ? tt1 : Next(tt2);
1026 if(!ToSwap) tt1 = Next(tt2);
1027 }
1028 else { // changement de sens
1029 if (ForDebugging) cout << "changement de sens" << endl;
1030 ret = -1;
1031 Exchange(pva,pvb);
1032 Exchange(detsa,detsb);
1033 Exchange(dets1,dets2);
1034 Exchange(tt1,tt2);
1035 dets1=-dets1;
1036 dets2=-dets2;
1037 detsa=-detsa;
1038 detsb=-detsb;
1039
1040 if (ToSwap)
1041 if (dets2 < 0) {// haut
1042 dets1 = (ToSwap ? dets1 : detsa) ;
1043 detsa = dets2;
1044 tt1 = Previous(tt2) ;}
1045 else if (dets2 > 0){// bas
1046 dets1 = (ToSwap ? dets1 : detsb) ;
1047 detsb = dets2;
1048 if(!ToSwap) tt1 = Next(tt2);
1049 }
1050 else {// on a fin ???
1051 tt1 = Next(tt2);
1052 ret =0;}
1053
1054 }
1055 return ret;
1056 }
1057
ForceEdge(Vertex & a,Vertex & b,TriangleAdjacent & taret)1058 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
1059 {
1060 #ifdef DEBUG
1061 restart: // for debug
1062 #endif
1063 int NbSwap =0;
1064 assert(a.t && b.t); // the 2 vertex is in a mesh
1065 int k=0;
1066 taret=TriangleAdjacent(0,0); // erreur
1067
1068 TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]);
1069 Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
1070 // we turn around a in the direct sens
1071
1072 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1;
1073 if(v2) // normal case
1074 det2 = det(*v2,a,b);
1075 else { // no chance infini vertex try the next
1076 tta= Previous(Adj(tta));
1077 v2 = tta.EdgeVertex(0);
1078 vbegin =v2;
1079 assert(v2);
1080 det2 = det(*v2,a,b);
1081 // cout << " No Change try the next" << endl;
1082 }
1083
1084 #ifdef DRAWING1
1085 a.MoveTo();b.LineTo();
1086 #endif
1087
1088 while (v2 != &b) {
1089 TriangleAdjacent tc = Previous(Adj(tta));
1090 v1 = v2;
1091 v2 = tc.EdgeVertex(0);
1092 det1 = det2;
1093 #ifdef DEBUG
1094 assert( v1 == tta.EdgeVertex(0));
1095 assert( &a == tc.EdgeVertex(1) );
1096 #endif
1097 det2 = v2 ? det(*v2,a,b): det2;
1098
1099 if((det1 < 0) && (det2 >0)) {
1100 // try to force the edge
1101 Vertex * va = &a, *vb = &b;
1102 tc = Previous(tc);
1103 assert ( v1 && v2);
1104 Icoor2 detss = 0,l=0,ks;
1105 // cout << "Real ForcingEdge " << *va << *vb << detss << endl;
1106 #ifdef DEBUG
1107 Icoor2 dettt1 = det(*v1,a,b);
1108 Icoor2 dettt2 = det(*v2,a,b);
1109
1110 if (!(dettt1==det1 && dettt2==det2))
1111 {
1112 assert(ForDebugging==0);
1113 ForDebugging=1;
1114 goto restart;
1115 }
1116
1117 #endif
1118 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap)))
1119 if(l++ > 10000000) {
1120 cerr << " Loop in forcing Egde AB"
1121 <<"\n vertex A " << a
1122 <<"\n vertex B " << b
1123 <<"\n nb de swap " << NbSwap
1124 <<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl;
1125
1126 if ( CurrentTh )
1127 cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl;
1128 #ifdef DEBUG
1129 ForDebugging = 1;
1130 #endif
1131 #ifdef DRAWING1
1132 if ( CurrentTh ) {
1133 reffecran();
1134 couleur(6);
1135 CurrentTh->Draw();
1136 couleur(1);
1137 penthickness(10);
1138 a.MoveTo();b.LineTo();
1139 penthickness(1);
1140 CurrentTh->inquire();
1141 couleur(6);
1142 l=0;
1143 reffecran();
1144 while (ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap) && (l++ < 1000))
1145 cerr << " " << CurrentTh->Number(tc.EdgeVertex(0))<<" " <<CurrentTh->Number(tc.EdgeVertex(1)) << " ";
1146 }
1147 #endif
1148 MeshError(990);
1149 }
1150 Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
1151 if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) {
1152 tc.SetLock();
1153 a.Optim(1,0);
1154 b.Optim(1,0);
1155 taret = tc;
1156 return NbSwap;
1157 }
1158 else
1159 {
1160 taret = tc;
1161 return -2; // error boundary is crossing
1162 /* cerr << "Fatal Error boundary is crossing ";
1163 if(CurrentTh)
1164 {
1165 cerr << " edge: [" << CurrentTh->Number(a) << ", " << CurrentTh->Number(b) << " ] and [ ";
1166 cerr << CurrentTh->Number(aa) << " " << CurrentTh->Number(bb) << " ] " << endl;
1167 }
1168 MeshError(991);
1169 */
1170 }
1171 }
1172 tta = tc;
1173 assert(k++<2000);
1174 if ( vbegin == v2 ) return -1;// error
1175 }
1176
1177 tta.SetLock();
1178 taret=tta;
1179 a.Optim(1,0);
1180 b.Optim(1,0);
1181 return NbSwap;
1182 }
1183
1184
swap(Int2 a,int koption)1185 int Triangle::swap(Int2 a,int koption){
1186 #ifdef DEBUG
1187 if(a &4 ) return 0;// arete lock
1188 int munswap1 = a/4;
1189 a &=3;
1190 #else
1191 if(a/4 !=0) return 0;// arete lock or MarkUnSwap
1192 #endif
1193
1194 register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent
1195 register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles
1196 #ifdef DEBUG
1197 if(a2 & 4) return 0; // arete lock
1198 int munswap2 = a2/4;
1199 a2 &= 3;
1200 #else
1201 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
1202 #endif
1203
1204 register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]];
1205 register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]];
1206 register Vertex *s1=t1->ns[OppositeVertex[a1]];
1207 register Vertex *s2=t2->ns[OppositeVertex[a2]];
1208
1209 #ifdef DEBUG
1210 assert ( a >= 0 && a < 3 );
1211 #endif
1212
1213 Icoor2 det1=t1->det , det2=t2->det ;
1214 Icoor2 detT = det1+det2;
1215 Icoor2 detA = Abs(det1) + Abs(det2);
1216 Icoor2 detMin = Min(det1,det2);
1217
1218 int OnSwap = 0;
1219 // si 2 triangle infini (bord) => detT = -2;
1220 if (sa == 0) {// les deux triangles sont frontieres
1221 det2=bamg::det(s2->i,sb->i,s1->i);
1222 OnSwap = det2 >0;}
1223 else if (sb == 0) { // les deux triangles sont frontieres
1224 det1=bamg::det(s1->i,sa->i,s2->i);
1225 OnSwap = det1 >0;}
1226 else if(( s1 != 0) && (s2 != 0) ) {
1227 det1 = bamg::det(s1->i,sa->i,s2->i);
1228 det2 = detT - det1;
1229 OnSwap = (Abs(det1) + Abs(det2)) < detA;
1230
1231 Icoor2 detMinNew=Min(det1,det2);
1232 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test
1233 if (! OnSwap &&(detMinNew>0)) {
1234 OnSwap = detMin ==0;
1235 if (! OnSwap) {
1236 int kopt = koption;
1237 while (1)
1238 if(kopt) {
1239 // critere de Delaunay pure isotrope
1240 register Icoor2 xb1 = sb->i.x - s1->i.x,
1241 x21 = s2->i.x - s1->i.x,
1242 yb1 = sb->i.y - s1->i.y,
1243 y21 = s2->i.y - s1->i.y,
1244 xba = sb->i.x - sa->i.x,
1245 x2a = s2->i.x - sa->i.x,
1246 yba = sb->i.y - sa->i.y,
1247 y2a = s2->i.y - sa->i.y;
1248 register double
1249 cosb12 = double(xb1*x21 + yb1*y21),
1250 cosba2 = double(xba*x2a + yba*y2a) ,
1251 sinb12 = double(det2),
1252 sinba2 = double(t2->det);
1253
1254
1255 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2)
1256 OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12);
1257 // if(CurrentTh)
1258 // cout << "swap " << CurrentTh->Number(sa) << " " << CurrentTh->Number(sb) << " " ;
1259 // cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb12
1260 // << " Onswap = " << OnSwap << endl;
1261 break;
1262 }
1263 else
1264 {
1265 // critere de Delaunay anisotrope
1266 Real8 som;
1267 I2 AB=(I2) *sb - (I2) *sa;
1268 I2 MAB2=((I2) *sb + (I2) *sa);
1269 R2 MAB(MAB2.x*0.5,MAB2.y*0.5);
1270 I2 A1=(I2) *s1 - (I2) *sa;
1271 I2 D = (I2) * s1 - (I2) * sb ;
1272 R2 S2(s2->i.x,s2->i.y);
1273 R2 S1(s1->i.x,s1->i.y);
1274 {
1275 Metric M=s1->m;
1276 R2 ABo = M.Orthogonal(AB);
1277 R2 A1o = M.Orthogonal(A1);
1278 // (A+B)+ x ABo = (S1+B)/2+ y A1
1279 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
1280 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
1281 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
1282 if (Abs(d) > dd*1.e-3) {
1283 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
1284 som = M(C - S2)/M(C - S1);
1285 } else
1286 {kopt=1;continue;}
1287
1288 }
1289 {
1290 Metric M=s2->m;
1291 R2 ABo = M.Orthogonal(AB);
1292 R2 A1o = M.Orthogonal(A1);
1293 // (A+B)+ x ABo = (S1+B)/2+ y A1
1294 // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
1295 double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
1296 double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
1297 if(Abs(d) > dd*1.e-3) {
1298 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
1299 som += M(C - S2)/M(C - S1);
1300 } else
1301 {kopt=1;continue;}
1302 }
1303 OnSwap = som < 2;
1304 break;
1305 }
1306
1307 } // OnSwap
1308 } // (! OnSwap &&(det1 > 0) && (det2 > 0) )
1309 }
1310 #ifdef DEBUG1
1311 if (OnSwap && ( munswap1 || munswap2)) {
1312 cout << " erreur Mark unswap T " << CurrentTh->Number(t1) << " " << CurrentTh->Number(t2) << endl
1313 << *t1 << endl
1314 << *t2 << endl;
1315 return 0;
1316 }
1317 #endif
1318 if( OnSwap )
1319 bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2);
1320 else {
1321 NbUnSwap ++;
1322 t1->SetMarkUnSwap(a1);
1323 }
1324 return OnSwap;
1325 }
1326
Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart,Real8 omega)1327 Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart ,Real8 omega)
1328 {
1329 #ifdef DEBUG
1330 register Int4 NbSwap =0;
1331 #endif
1332 register Vertex * s = this;
1333 Vertex &vP = *s,vPsave=vP;
1334 // if (vP.on) return 0;// Don't move boundary vertex
1335
1336 register Triangle * tbegin= t , *tria = t , *ttc;
1337
1338 register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc;
1339 R2 P(s->r),PNew(0,0);
1340 // cout << BTh.quadtree << " " << BTh.quadtree->root << endl;
1341 // assert(BTh.quadtree && BTh.quadtree->root);
1342 do {
1343 k++;
1344
1345 #ifdef DEBUG
1346 assert( s == & (*tria)[VerticesOfTriangularEdge[j][1]] );
1347 assert( tria->det >0);
1348 #endif
1349 if (!tria->Hidden(j))
1350 {
1351 Vertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];
1352
1353 R2 Q = vQ,QP(P-Q);
1354 Real8 lQP = LengthInterpole(vP,vQ,QP);
1355 PNew += Q+QP/Max(lQP,1e-20);
1356 kk ++;
1357 }
1358 ttc = tria->TriangleAdj(j);
1359 jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
1360 tria = ttc;
1361 j = NextEdge[jc];
1362 assert(k<2000);
1363 } while ( tbegin != tria);
1364 if (kk<4) return 0;
1365 PNew = PNew/(Real8)kk;
1366 R2 Xmove((PNew-P)*omega);
1367 PNew = P+Xmove;
1368 Real8 delta=Norme2_2(Xmove);
1369
1370
1371 //
1372 Icoor2 deta[3];
1373 I2 IBTh = BTh.toI2(PNew);
1374
1375 tstart=BTh.FindTriangleContening(IBTh,deta,tstart);
1376
1377 if (tstart->det <0)
1378 { // outside
1379 double ba,bb;
1380 TriangleAdjacent edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ;
1381 tstart = edge;
1382 vP.m= Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));
1383 }
1384 else
1385 { // inside
1386 Real8 aa[3];
1387 Real8 s = deta[0]+deta[1]+deta[2];
1388 aa[0]=deta[0]/s;
1389 aa[1]=deta[1]/s;
1390 aa[2]=deta[2]/s;
1391 vP.m = Metric(aa,(*tstart)[0],(*tstart)[1],(*tstart)[2]);
1392 }
1393
1394 // recompute the det of the triangle
1395 vP.r = PNew;
1396
1397 vP.i = Th.toI2(PNew);
1398
1399 Vertex vPnew = vP;
1400
1401 int ok=1;
1402 int loop=1;
1403 k=0;
1404 while (ok)
1405 {
1406 ok =0;
1407 do {
1408 k++;
1409 double detold = tria->det;
1410 tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]);
1411 if (loop)
1412 {
1413 Vertex *v0,*v1,*v2,*v3;
1414 if (tria->det<0) ok =1;
1415 else if (tria->Quadrangle(v0,v1,v2,v3))
1416 {
1417 vP = vPsave;
1418 Real8 qold =QuadQuality(*v0,*v1,*v2,*v3);
1419 vP = vPnew;
1420 Real8 qnew = QuadQuality(*v0,*v1,*v2,*v3);
1421 if (qnew<qold) ok = 1;
1422 }
1423 else if ( (double)tria->det < detold/2 ) ok=1;
1424
1425 }
1426 tria->SetUnMarkUnSwap(0);
1427 tria->SetUnMarkUnSwap(1);
1428 tria->SetUnMarkUnSwap(2);
1429 ttc = tria->TriangleAdj(j);
1430 jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
1431 tria = ttc;
1432 j = NextEdge[jc];
1433 assert(k<2000);
1434 } while ( tbegin != tria);
1435 if (ok && loop) vP=vPsave; // no move
1436 loop=0;
1437 }
1438 return delta;
1439 }
1440
1441
1442
1443
Add(Vertex & s,Triangle * t,Icoor2 * det3)1444 void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
1445 {
1446 // -------------------------------------------
1447 // s2
1448 // !
1449 // /|\ !
1450 // / | \ !
1451 // / | \ !
1452 // tt1 / | \ tt0 !
1453 // / |s \ !
1454 // / . \ !
1455 // / . ` \ !
1456 // / . ` \ !
1457 // ---------------- !
1458 // s0 tt2 s1
1459 //--------------------------------------------
1460
1461 Triangle * tt[3]; // the 3 new Triangles
1462 Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2];
1463 Icoor2 det3local[3];
1464 int infv = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0;
1465 // infv = ordre of the infini vertex (null)
1466 register int nbd0 =0; // number of zero det3
1467 register int izerodet=-1,iedge; // izerodet = egde contening the vertex s
1468 Icoor2 detOld = t->det;
1469
1470 if ( ( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) )
1471 {
1472 cerr << " infv " << infv << " det = " << detOld << endl;
1473 cerr << Number(s) << " "<< Number(s0) << " "
1474 << Number(s1) << " " << Number(s2) << endl;
1475 MeshError(3);
1476 }
1477
1478 // if det3 do not exist then constuct det3
1479 if (!det3) {
1480 det3 = det3local; // alloc
1481 if ( infv<0 ) {
1482 det3[0]=bamg::det(s ,s1,s2);
1483 det3[1]=bamg::det(s0,s ,s2);
1484 det3[2]=bamg::det(s0,s1,s );}
1485 else {
1486 // one of &s1 &s2 &s0 is NULL so (&si || &sj) <=> !&sk
1487 det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ;
1488 det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ;
1489 det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}}
1490
1491
1492 if (!det3[0]) izerodet=0,nbd0++;
1493 if (!det3[1]) izerodet=1,nbd0++;
1494 if (!det3[2]) izerodet=2,nbd0++;
1495
1496 if (nbd0 >0 ) // point s on a egde or on a vertex
1497 if (nbd0 == 1) {
1498 iedge = OppositeEdge[izerodet];
1499 TriangleAdjacent ta = t->Adj(iedge);
1500
1501 #ifdef DEBUG1
1502 cout << " the point " << Number(s) << " is the edge " << izerodet
1503 << " of " << Number(t) << " det3 = "
1504 << det3[0] << " " << det3[1] << " " << det3[2] << " " << endl;
1505 cout << " ta = " << ta << "ta->det =" << ((Triangle*) ta)->det
1506 << " "<< t->det<< endl;
1507 #endif
1508
1509 // the point is on the edge
1510 // if the point is one the boundary
1511 // add the point in outside part
1512 if ( t->det >=0) { // inside triangle
1513 if ((( Triangle *) ta)->det < 0 ) {
1514 // add in outside triangle
1515 Add(s,( Triangle *) ta);
1516 return;}
1517 }}
1518 else {
1519 cerr << " bug " << nbd0 <<endl;
1520 cerr << " Bug double points in " << endl ;
1521 cerr << " s = " << Number(s) << " " << s << endl;
1522 cerr << " s0 = "<< Number(s0) << " " << s0 << endl;
1523 cerr << " s1 = "<< Number(s1) << " " << s1 << endl;
1524 cerr << " s2 = "<< Number(s2) << " " << s2 << endl;
1525 MeshError(5,this);}
1526
1527 // remove de MarkUnSwap edge
1528 t->SetUnMarkUnSwap(0);
1529 t->SetUnMarkUnSwap(1);
1530 t->SetUnMarkUnSwap(2);
1531
1532 tt[0]= t;
1533 tt[1]= &triangles[nbt++];
1534 tt[2]= &triangles[nbt++];
1535
1536 if (nbt>nbtx) {
1537 cerr << " No enougth triangles " << endl;
1538 MeshError(999,this);
1539 }
1540
1541 *tt[1]= *tt[2]= *t;
1542 // gestion of the link
1543 tt[0]->link=tt[1];
1544 tt[1]->link=tt[2];
1545
1546 (* tt[0])(OppositeVertex[0])=&s;
1547 (* tt[1])(OppositeVertex[1])=&s;
1548 (* tt[2])(OppositeVertex[2])=&s;
1549
1550 tt[0]->det=det3[0];
1551 tt[1]->det=det3[1];
1552 tt[2]->det=det3[2];
1553
1554 // update adj des triangles externe
1555 tt[0]->SetAdjAdj(0);
1556 tt[1]->SetAdjAdj(1);
1557 tt[2]->SetAdjAdj(2);
1558 // update des adj des 3 triangle interne
1559 const int i0 = 0;
1560 const int i1= NextEdge[i0];
1561 const int i2 = PreviousEdge[i0];
1562
1563 tt[i0]->SetAdj2(i2,tt[i2],i0);
1564 tt[i1]->SetAdj2(i0,tt[i0],i1);
1565 tt[i2]->SetAdj2(i1,tt[i1],i2);
1566
1567 tt[0]->SetTriangleContainingTheVertex();
1568 tt[1]->SetTriangleContainingTheVertex();
1569 tt[2]->SetTriangleContainingTheVertex();
1570
1571
1572 // swap if the point s is on a edge
1573 if(izerodet>=0) {
1574 // cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl;
1575 int rswap =tt[izerodet]->swap(iedge);
1576
1577 if (!rswap)
1578 {
1579 cout << " Pb swap the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl;
1580 #ifdef DRAWING
1581 if( CurrentTh && withrgraphique)
1582 {
1583 reffecran();
1584
1585 DrawMark(s.r);
1586 CurrentTh->inquire();
1587 DrawMark(s.r);
1588 rattente(1);
1589 }
1590 #endif
1591 }
1592 assert(rswap);
1593 }
1594
1595 #ifdef DEBUG
1596 tt[0]->check();
1597 tt[1]->check();
1598 tt[2]->check();
1599 #endif
1600 #ifdef DRAWING1
1601 tt[0]->Draw();
1602 tt[1]->Draw();
1603 tt[2]->Draw();
1604 #endif
1605
1606 }
1607
1608
SplitInternalEdgeWithBorderVertices()1609 Int4 Triangles::SplitInternalEdgeWithBorderVertices()
1610 {
1611 Int4 NbSplitEdge=0;
1612 SetVertexFieldOn();
1613 Int4 it;
1614 Int4 nbvold=nbv;
1615 for (it=0;it<nbt;it++)
1616 {
1617 Triangle &t=triangles[it];
1618 if (t.link)
1619 for (int j=0;j<3;j++)
1620 if(!t.Locked(j) && !t.Hidden(j)){
1621 Triangle &tt = *t.TriangleAdj(j);
1622 if ( &tt && tt.link && it < Number(tt))
1623 { // an internal edge
1624 Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];
1625 Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];
1626 if (v0.on && v1.on)
1627 {
1628 R2 P= ((R2) v0 + (R2) v1)*0.5;
1629 if ( nbv<nbvx) {
1630 vertices[nbv].r = P;
1631 vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);
1632 vertices[nbv].ReferenceNumber=0;
1633 vertices[nbv].DirOfSearch = NoDirOfSearch ;
1634 }
1635 NbSplitEdge++;
1636 if (verbosity>7)
1637 cout <<" Internal edge with two vertices on boundary"
1638 << Number(v0) << " " << Number(v1) << " by " << endl;
1639 }
1640 }
1641 }
1642 }
1643 ReMakeTriangleContainingTheVertex();
1644 if (nbvold!=nbv)
1645 {
1646 Int4 iv = nbvold;
1647 Int4 NbSwap = 0;
1648 Icoor2 dete[3];
1649 for (Int4 i=nbvold;i<nbv;i++)
1650 {// for all the new point
1651 Vertex & vi = vertices[i];
1652 vi.i = toI2(vi.r);
1653 vi.r = toR2(vi.i);
1654 // if (!quadtree->ToClose(vi,seuil,hi,hj)) {
1655 // a good new point
1656 vi.ReferenceNumber=0;
1657 vi.DirOfSearch =NoDirOfSearch;
1658 // cout << " Add " << Number(vi) << " " << vi
1659 // << " " << Number(vi) << " <--> " << Number(vi) <<endl;
1660 Triangle *tcvi = FindTriangleContening(vi.i,dete);
1661 if (tcvi && !tcvi->link) {
1662 cout << i << " PB insert point " << Number(vi) << vi << Number(vi)
1663 << " tcvi = " << tcvi << " " << tcvi->link << endl;
1664 cout << (*tcvi)[1] << (*tcvi)[2] << endl;
1665 tcvi = FindTriangleContening(vi.i,dete);
1666 cout << (*tcvi)[1] << (*tcvi)[2] << endl;
1667 #ifdef DRAWING1
1668 inquire();
1669 penthickness(5);
1670 DrawMark(vi.r);
1671 penthickness(1);
1672 inquire();
1673 #endif
1674
1675 MeshError(1001,this);
1676 }
1677
1678
1679 quadtree->Add(vi);
1680 #ifdef DRAWING1
1681 DrawMark(vi.r);
1682 #endif
1683 assert (tcvi && tcvi->det >= 0) ;// internal
1684 Add(vi,tcvi,dete);
1685 NbSwap += vi.Optim(1);
1686 iv++;
1687 // }
1688 }
1689 if (verbosity>3)
1690 {
1691 cout << " Nb Of New Point " << iv ;
1692 cout << " Nb swap = " << NbSwap << " to split internal edges with border vertices" ;}
1693
1694 nbv = iv;
1695 }
1696 if (NbSplitEdge > nbv-nbvold)
1697 cout << " Warning not enough vertices to split all internal edges " << endl
1698 << " we lost " << NbSplitEdge - ( nbv-nbvold) << " Edges Sorry " << endl;
1699 if (verbosity>2)
1700 cout << "SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge << endl;
1701 return NbSplitEdge;
1702 }
InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)1703 Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)
1704 {
1705 Real8 seuil= 1.414/2 ;// for two close point
1706 Int4 i;
1707 // insertion part ---
1708
1709 const Int4 nbvnew = nbv-nbvold;
1710 if (verbosity>5)
1711 cout << " Try to Insert the " <<nbvnew<< " new points " << endl;
1712 Int4 NbSwap=0;
1713 Icoor2 dete[3];
1714
1715 // construction d'un ordre aleatoire
1716 if (! nbvnew)
1717 return 0;
1718 if (nbvnew) {
1719 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;
1720 Int4 k3 = rand()%nbvnew ;
1721 for (Int4 is3=0; is3<nbvnew; is3++) {
1722 register Int4 j = nbvold +(k3 = (k3 + PrimeNumber)% nbvnew);
1723 register Int4 i = nbvold+is3;
1724 ordre[i]= vertices + j;
1725 ordre[i]->ReferenceNumber=i;
1726 }
1727 // be carefull
1728 Int4 iv = nbvold;
1729 for (i=nbvold;i<nbv;i++)
1730 {// for all the new point
1731 Vertex & vi = *ordre[i];
1732 vi.i = toI2(vi.r);
1733 vi.r = toR2(vi.i);
1734 Real4 hx,hy;
1735 vi.m.Box(hx,hy);
1736 Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
1737 if (!quadtree->ToClose(vi,seuil,hi,hj))
1738 {
1739 // a good new point
1740 Vertex & vj = vertices[iv];
1741 Int4 j = vj.ReferenceNumber;
1742 assert( &vj== ordre[j]);
1743 if(i!=j)
1744 { // for valgring
1745 Exchange(vi,vj);
1746 Exchange(ordre[j],ordre[i]);
1747 }
1748 vj.ReferenceNumber=0;
1749 // cout << " Add " << Number(vj) << " " << vj
1750 // << " " << Number(vi) << " <--> " << Number(vj) <<endl;
1751 Triangle *tcvj = FindTriangleContening(vj.i,dete);
1752 if (tcvj && !tcvj->link)
1753 {
1754 cerr << i << " PB insert point " << Number(vj) << vj << Number(vi)
1755 << " tcvj = " << tcvj << " " << tcvj->link << endl;
1756 cerr << (*tcvj)[1] << (*tcvj)[2] << endl;
1757 tcvj = FindTriangleContening(vj.i,dete);
1758 cout << (*tcvj)[1] << (*tcvj)[2] << endl;
1759 #ifdef DRAWING1
1760 inquire();
1761 penthickness(5);
1762 DrawMark(vj.r);
1763 penthickness(1);
1764
1765 inquire();
1766 #endif
1767
1768 MeshError(1001,this);
1769 }
1770
1771
1772 quadtree->Add(vj);
1773 #ifdef DRAWING1
1774 DrawMark(vj.r);
1775 #endif
1776 assert (tcvj && tcvj->det >= 0) ;// internal
1777 Add(vj,tcvj,dete);
1778 NbSwap += vj.Optim(1);
1779 iv++;
1780 }
1781 }
1782 if (verbosity>3) {
1783 cout << " Nb Of New Point " << iv << " Nb Of To close Points " << nbv-iv ;
1784 cout << " Nb swap = " << NbSwap << " after " ;}
1785
1786 nbv = iv;
1787 }
1788
1789 #ifdef DRAWING1
1790 inquire();
1791 #endif
1792
1793 for (i=nbvold;i<nbv;i++)
1794 NbSwap += vertices[i].Optim(1);
1795 if (verbosity>3)
1796 cout << " NbSwap = " << NbSwap << endl;
1797
1798
1799 NbTSwap += NbSwap ;
1800 #ifdef DEBUG
1801 {
1802 Int4 NbErr=0;
1803 Int4 i;
1804 for (i=0;i<nbt;i++)
1805 if (triangles[i].link)
1806 {
1807 double dd =Det(triangles[i][1].r-triangles[i][0].r,triangles[i][2].r-triangles[i][0].r);
1808 if(dd <=0)
1809 {
1810 NbErr++;
1811 cerr << " det triangle i " << i << " = " << triangles[i].det ;
1812 cerr << " det triangle " << dd ;
1813 cerr << " Les trois sommets " ;
1814 cerr << Number(triangles[i][0]) << " " << Number(triangles[i][1]) << " "
1815 << Number(triangles[i][2]) << endl;
1816 cerr << "I2 " <<triangles[i][0].r << triangles[i][1].r << triangles[i][2].r << endl;
1817 cerr << "R2 " << triangles[i][0].i << triangles[i][1].i << triangles[i][2].i << endl;
1818 cerr << "I2-R2 =" <<toR2(triangles[i][0].i)-triangles[i][0].r
1819 << toR2(triangles[i][1].i)-triangles[i][1].r
1820 << toR2(triangles[i][2].i)-triangles[i][2].r << endl;
1821 }
1822 }
1823 if(NbErr) {
1824 #ifdef DRAWING
1825 Int4 kkk=0;
1826 // UnMarkUnSwapTriangle();
1827 // for (i=0;i<nbv;i++)
1828 // kkk += vertices[i].Optim(0);
1829 if(verbosity>3)
1830 cout << " Nb of swap louche " << kkk << endl;
1831 if(kkk) {
1832 for (i=0;i<nbt;i++)
1833 if (triangles[i].link)
1834 {
1835 double dd =Det(triangles[i][1].r-triangles[i][0].r,triangles[i][2].r-triangles[i][0].r);
1836 if(dd <=0)
1837 {
1838 NbErr++;
1839 cerr << " xxxdet triangle i " << i << " = " << triangles[i].det ;
1840 cerr << " xxxdet triangle " << dd ;
1841 cerr << " xxxLes trois sommets " ;
1842 cerr << Number(triangles[i][0]) << " " << Number(triangles[i][1]) << " "
1843 << Number(triangles[i][2]) << endl;
1844 cerr << triangles[i][0].r << triangles[i][1].r << triangles[i][2].r << endl;
1845 cerr << triangles[i][0].i << triangles[i][1].i << triangles[i][2].i << endl;
1846 }
1847 } }
1848 inquire();
1849 #endif
1850 // MeshError(11);
1851 }
1852
1853 }
1854 #endif
1855 return nbv-nbvold;
1856
1857 }
NewPoints(Triangles & Bh,int KeepBackVertex)1858 void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex)
1859 { // Triangles::NewPoints
1860 Int4 nbtold(nbt),nbvold(nbv);
1861 if (verbosity>2)
1862 cout << " -- Triangles::NewPoints ";
1863 if (verbosity>3)cout << " nbv (in) on Boundary = " << nbv <<endl;
1864 Int4 i,k;
1865 int j;
1866 Int4 *first_np_or_next_t = new Int4[nbtx];
1867 Int4 NbTSwap =0;
1868 // insert old point
1869 nbtold = nbt;
1870
1871 if (KeepBackVertex && (&Bh != this) && (nbv+Bh.nbv< nbvx))
1872 {
1873 // Bh.SetVertexFieldOn();
1874 for (i=0;i<Bh.nbv;i++)
1875 {
1876 Vertex & bv = Bh[i];
1877 if (!bv.on) {
1878 vertices[nbv].r = bv.r;
1879 vertices[nbv++].m = bv.m;}
1880 }
1881 int nbv1=nbv;
1882 Bh.ReMakeTriangleContainingTheVertex();
1883 InsertNewPoints(nbvold,NbTSwap) ;
1884 if (verbosity>2)
1885 cout << " (Nb of Points from background mesh = "
1886 << nbv-nbvold << " / " << nbv1-nbvold << ")" << endl;
1887 }
1888 else
1889 Bh.ReMakeTriangleContainingTheVertex();
1890
1891 Triangle *t;
1892 // generation of the list of next Triangle
1893 // at 1 time we test all the triangles
1894 Int4 Headt =0,next_t;
1895 for(i=0;i<nbt;i++)
1896 first_np_or_next_t[i]=-(i+1);
1897 // end list i >= nbt
1898 // the list of test triangle is
1899 // the next traingle on i is -first_np_or_next_t[i]
1900 int iter=0;
1901 // Big loop
1902 do {
1903 iter++;
1904 nbtold = nbt;
1905 nbvold = nbv;
1906 #ifdef DRAWING1
1907 inquire();
1908 #endif
1909
1910 // default size of IntersectionTriangle
1911
1912 i=Headt;
1913 next_t=-first_np_or_next_t[i];
1914 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])
1915 { // for each triangle t
1916 // we can change first_np_or_next_t[i]
1917 // cout << " Do the triangle " << i << " Next_t=" << next_t << endl;
1918 assert(i>=0 && i < nbt );
1919 first_np_or_next_t[i] = iter;
1920 for(j=0;j<3;j++)
1921 { // for each edge
1922 TriangleAdjacent tj(t,j);
1923 Vertex & vA = * tj.EdgeVertex(0);
1924 Vertex & vB = * tj.EdgeVertex(1);
1925
1926 if (!t->link) continue;// boundary
1927 if (t->det <0) continue;
1928 if (t->Locked(j)) continue;
1929
1930 TriangleAdjacent tadjj = t->Adj(j);
1931 Triangle * ta= tadjj;
1932
1933 if (ta->det <0) continue;
1934
1935 R2 A = vA;
1936 R2 B = vB;
1937
1938 k=Number(ta);
1939
1940 if(first_np_or_next_t[k]==iter) // this edge is done before
1941 continue; // next edge of the triangle
1942
1943 //const Int4 NbvOld = nbv;
1944 lIntTria.SplitEdge(Bh,A,B);
1945 lIntTria.NewPoints(vertices,nbv,nbvx);
1946 } // end loop for each edge
1947
1948 }// for triangle
1949
1950 #ifdef DRAWING1
1951 cout << " -------------------------------------------- " << endl;
1952 inquire();
1953 reffecran();
1954 Draw();
1955 penthickness(2);
1956 #endif
1957
1958 if (!InsertNewPoints(nbvold,NbTSwap))
1959 break;
1960
1961 for (i=nbtold;i<nbt;i++)
1962 first_np_or_next_t[i]=iter;
1963
1964 Headt = nbt; // empty list
1965 for (i=nbvold;i<nbv;i++)
1966 { // for all the triangle contening the vertex i
1967 Vertex * s = vertices + i;
1968 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
1969 Triangle * tbegin= (Triangle*) ta;
1970 Int4 kt;
1971 do {
1972 kt = Number((Triangle*) ta);
1973 if (first_np_or_next_t[kt]>0)
1974 first_np_or_next_t[kt]=-Headt,Headt=kt;
1975 assert( ta.EdgeVertex(0) == s);
1976 ta = Next(Adj(ta));
1977 } while ( (tbegin != (Triangle*) ta));
1978 }
1979
1980 } while (nbv!=nbvold);
1981
1982 delete [] first_np_or_next_t;
1983
1984 Int4 NbSwapf =0,NbSwp;
1985
1986 // bofbof
1987
1988
1989 NbSwp = NbSwapf;
1990 for (i=0;i<nbv;i++)
1991 NbSwapf += vertices[i].Optim(0);
1992 /*
1993 for (i=0;i<nbv;i++)
1994 NbSwapf += vertices[i].Optim(0);
1995 for (i=0;i<nbv;i++)
1996 NbSwapf += vertices[i].Optim(0);
1997 for (i=0;i<nbv;i++)
1998 NbSwapf += vertices[i].Optim(0);
1999 for (i=0;i<nbv;i++)
2000 NbSwapf += vertices[i].Optim(0);
2001 */
2002 NbTSwap += NbSwapf ;
2003 if (verbosity>3) cout << " " ;
2004 if (verbosity>2)
2005 cout << " Nb Of Vertices =" << nbv << " Nb of triangles = " << nbt-NbOutT
2006 << " NbSwap final = " << NbSwapf << " Nb Total Of Swap = " << NbTSwap << endl;
2007
2008
2009 }
2010
NewPointsOld(Triangles & Bh)2011 void Triangles::NewPointsOld(Triangles & Bh)
2012 { // Triangles::NewPointsOld
2013 Real8 seuil= 0.7 ;// for two neart point
2014 if (verbosity>1)
2015 cout << " begin : Triangles::NewPointsOld " << endl;
2016 Int4 i,k;
2017 int j;
2018 Int4 BeginNewPoint[3];
2019 Int4 EndNewPoint[3];
2020 #ifdef TRACETRIANGLE
2021 Int4 trace=0;
2022 #endif
2023 int step[3];
2024 Int4 *first_np_or_next_t = new Int4[nbtx];
2025 Int4 ColorEdge[3];
2026 Int4 color=-1;
2027 Triangle *t;
2028 // generation of the list of next Triangle
2029 // at 1 time we test all the triangles
2030 Int4 Headt =0,next_t;
2031 for(i=0;i<nbt;i++)
2032 first_np_or_next_t[i]=-(i+1);
2033 // end list i >= nbt
2034 // the list of test triangle is
2035 // the next Triangle on i is -first_np_or_next_t[i]
2036 Int4 nbtold,nbvold;
2037
2038 // Big loop
2039 do {
2040 nbtold = nbt;
2041 nbvold = nbv;
2042 #ifdef DRAWING1
2043 inquire();
2044 #endif
2045
2046 // default size of IntersectionTriangle
2047
2048 i=Headt;
2049 next_t=-first_np_or_next_t[i];
2050 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])
2051 { // for each triangle t
2052 // we can change first_np_or_next_t[i]
2053 #ifdef TRACETRIANGLE
2054 trace = TRACETRIANGLE <0 ? 1 : i == TRACETRIANGLE;
2055 #endif
2056 // cout << " Do the triangle " << i << " Next_t=" << next_t << endl;
2057 assert(i>=0 && i < nbt );
2058 first_np_or_next_t[i] = nbv; // to save the fist new point of triangle
2059 for(j=0;j<3;j++)
2060 { // for each edge
2061 TriangleAdjacent tj(t,j);
2062 // color++;// the color is 3i+j
2063 color = 3*i + j ;;
2064 ColorEdge[j]=color;
2065 BeginNewPoint[j]=nbv;
2066 EndNewPoint[j]=nbv-1;
2067 step[j]=1;// right sens
2068 Vertex & vA = * tj.EdgeVertex(0);
2069 Vertex & vB = * tj.EdgeVertex(1);
2070
2071 #ifdef TRACETRIANGLE
2072 if(trace) {
2073 cout << " i " << Number(vA) <<" j "<< Number(vB)
2074 << " " << t->Locked(j) ;
2075 }
2076 #endif
2077 if (!t->link) continue;// boundary
2078 if (t->det <0) continue;
2079 if (t->Locked(j)) continue;
2080
2081 TriangleAdjacent tadjj = t->Adj(j);
2082 Triangle * ta= tadjj;
2083
2084 if (ta->det <0) continue;
2085
2086 R2 A = vA;
2087 R2 B = vB;
2088
2089 k=Number(ta);
2090 // the 2 opposite vertices
2091 const Vertex & vC1 = *tj.OppositeVertex();
2092 const Vertex & vC2 = *tadjj.OppositeVertex();
2093
2094 #ifdef TRACETRIANGLE
2095 trace = trace || k == TRACETRIANGLE;
2096 if(trace) {
2097 cout << "Test Arete " << i << " AB = " << A << B
2098 << "i " <<Number(vA)<< "j" <<Number(vB);
2099 cout << " link" <<(int)t->link << " ta=" << Number( ta)
2100 << " det " <<ta->det ;
2101 cout << " hA = " <<vA.m.h << " hB = " <<vB.m.h ;
2102 cout << " loc " << ta->Locked(j) << endl;
2103 }
2104 #endif
2105
2106 if(first_np_or_next_t[k]>0) { // this edge is done before
2107 // find the color of the edge and begin , end of newpoint
2108 register int kk = t->NuEdgeTriangleAdj(j);
2109 assert ((*t)(VerticesOfTriangularEdge[j][0]) == (*ta)(VerticesOfTriangularEdge[kk][1]));
2110 assert ((*t)(VerticesOfTriangularEdge[j][1]) == (*ta)(VerticesOfTriangularEdge[kk][0]));
2111 register Int4 kolor =3*k + kk;
2112 ColorEdge[j]=kolor;
2113 register Int4 kkk= 1;
2114 step[j]=-1;// other sens
2115 BeginNewPoint[j]=0;
2116 EndNewPoint[j]=-1; // empty list
2117 for (Int4 iv=first_np_or_next_t[k];iv<nbv;iv++)
2118 if (vertices[iv].color > kolor)
2119 break; // the color is passed
2120 else if (vertices[iv].color == kolor) {
2121 EndNewPoint[j]=iv;
2122 if (kkk) // one time test
2123 kkk=0,BeginNewPoint[j]=iv;}
2124 continue; // next edge of the triangle
2125 } // end if( k < i)
2126
2127
2128 #ifdef DRAWING1
2129 penthickness(2); Move(A);Line(B); penthickness(1);
2130 #endif
2131 const Int4 NbvOld = nbv;
2132 lIntTria.SplitEdge(Bh,A,B);
2133 // Int4 NbvNp =
2134 lIntTria.NewPoints(vertices,nbv,nbvx);
2135 Int4 nbvNew=nbv;
2136 nbv = NbvOld;
2137 for (Int4 iv=NbvOld;iv<nbvNew;iv++) {
2138 vertices[nbv].color = color;
2139 vertices[nbv].ReferenceNumber=nbv;// circular link
2140 R2 C = vertices[iv].r;
2141 vertices[nbv].r = C;
2142 vertices[nbv].m = vertices[iv].m ;
2143 // test if the new point is not to close to the 2 opposite vertex
2144 R2 CC1 = C-vC1 , CC2 = C-vC2;
2145 if ( ( (vC1.m(CC1) + vertices[nbv].m(CC1)) > seuil)
2146 && ( (vC2.m(CC2) + vertices[nbv].m(CC2)) > seuil) )
2147 nbv++;
2148 }
2149
2150 EndNewPoint[j] = nbv-1;
2151 } // end loop for each edge
2152
2153 #ifdef TRACETRIANGLE
2154 if(trace) {
2155 // verification des point cree
2156 cout << "\n ------------ " << t->link << " " << t->det
2157 << " b " << BeginNewPoint[0] << " " << BeginNewPoint[1]
2158 << " " << BeginNewPoint[2] << " "
2159 << " e " << EndNewPoint[0] << " " << EndNewPoint[1]
2160 << " " << EndNewPoint[2] << " "
2161 << " s " << step[0] << " " << step[1] << " " << step[2] << " "
2162 << endl;
2163 }
2164 #endif
2165
2166 if (!t->link) continue;// boundary
2167 if (t->det<=0) continue;// outside
2168 // continue;
2169 for(int j0=0;j0<3;j0++)
2170 for (Int4 i0= BeginNewPoint[j0]; i0 <= EndNewPoint[j0];i0++)
2171 {
2172 // find the neart point one the opposite edge
2173 // to compute i1
2174 Vertex & vi0 = vertices[i0];
2175 int kstack = 0;
2176 Int4 stack[10];
2177 // Int4 savRef[10];
2178 int j1 = j0;
2179 while (j0 != (j1 = NextEdge[j1])) {//loop on the 2 other edge
2180 // computation of the intersection of edge j1 and DOrto
2181 // take the good sens
2182
2183 if (BeginNewPoint[j1]> EndNewPoint[j1])
2184 continue; //
2185 else if (EndNewPoint[j1] - BeginNewPoint[j1] <1) {
2186 for (Int4 ii1= BeginNewPoint[j1];ii1<=EndNewPoint[j1];ii1++)
2187 stack[kstack++] = ii1;
2188 continue;}
2189
2190
2191 int k0,k1;
2192 if (step[j1]<0) k0=1,k1=0; // reverse
2193 else k0=0,k1=1;
2194 R2 V10 = (R2)(*t)[VerticesOfTriangularEdge[j1][k0]];
2195 R2 V11 = (R2)(*t)[VerticesOfTriangularEdge[j1][k1]];
2196 R2 D = V11-V10;
2197 Real8 c0 = vi0.m(D,(R2) vi0);
2198
2199 Real8 c10 = vi0.m(D,V10);
2200 Real8 c11 = vi0.m(D,V11);
2201
2202 Real8 s;
2203 //cout << " --i0 = " << i0 << D << V10 << V11 << endl ;
2204 //cout << " c10 " << c10 << " c0 " << c0 << " c11 " << c11 << endl;
2205 if (( c10 < c0 ) && (c0 < c11))
2206 s = (c11-c0)/(c11-c10);
2207 else if (( c11 < c0 ) && (c0 < c10))
2208 s = (c11-c0) /(c11-c10);
2209 else break;
2210 R2 VP = V10*s + V11*(1-s);
2211 int sss = (c11-c10) >0 ? 1 : -1;
2212 #ifdef DRAWING1
2213 penthickness(2);
2214 Move((R2) vi0);
2215 Line(VP);
2216 penthickness(1);
2217 #endif
2218 // find the 2 point by dichotomie
2219 //cout << " t =" << Number(t) << " c0 " << c0 ;
2220 Int4 ii0 = BeginNewPoint[j1];
2221 Int4 ii1 = EndNewPoint[j1];
2222 Real8 ciii=-1,cii0=-1,cii1=-1 ;
2223 if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 )
2224 stack[kstack++] = ii0;//cout << " add+0 " << ii0;
2225 else if ( sss * ((cii1= vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 )
2226 stack[kstack++] = ii1;//cout << " add+1 " << ii1;
2227 else {
2228 while ((ii1-ii0)> 1) {
2229 Int4 iii = (ii0+ii1)/2;
2230 ciii = vi0.m( D ,(R2) vertices[iii]);
2231 //cout << " (iii = " << iii << " " << ciii << ") ";
2232 if ( sss * (ciii - c0) <0 ) ii0 = iii;
2233 else ii1 = iii;}
2234 stack[kstack++] = ii0;// cout << " add0 " << ii0;
2235 if (ii1 != ii0) stack[kstack++] = ii1;//cout << " add1 " << ii1;
2236 }
2237 #ifdef DEBUG2
2238 cout << "ii1 = " << ii1
2239 << " ii0 = " << ii0 << endl;
2240 cout << " cccc = " << cii0 << " " << ciii
2241 << " " << cii1 << " sss=" << sss << endl;
2242 #endif
2243 if (kstack >5) // bug ?
2244 cout << "NewPoints: bug????? " << kstack << " stack " << stack[kstack]<< endl;
2245 }
2246
2247 stack[kstack++] = -1; // to stop
2248 Int4 i1;
2249 kstack =0;
2250 while( (i1=stack[kstack++]) >= 0)
2251 { // the two parameter is i0 and i1
2252 assert(i1 < nbv && i1 >= 0);
2253 assert(i0 < nbv && i0 >= 0);
2254 assert(i1 != i0);
2255 R2 v01 = (R2) vertices[i1]- (R2) vertices[i0];
2256 Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01));
2257
2258
2259 #ifdef DRAWING1
2260 Move(vertices[i0].r);
2261 Line(vertices[i1].r);
2262 #endif
2263 #ifdef TRACETRIANGLE
2264 if(trace) {
2265 cout << "\n test j" << j <<" " << i0
2266 << " " << i1 << " d01=" << d01 <<endl;}
2267 #endif
2268 assert (i0 >= nbvold);
2269 assert (i1 >= nbvold);
2270 assert(i0 != i1);
2271 if (d01 == 0)
2272 break;
2273 if ( d01 < seuil)
2274 if (i1<nbvold) {
2275 // remove all the points i0;
2276 register Int4 ip,ipp;
2277 for (ip=i0;i0 != (ipp = vertices[ip].ReferenceNumber);ip=ipp)
2278 vertices[ip].ReferenceNumber = -1;// mark remove
2279 vertices[ip].ReferenceNumber = -1;// mark remove
2280 }
2281 else {
2282 // remove on of two points
2283 register Int4 ip0, ip1, ipp0,ipp1;
2284 register int kk0=1,kk1=1;
2285 // count the number of common points to compute weight w0,w1
2286 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) kk0++;
2287 for (ip1=i1;i1 != (ipp1 = vertices[ip1].ReferenceNumber);ip1=ipp1) kk1++;
2288
2289 register Real8 w0 = ((Real8) kk0)/(kk0+kk1);
2290 register Real8 w1 = ((Real8) kk1)/(kk0+kk1);
2291
2292 // make a circular link
2293 Exchange(vertices[i0].ReferenceNumber,vertices[i1].ReferenceNumber);
2294 // the new coordinate
2295 R2 C //= vertices[i0] ;
2296 = vertices[i0].r *w0 + vertices[i1].r* w1;
2297
2298 #ifdef TRACETRIANGLE
2299 if(trace) {
2300 cout << "\n ref = "<< vertices[i0].ref << " " <<vertices[i1].ref <<endl;
2301 }
2302 #endif
2303 #ifdef DRAWING1
2304 Move(vertices[i0].r);
2305 Line(vertices[i1].r);
2306 DrawMark(C);
2307 #endif
2308 // update the new point points of the list
2309 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0)
2310 vertices[ip0].r = C;
2311 vertices[ip0].r = C;
2312 }
2313 }
2314 } // for (i0= ....
2315 }// for triangle
2316
2317 #ifdef DRAWING1
2318 cout << " -------------------------------------------- " << endl;
2319 inquire();
2320 reffecran();
2321 Draw();
2322 penthickness(2);
2323 #endif
2324
2325 // remove of all the double points
2326
2327 Int4 ip,ipp,kkk=nbvold;
2328 for (i=nbvold;i<nbv;i++)
2329 if (vertices[i].ReferenceNumber>=0) {// good points
2330 // cout <<" i = " << i ;
2331 for (ip=i;i != (ipp = vertices[ip].ReferenceNumber);ip=ipp)
2332 vertices[ip].ReferenceNumber = -1;// mark remove
2333 vertices[ip].ReferenceNumber = -1;// mark remove
2334 // cout << i << " ---> " << kkk << endl;
2335 vertices[kkk] = vertices[i];
2336 vertices[kkk].i = toI2(vertices[kkk].r);
2337 #ifdef DRAWING1
2338 DrawMark(vertices[kkk]);
2339 #endif
2340 vertices[kkk++].ReferenceNumber = 0;
2341
2342 }
2343 #ifdef DRAWING1
2344 penthickness(1);
2345 #endif
2346
2347 // insertion part ---
2348
2349 const Int4 nbvnew = kkk-nbvold;
2350
2351 cout << " Remove " << nbv - kkk << " to close vertex " ;
2352 cout << " and Insert the " <<nbvnew<< " new points " << endl;
2353 nbv = kkk;
2354 Int4 NbSwap=0;
2355 Icoor2 dete[3];
2356
2357 // construction d'un ordre aleatoire
2358 if (! nbvnew)
2359 break;
2360 if (nbvnew) {
2361 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;
2362 Int4 k3 = rand()%nbvnew ;
2363 for (Int4 is3=0; is3<nbvnew; is3++)
2364 ordre[nbvold+is3]= &vertices[nbvold +(k3 = (k3 + PrimeNumber)% nbvnew)];
2365
2366 for (i=nbvold;i<nbv;i++)
2367 { Vertex * vi = ordre[i];
2368 Triangle *tcvi = FindTriangleContening(vi->i,dete);
2369 // Vertex * nv = quadtree->NearestVertex(vi->i.x,vi->i.y);
2370 // cout << " Neart Vertex of " << Number(vi)<< vi->i << " is "
2371 // << Number(nv) << nv->i << endl;
2372 // Int4 kt = Number(tcvi);
2373 //
2374
2375 quadtree->Add(*vi); //
2376 #ifdef DRAWING1
2377 DrawMark(vi->r);
2378 #endif
2379 assert (tcvi->det >= 0) ;// internal
2380 Add(*vi,tcvi,dete),NbSwap += vi->Optim(1);
2381 }
2382 }
2383 cout << " Nb swap = " << NbSwap << " after " ;
2384 #ifdef DRAWING1
2385 inquire();
2386 #endif
2387
2388 for (i=nbvold;i<nbv;i++)
2389 NbSwap += vertices[i].Optim(1);
2390 cout << NbSwap << endl;
2391
2392 for (i=nbtold;i<nbt;i++)
2393 first_np_or_next_t[i]=1;
2394
2395 Headt = nbt; // empty list
2396 for (i=nbvold;i<nbv;i++)
2397 { // for all the triangle contening the vertex i
2398 Vertex * s = vertices + i;
2399 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
2400 Triangle * tbegin= (Triangle*) ta;
2401 Int4 kt;
2402 do {
2403 kt = Number((Triangle*) ta);
2404 if (first_np_or_next_t[kt]>0)
2405 first_np_or_next_t[kt]=-Headt,Headt=kt;
2406 assert( ta.EdgeVertex(0) == s);
2407 ta = Next(Adj(ta));
2408 } while ( (tbegin != (Triangle*) ta));
2409 }
2410
2411
2412 } while (nbv!=nbvold);
2413 delete [] first_np_or_next_t;
2414 #ifdef DEBUG
2415 int nberr=0;
2416 for (int it=0;it<nbt;it++)
2417 if(triangles[it].det==0)
2418 if(nberr++<10) cerr << "Bug Triangle nul" << it << triangles[it] << endl;
2419 if(nberr) MeshError(992,this);
2420 #endif
2421 cout << " end : Triangles::NewPoints old nbv=" << nbv << endl;
2422
2423 }
2424
Insert()2425 void Triangles::Insert()
2426 {
2427 if (verbosity>2) cout << " -- Insert initial " << nbv << " vertices " << endl ;
2428 Triangles * OldCurrentTh =CurrentTh;
2429
2430 CurrentTh=this;
2431 double time0=CPUtime(),time1,time2,time3;
2432 SetIntCoor();
2433 Int4 i;
2434 for (i=0;i<nbv;i++)
2435 ordre[i]= &vertices[i] ;
2436
2437 // construction d'un ordre aleatoire
2438 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;
2439 Int4 k3 = rand()%nbv ;
2440 for (int is3=0; is3<nbv; is3++)
2441 ordre[is3]= &vertices[k3 = (k3 + PrimeNumber)% nbv];
2442
2443
2444
2445
2446 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
2447 if ( ++i >= nbv) {
2448 cerr << " All the vertices are aline " << endl;
2449 MeshError(998,this); }
2450
2451 // echange i et 2 dans ordre afin
2452 // que les 3 premiers ne soit pas aligne
2453 Exchange( ordre[2], ordre[i]);
2454
2455 // on ajoute un point a l'infini pour construire le maillage
2456 // afin d'avoir une definition simple des aretes frontieres
2457 nbt = 2;
2458
2459
2460 // on construit un maillage trivale forme
2461 // d'une arete et de 2 triangles
2462 // construit avec le 2 aretes orientes et
2463 Vertex * v0=ordre[0], *v1=ordre[1];
2464
2465 triangles[0](0) = 0; // sommet pour infini
2466 triangles[0](1) = v0;
2467 triangles[0](2) = v1;
2468
2469 triangles[1](0) = 0;// sommet pour infini
2470 triangles[1](2) = v0;
2471 triangles[1](1) = v1;
2472 const int e0 = OppositeEdge[0];
2473 const int e1 = NextEdge[e0];
2474 const int e2 = PreviousEdge[e0];
2475 triangles[0].SetAdj2(e0, &triangles[1] ,e0);
2476 triangles[0].SetAdj2(e1, &triangles[1] ,e2);
2477 triangles[0].SetAdj2(e2, &triangles[1] ,e1);
2478
2479 triangles[0].det = -1; // faux triangles
2480 triangles[1].det = -1; // faux triangles
2481
2482 triangles[0].SetTriangleContainingTheVertex();
2483 triangles[1].SetTriangleContainingTheVertex();
2484
2485 triangles[0].link=&triangles[1];
2486 triangles[1].link=&triangles[0];
2487
2488 #ifdef DEBUG
2489 triangles[0].check();
2490 triangles[1].check();
2491 #endif
2492 // nbtf = 2;
2493 if ( !quadtree ) quadtree = new QuadTree(this,0);
2494 quadtree->Add(*v0);
2495 quadtree->Add(*v1);
2496
2497 // on ajoute les sommets un � un
2498 Int4 NbSwap=0;
2499
2500 time1=CPUtime();
2501
2502 if (verbosity>3) cout << " -- Begin of insertion process " << endl;
2503
2504 for (Int4 icount=2; icount<nbv; icount++) {
2505 Vertex *vi = ordre[icount];
2506 // cout << " Insert " << Number(vi) << endl;
2507 Icoor2 dete[3];
2508 Triangle *tcvi = FindTriangleContening(vi->i,dete);
2509 quadtree->Add(*vi);
2510 Add(*vi,tcvi,dete);
2511 NbSwap += vi->Optim(1,0);
2512 #ifdef DRAWING1
2513 inquire();
2514 #endif
2515
2516 }// fin de boucle en icount
2517 time2=CPUtime();
2518 if (verbosity>3)
2519 cout << " NbSwap of insertion " << NbSwap
2520 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv
2521 << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv "
2522 << (float)NbUnSwap /(float) nbv
2523 <<endl;
2524 NbUnSwap = 0;
2525 // construction d'un ordre aleatoire
2526 // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ;
2527 #ifdef NBLOOPOPTIM
2528
2529 k3 = rand()%nbv ;
2530 for (int is4=0; is4<nbv; is4++)
2531 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv];
2532
2533 double timeloop = time2 ;
2534 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++)
2535 {
2536 double time000 = timeloop;
2537 Int4 NbSwap = 0;
2538 for (int is1=0; is1<nbv; is1++)
2539 NbSwap += ordre[is1]->Optim(0,0);
2540 timeloop = CPUtime();
2541 if (verbosity>3)
2542 cout << " Optim Loop "<<Nbloop<<" NbSwap: " << NbSwap
2543 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv
2544 << " CPU=" << timeloop - time000 << " s, "
2545 << " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv
2546 << endl;
2547 NbUnSwap = 0;
2548 if(!NbSwap) break;
2549 }
2550 ReMakeTriangleContainingTheVertex();
2551 // because we break the TriangleContainingTheVertex
2552 #endif
2553 #ifdef DEBUG
2554 int nberr=0;
2555 for (int it=0;it<nbt;it++)
2556 if(triangles[it].det==0)
2557 if(nberr++<10) cerr << "Bug Triangle nul" << it << triangles[it] << endl;
2558 if(nberr) MeshError(991,this);
2559 #endif
2560 time3=CPUtime();
2561 if (verbosity>4)
2562 cout << " init " << time1 - time0 << " initialisation, "
2563 << time2 - time1 << "s, insert point "
2564 << time3 -time2 << "s, optim " << endl
2565 << " Init Total Cpu Time = " << time3 - time0 << "s " << endl;
2566
2567 #ifdef DRAWING1
2568 inquire();
2569 #endif
2570 CurrentTh=OldCurrentTh;
2571 }
2572
2573
ForceBoundary()2574 void Triangles::ForceBoundary()
2575 {
2576 if (verbosity > 2)
2577 cout << " -- ForceBoundary nb of edge " << nbe << endl;
2578 int k=0;
2579 Int4 nbfe=0,nbswp=0,Nbswap=0;
2580 for (Int4 t = 0; t < nbt; t++)
2581 if (!triangles[t].det)
2582 k++,cerr << " det T" << t << " = " << 0 << endl;
2583 if (k!=0) {
2584 cerr << " ther is " << k << " triangles of mes = 0 " << endl;
2585 MeshError(11,this);}
2586
2587 TriangleAdjacent ta(0,0);
2588 for (Int4 i = 0; i < nbe; i++)
2589 {
2590 nbswp = ForceEdge(edges[i][0],edges[i][1],ta);
2591
2592 if ( nbswp < 0) k++;
2593 else Nbswap += nbswp;
2594 if (nbswp) nbfe++;
2595 if ( nbswp < 0 && k < 5)
2596 {
2597 cerr << " Missing Edge " << i << " v0 = " << Number(edges[i][0]) << edges[i][0].r
2598 <<" v1= " << Number(edges[i][1]) << edges[i][1].r << " " << edges[i].on->Cracked() << " " << (Triangle *) ta ;
2599 if(ta.t)
2600 {
2601 Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1);
2602 cerr << " crossing with [" << Number(aa) << ", " << Number(bb) << "]\n";
2603 }
2604 else cerr << endl;
2605
2606 }
2607 if ( nbswp >=0 && edges[i].on->Cracked())
2608 ta.SetCracked();
2609 }
2610
2611
2612 if (k!=0) {
2613 cerr << " they is " << k << " lost edges " << endl;
2614 cerr << " The boundary is crossing may be!" << endl;
2615 MeshError(10,this);
2616 }
2617 for (Int4 j=0;j<nbv;j++)
2618 Nbswap += vertices[j].Optim(1,0);
2619 if (verbosity > 3)
2620 cout << " Nb of inforced edge = " << nbfe << " Nb of Swap " << Nbswap << endl;
2621
2622 }
2623
FindSubDomain(int OutSide=0)2624 void Triangles::FindSubDomain(int OutSide=0)
2625 {
2626 //#define DRAWING1
2627
2628 if (verbosity >2)
2629 {
2630 if (OutSide)
2631 cout << " -- Find all external sub-domain ";
2632 else
2633 cout << " -- Find all internal sub-domain ";
2634 if(verbosity>99)
2635 {
2636
2637 for(int i=0;i<nbt;++i)
2638 cout << i<< " " << triangles[i] << endl;
2639 }
2640
2641 }
2642 // if (verbosity > 4) cout << " OutSide=" << OutSide << endl;
2643 short * HeapArete = new short[nbt];
2644 Triangle ** HeapTriangle = new Triangle* [nbt];
2645 Triangle *t,*t1;
2646 Int4 k,it;
2647
2648 for (Int4 itt=0;itt<nbt;itt++)
2649 triangles[itt].link=0; // par defaut pas de couleur
2650 #ifdef DRAWING1
2651 reffecran();
2652 #endif
2653
2654 Int4 NbSubDomTot =0;
2655 for ( it=0;it<nbt;it++) {
2656 if ( ! triangles[it].link ) {
2657 t = triangles + it;
2658 NbSubDomTot++;; // new composante connexe
2659 Int4 i = 0; // niveau de la pile
2660 t->link = t ; // sd forme d'un triangle cicular link
2661 #ifdef DRAWING1
2662 t->Draw(NbSubDomTot-1);
2663 #endif
2664
2665
2666 HeapTriangle[i] =t ;
2667 HeapArete[i] = 3;
2668
2669 while (i >= 0) // boucle sur la pile
2670 { while ( HeapArete[i]--) // boucle sur les 3 aretes
2671 {
2672 int na = HeapArete[i];
2673 Triangle * tc = HeapTriangle[i]; // triangle courant
2674 if( ! tc->Locked(na)) // arete non frontiere
2675 {
2676 Triangle * ta = tc->TriangleAdj(na) ; // n� triangle adjacent
2677 if (ta->link == 0 ) // non deja chainer => on enpile
2678 {
2679 i++;
2680 #ifdef DRAWING1
2681 ta->Draw(NbSubDomTot-1);
2682 #endif
2683 ta->link = t->link ; // on chaine les triangles
2684 t->link = ta ; // d'un meme sous domaine
2685 HeapArete[i] = 3; // pour les 3 triangles adjacents
2686 HeapTriangle[i] = ta;
2687 }}
2688 } // deplie fin de boucle sur les 3 adjacences
2689 i--;
2690 }
2691 }
2692 }
2693
2694 // supression de tous les sous domaine infini <=> contient le sommet NULL
2695 it =0;
2696 NbOutT = 0;
2697 while (it<nbt) {
2698 if (triangles[it].link)
2699 {
2700 if (!( triangles[it](0) && triangles[it](1) && triangles[it](2) ))
2701 {
2702 // infini triangle
2703 NbSubDomTot --;
2704 // cout << " triangle infini " << it << triangles[it] << endl;
2705 t=&triangles[it];
2706 NbOutT--; // on fait un coup de trop.
2707 while (t){ // cout << Number(t) << " " << endl;
2708 NbOutT++;
2709 t1=t;
2710 t=t->link;
2711 t1->link=0;}//while (t)
2712 }
2713 }
2714 it++;} // end while (it<nbt)
2715 if (nbt == NbOutT || !NbSubDomTot)
2716 {
2717 cout << "\n error : " << NbOutT << " " << NbSubDomTot <<" " << nbt << endl;
2718 cerr << "Error: The boundary is not close => All triangles are outside " << endl;
2719 MeshError(888,this);
2720 }
2721
2722 delete [] HeapArete;
2723 delete [] HeapTriangle;
2724
2725
2726 if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains )
2727 { // No geom sub domain
2728 Int4 i;
2729 if (subdomains) delete [] subdomains;
2730 subdomains = new SubDomain[ NbSubDomTot];
2731 NbSubDomains= NbSubDomTot;
2732 for ( i=0;i<NbSubDomains;i++) {
2733 subdomains[i].head=NULL;
2734 subdomains[i].ref=i+1;
2735 }
2736 Int4 * mark = new Int4[nbt];
2737 for (it=0;it<nbt;it++)
2738 mark[it]=triangles[it].link ? -1 : -2;
2739
2740 it =0;
2741 k = 0;
2742 while (it<nbt) {
2743 if (mark[it] == -1) {
2744 t1 = & triangles[it];
2745 t = t1->link;
2746 mark[it]=k;
2747 #ifdef DRAWING1
2748 t1->Draw(k);
2749 #endif
2750 subdomains[k].head = t1;
2751 // cout << " New -- " << Number(t1) << " " << it << endl;
2752 do {// cout << " k " << k << " " << Number(t) << endl;
2753 mark[Number(t)]=k;
2754 #ifdef DRAWING1
2755 t->Draw(k);
2756 #endif
2757 t=t->link;
2758 } while (t!=t1);
2759 #ifdef DRAWING1
2760 t1->Draw(k);
2761 #endif
2762 mark[it]=k++;}
2763 // else if(mark[it] == -2 ) triangles[it].Draw(999);
2764 it++;} // end white (it<nbt)
2765 assert(k== NbSubDomains);
2766 if(OutSide)
2767 {
2768 // to remove all the sub domain by parity adjacents
2769 // because in this case we have only the true boundary edge
2770 // so teh boundary is manifold
2771 Int4 nbk = NbSubDomains;
2772 while (nbk)
2773 for (it=0;it<nbt && nbk ;it++)
2774 for (int na=0;na<3 && nbk ;na++)
2775 {
2776 Triangle *ta = triangles[it].TriangleAdj(na);
2777 Int4 kl = ta ? mark[Number(ta)] : -2;
2778 Int4 kr = mark[it];
2779 if(kr !=kl) {
2780 //cout << kl << " " << kr << " rl " << subdomains[kl].ref
2781 // << " rr " << subdomains[kr].ref ;
2782 if (kl >=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0)
2783 nbk--,subdomains[kr].ref=subdomains[kl].ref-1;
2784 if (kr >=0 && subdomains[kr].ref <0 && kl >=0 && subdomains[kl].ref>=0)
2785 nbk--,subdomains[kl].ref=subdomains[kr].ref-1;
2786 if(kr<0 && kl >=0 && subdomains[kl].ref>=0)
2787 nbk--,subdomains[kl].ref=-1;
2788 if(kl<0 && kr >=0 && subdomains[kr].ref>=0)
2789 nbk--,subdomains[kr].ref=-1;
2790 // cout << " after \t "
2791 // << kl << subdomains[kl].ref << " rr " << kr
2792 // << subdomains[kr].ref << endl;
2793 }
2794 }
2795 // cout << subdomains[0].ref << subdomains[1].ref << endl;
2796 Int4 j=0;
2797 for ( i=0;i<NbSubDomains;i++)
2798 if((-subdomains[i].ref) %2) { // good
2799 //cout << " sudom ok = " << i << " " << subdomains[i].ref
2800 // << " " << (-subdomains[i].ref) %2 << endl;
2801 if(i != j)
2802 Exchange(subdomains[i],subdomains[j]);
2803 j++;}
2804 else
2805 { //cout << " remove sub domain " << i << endl;
2806 t= subdomains[i].head;
2807 while (t){// cout << Number(t) << " " << endl;
2808 NbOutT++;
2809 t1=t;
2810 t=t->link;
2811 t1->link=0;}//while (t)
2812 }
2813
2814 if(verbosity>4)
2815 cout << " Number of remove sub domain (OutSideMesh) =" << NbSubDomains-j << endl;
2816 NbSubDomains=j;
2817 }
2818
2819 delete [] mark;
2820
2821 }
2822 else
2823 { // find the head for all sub domaine
2824 if (Gh.NbSubDomains != NbSubDomains && subdomains)
2825 delete [] subdomains, subdomains=0;
2826 if (! subdomains )
2827 subdomains = new SubDomain[ Gh.NbSubDomains];
2828 NbSubDomains =Gh.NbSubDomains;
2829 if(verbosity>4)
2830 cout << " find the " << NbSubDomains << " sub domain " << endl;
2831 Int4 err=0;
2832 ReMakeTriangleContainingTheVertex();
2833 Int4 * mark = new Int4[nbt];
2834 Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge();
2835
2836 for (it=0;it<nbt;it++)
2837 mark[it]=triangles[it].link ? -1 : -2;
2838 Int4 inew =0;
2839 for (Int4 i=0;i<NbSubDomains;i++)
2840 {
2841 GeometricalEdge &eg = *Gh.subdomains[i].edge;
2842 subdomains[i].ref = Gh.subdomains[i].ref;
2843 int ssdlab = subdomains[i].ref;
2844 // by carefull is not easy to find a edge create from a GeometricalEdge
2845 // see routine MakeGeometricalEdgeToEdge
2846 Edge &e = *GeometricalEdgetoEdge[Gh.Number(eg)];
2847 assert(&e);
2848 Vertex * v0 = e(0),*v1 = e(1);
2849 Triangle *t = v0->t;
2850 int sens = Gh.subdomains[i].sens;
2851 // test if ge and e is in the same sens
2852 // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref
2853 // << " ref edge =" << eg.ref << " sens " << sens ;
2854 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0)
2855 sens = -sens ;
2856 subdomains[i].sens = sens;
2857 subdomains[i].edge = &e;
2858 // cout << " sens " << sens << " in geom " << eg[0].r << eg[1].r << " in mesh " << e[0].r << e[1].r << endl;
2859 // cout << " v0 , v1 = " << Number(v0) << " " << Number(v1) << endl;
2860 assert(t && sens);
2861
2862 TriangleAdjacent ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges
2863
2864 while (1)
2865 {
2866 assert( v0 == ta.EdgeVertex(1) );
2867 // cout << " recherche " << Number( ta.EdgeVertex(0)) << endl;
2868 if (ta.EdgeVertex(0) == v1) { // ok we find the edge
2869 if (sens>0)
2870 subdomains[i].head=t=Adj(ta);
2871 else
2872 subdomains[i].head=t=ta;
2873 //cout << " triangle =" << Number(t) << " = " << (*t)[0].r << (*t)[1].r << (*t)[2].r << endl;
2874 if(t<triangles || t >= triangles+nbt || t->det < 0
2875 || t->link == 0) // Ajoute aout 200
2876 {
2877 cerr << " Error in the def of sub domain "<<i
2878 << " form border " << NbSubDomains - i << "/" << NbSubDomains
2879 << ": Bad sens " << Gh.Number(eg) <<" "<< sens << endl;
2880 err++;
2881 break;}
2882 Int4 it = Number(t);
2883 if (mark[it] >=0) {
2884 if(verbosity>10)
2885 cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref
2886 << " is previouly defined with " <<mark[it] << " ref = " << subdomains[mark[it]].ref
2887 << " skip this def " << endl;
2888 break;}
2889 if(i != inew)
2890 Exchange(subdomains[i],subdomains[inew]);
2891 inew++;
2892 Triangle *tt=t;
2893 Int4 kkk=0;
2894 do
2895 {
2896 kkk++;
2897 assert(mark[Number(tt)]<0);
2898 #ifdef DRAWING1
2899 tt->Draw(i);
2900 #endif
2901 mark[Number(tt)]=i;
2902 tt=tt->link;
2903 } while (tt!=t);
2904 if(verbosity>7)
2905 cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl;
2906 break;}
2907 ta = Previous(Adj(ta));
2908 if(t == (Triangle *) ta) {
2909 err++;
2910 cerr << " Error in the def of sub domain " << i
2911 << " edge=" << Gh.Number(eg) << " " << sens << endl;
2912 break;}
2913 // cout << " NB of remove subdomain " << NbSubDomTot-NbSubDomains<< endl;
2914
2915 }
2916
2917 }
2918 if (err) MeshError(777,this);
2919
2920 if (inew < NbSubDomains) {
2921 if (verbosity>5)
2922 cout << " Warning: We remove " << NbSubDomains-inew << " SubDomains " << endl;
2923 NbSubDomains=inew;}
2924
2925
2926 for (it=0;it<nbt;it++)
2927 if ( mark[it] ==-1 )
2928 NbOutT++,triangles[it].link =0;
2929 delete [] GeometricalEdgetoEdge;
2930 delete [] mark;
2931
2932 }
2933
2934 #ifdef DRAWING1
2935 inquire();
2936 #endif
2937 NbOutT=0;
2938 for (it=0;it<nbt;it++)
2939 if(!triangles[it].link) NbOutT++;
2940 if (verbosity> 4)
2941 cout << " " ;
2942 if (verbosity> 2)
2943 cout << " Nb of Sub borned Domain = " << NbSubDomTot << " NbOutTriangles = " << NbOutT <<endl;
2944 #ifdef DRAWING1
2945 inquire();
2946 #endif
2947
2948 //#undef DRAWING1
2949
2950
2951 }
ReNumberingVertex(Int4 * renu)2952 void Triangles::ReNumberingVertex(Int4 * renu)
2953 {
2954 // warning be carfull because pointeur
2955 // from on mesh to over mesh
2956 // -- so do ReNumbering a the beginning
2957 Vertex * ve = vertices+nbv;
2958 Int4 it,ie,i;
2959
2960 for ( it=0;it<nbt;it++)
2961 triangles[it].ReNumbering(vertices,ve,renu);
2962
2963 for ( ie=0;ie<nbe;ie++)
2964 edges[ie].ReNumbering(vertices,ve,renu);
2965
2966 for (i=0;i< NbVerticesOnGeomVertex;i++)
2967 {
2968 Vertex *v = VerticesOnGeomVertex[i].mv;
2969 if (v>=vertices && v < ve)
2970 VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];
2971 }
2972
2973 for (i=0;i< NbVerticesOnGeomEdge;i++)
2974 {
2975 Vertex *v =VerticesOnGeomEdge[i].mv;
2976 if (v>=vertices && v < ve)
2977 VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)];
2978 }
2979
2980 for (i=0;i< NbVertexOnBThVertex;i++)
2981 {
2982 Vertex *v=VertexOnBThVertex[i].v;
2983 if (v>=vertices && v < ve)
2984 VertexOnBThVertex[i].v=vertices+renu[Number(v)];
2985 }
2986
2987 for (i=0;i< NbVertexOnBThEdge;i++)
2988 {
2989 Vertex *v=VertexOnBThEdge[i].v;
2990 if (v>=vertices && v < ve)
2991 VertexOnBThEdge[i].v=vertices+renu[Number(v)];
2992 }
2993
2994 // move the Vertices without a copy of the array
2995 // be carefull not trivial code
2996 Int4 j;
2997 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
2998 if (renu[it] >= 0) // a new sub cycle
2999 {
3000 i=it;
3001 Vertex ti=vertices[i],tj;
3002 while ( (j=renu[i]) >= 0)
3003 { // i is old, and j is new
3004 renu[i] = -1-renu[i]; // mark
3005 tj = vertices[j]; // save new
3006 vertices[j]= ti; // new <- old
3007 i=j; // next
3008 ti = tj;
3009 }
3010 }
3011 if (quadtree)
3012 { delete quadtree;
3013 quadtree = new QuadTree(this);
3014 }
3015 for ( it=0;it<nbv;it++)
3016 renu[i]= -renu[i]-1;
3017
3018 }
ReNumberingTheTriangleBySubDomain(bool justcompress)3019 void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress)
3020 {
3021 Int4 *renu= new Int4[nbt];
3022 register Triangle *t0,*t,*te=triangles+nbt;
3023 register Int4 k=0,it,i,j;
3024
3025 for ( it=0;it<nbt;it++)
3026 renu[it]=-1; // outside triangle
3027 for ( i=0;i<NbSubDomains;i++)
3028 {
3029 t=t0=subdomains[i].head;
3030 assert(t0); // no empty sub domain
3031 do {
3032 Int4 kt = Number(t);
3033 assert(kt>=0 && kt < nbt );
3034 assert(renu[kt]==-1);
3035 renu[kt]=k++;
3036 }
3037 while (t0 != (t=t->link));
3038 }
3039 if (verbosity>9)
3040 cout << " number of inside triangles " << k << " nbt = " << nbt << endl;
3041 // take is same numbering if possible
3042 if(justcompress)
3043 for ( k=0,it=0;it<nbt;it++)
3044 if(renu[it] >=0 )
3045 renu[it]=k++;
3046
3047 // put the outside triangles at the end
3048 for ( it=0;it<nbt;it++)
3049 if (renu[it]==-1)
3050 renu[it]=k++;
3051
3052 assert(k == nbt);
3053 // do the change on all the pointeur
3054 for ( it=0;it<nbt;it++)
3055 triangles[it].ReNumbering(triangles,te,renu);
3056
3057 for ( i=0;i<NbSubDomains;i++)
3058 subdomains[i].head=triangles+renu[Number(subdomains[i].head)];
3059
3060 // move the Triangles without a copy of the array
3061 // be carefull not trivial code
3062 for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
3063 if (renu[it] >= 0) // a new sub cycle
3064 {
3065 i=it;
3066 Triangle ti=triangles[i],tj;
3067 while ( (j=renu[i]) >= 0)
3068 { // i is old, and j is new
3069 renu[i] = -1; // mark
3070 tj = triangles[j]; // save new
3071 triangles[j]= ti; // new <- old
3072 i=j; // next
3073 ti = tj;
3074 }
3075 }
3076 delete [] renu;
3077 nt = nbt - NbOutT;
3078
3079 #ifdef DEBUG
3080 // verif
3081 for ( it=0;it<nbt;it++)
3082 triangles[it].check();
3083 #endif
3084 }
ConsRefTriangle(Int4 * reft) const3085 Int4 Triangles::ConsRefTriangle(Int4 *reft) const
3086 {
3087 assert(reft);
3088 register Triangle *t0,*t;
3089 register Int4 k=0, num;
3090 for (Int4 it=0;it<nbt;it++)
3091 reft[it]=-1; // outside triangle
3092 for (Int4 i=0;i<NbSubDomains;i++)
3093 {
3094 t=t0=subdomains[i].head;
3095 assert(t0); // no empty sub domain
3096 // register Int4 color=i+1;// because the color 0 is outside triangle
3097 do { k++;
3098 num = Number(t);
3099 assert(num>=0 &&num < nbt);
3100 reft[num]=i;
3101 // cout << Number(t0) << " " <<Number(t)<< " " << i << endl;
3102 }
3103 while (t0 != (t=t->link));
3104 }
3105 // NbOutT = nbt - k;
3106 if (verbosity>5)
3107 cout << " Nb of Sub Domain =" << NbSubDomains << " Nb of In Triangles " << k
3108 << " Nbt = " << nbt << " Out Triangles = " << nbt - k << endl;
3109
3110 return k;
3111
3112 }
3113
3114 /*
3115 void Triangles::ConsLinkTriangle()
3116 {
3117 for (Int4 i=0;i<NbSubDomains;i++)
3118 subdomains[i].head=0;
3119 register Triangle * t=triangles,*tend = triangles+nbt,*hst;
3120 for (;t<tend;t++)
3121 {
3122 if (t->link)
3123 {
3124 register Int4 color = t->color-1;
3125 assert(color<NbSubDomains && color>=0);
3126 if (hst=subdomains[color].head) {
3127 t->link=hst->link;
3128 hst->link=t; }
3129 else {
3130 subdomains[color].head = t;
3131 t->link=t;}// circular link
3132 }
3133 }
3134 {
3135 for (Int4 i=0;i<NbSubDomains;i++)
3136 assert(subdomains[i].head);
3137 }
3138
3139 }
3140 */
3141 /* void Triangles::RandomInit()
3142 {
3143 // Meshbegin = vertices;
3144 // Meshend = vertices + nbvx;
3145 nbv = nbvx;
3146 for (int i = 0; i < nbv; i++)
3147 {
3148 vertices[i].r.x= rand();
3149 vertices[i].r.y= rand();
3150 vertices[i].ref = 0;
3151 }
3152 }
3153 void Triangles::CubeInit(int n,int m)
3154 {
3155 // nbvx = n*m;
3156 // Meshbegin = vertices;
3157 // Meshend = vertices + nbvx;
3158 nbv = n*m;
3159 assert(nbv <= nbvx);
3160 int k =0;
3161 for (int i = 0; i < n; i++)
3162 for (int j = 0; j < m; j++)
3163 {
3164 vertices[k].r.x= i;
3165 vertices[k].r.y= j;
3166 vertices[k++].ref = 0;
3167 }
3168 }
3169 */
NearestVertex(Icoor1 i,Icoor1 j)3170 Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j)
3171 { return quadtree->NearestVertex(i,j); }
3172
3173
PreInit(Int4 inbvx,char * fname)3174 void Triangles::PreInit(Int4 inbvx,char *fname)
3175 {
3176 srand(19999999);
3177 OnDisk =0;
3178 NbRef=0;
3179 // allocGeometry=0;
3180 identity=0;
3181 NbOfTriangleSearchFind =0;
3182 NbOfSwapTriangle =0;
3183 nbiv=0;
3184 nbv=0;
3185 nbvx=inbvx;
3186 nbt=0;
3187 NbOfQuad = 0;
3188 nbtx=2*inbvx-2;
3189 NbSubDomains=0;
3190 NbVertexOnBThVertex=0;
3191 NbVertexOnBThEdge=0;
3192 VertexOnBThVertex=0;
3193 VertexOnBThEdge=0;
3194
3195 NbCrackedVertices=0;
3196 NbCrackedEdges =0;
3197 CrackedEdges =0;
3198 nbe = 0;
3199 name = fname ;
3200
3201 if (inbvx) {
3202 vertices=new Vertex[nbvx];
3203 assert(vertices);
3204 ordre=new Vertex* [this->nbvx];
3205 assert(ordre);
3206 triangles=new Triangle[nbtx];
3207 assert(triangles);}
3208 else {
3209 vertices=0;
3210 ordre=0;
3211 triangles=0;
3212 nbtx=0;
3213 }
3214 if ( name || inbvx)
3215 {
3216 time_t timer =time(0);
3217 char buf[70];
3218 strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer));
3219 counter++;
3220 char countbuf[30];
3221 sprintf(countbuf,"%d",counter);
3222 int lg =0 ;
3223 if (&BTh != this && BTh.name)
3224 lg = strlen(BTh.name)+4;
3225 identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2 + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)];
3226 identity[0]=0;
3227 if (lg)
3228 strcat(strcat(strcat(identity,"B="),BTh.name),", ");
3229
3230 if (Gh.name)
3231 strcat(strcat(identity,"G="),Gh.name);
3232 strcat(strcat(identity,";"),countbuf);
3233 strcat(identity,buf);
3234 // cout << "New MAILLAGE "<< identity << endl;
3235 }
3236
3237 quadtree=0;
3238 // edgescomponante=0;
3239 edges=0;
3240 VerticesOnGeomVertex=0;
3241 VerticesOnGeomEdge=0;
3242 NbVerticesOnGeomVertex=0;
3243 NbVerticesOnGeomEdge=0;
3244 // nbMaxIntersectionTriangles=0;
3245 // lIntTria;
3246 subdomains=0;
3247 NbSubDomains=0;
3248 // Meshbegin = vertices;
3249 // Meshend = vertices + nbvx;
3250 if (verbosity>98)
3251 cout << "Triangles::PreInit() " << nbvx << " " << nbtx
3252 << " " << vertices
3253 << " " << ordre << " " << triangles <<endl;
3254
3255 }
3256
GeomToTriangles1(Int4 inbvx,int KeepBackVertices)3257 void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
3258 {
3259 //#define DRAWING1
3260 Gh.NbRef++;// add a ref to Gh
3261
3262
3263 /*************************************************************************
3264 // methode in 2 step
3265 // 1 - compute the number of new edge to allocate
3266 // 2 - construct the edge
3267 remark:
3268 in this part we suppose to have a background mesh with the same
3269 geometry
3270
3271 To construct the discretisation of the new mesh we have to
3272 rediscretize the boundary of background Mesh
3273 because we have only the pointeur from the background mesh to the geometry.
3274 We need the abcisse of the background mesh vertices on geometry
3275 so a vertex is
3276 0 on GeometricalVertex ;
3277 1 on GeometricalEdge + abcisse
3278 2 internal
3279
3280 *************************************************************************/
3281 assert(&BTh.Gh == &Gh);
3282
3283 // if(verbosity==100) Gh.Write("/tmp/gg.gmsh");
3284 BTh.NbRef++; // add a ref to BackGround Mesh
3285 PreInit(inbvx);
3286 BTh.SetVertexFieldOn();
3287 #ifdef DRAWING
3288 if (withrgraphique)
3289 { BTh.InitDraw();
3290 reffecran();
3291 CurrentTh = this;}
3292 #endif
3293 int * bcurve = new int[Gh.NbOfCurves]; //
3294
3295 // we have 2 ways to make the loop
3296 // 1) on the geometry
3297 // 2) on the background mesh
3298 // if you do the loop on geometry, we don't have the pointeur on background,
3299 // and if you do the loop in background we have the pointeur on geometry
3300 // so do the walk on background
3301 // Int4 NbVerticesOnGeomVertex;
3302 // VertexOnGeom * VerticesOnGeomVertex;
3303 // Int4 NbVerticesOnGeomEdge;
3304 // VertexOnGeom * VerticesOnGeomEdge;
3305
3306
3307 NbVerticesOnGeomVertex=0;
3308 NbVerticesOnGeomEdge=0;
3309 //1 copy of the Required vertex
3310 int i;
3311 for ( i=0;i<Gh.nbv;i++)
3312 if (Gh[i].Required()) NbVerticesOnGeomVertex++;
3313
3314 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
3315 VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];
3316 //
3317 if( NbVerticesOnGeomVertex >= nbvx)
3318 {
3319 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl;
3320 MeshError(1,this);
3321 }
3322 assert(vertices);
3323 for (i=0;i<Gh.nbv;i++)
3324 if (Gh[i].Required()) {//Gh vertices Required
3325 vertices[nbv] = Gh[i];
3326 vertices[nbv].i = I2(0,0);
3327 Gh[i].to = vertices + nbv;// save Geom -> Th
3328 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
3329 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;
3330 nbv++;}
3331 else Gh[i].to=0;
3332 //
3333 for (i=0;i<BTh.NbVerticesOnGeomVertex;i++)
3334 {
3335 VertexOnGeom & vog = BTh.VerticesOnGeomVertex[i];
3336 if (vog.IsRequiredVertex()) {
3337 GeometricalVertex * gv = vog;
3338 Vertex *bv = vog;
3339 assert(gv->to);// use of Geom -> Th
3340 VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv);
3341 gv->to->m = bv->m; // for taking the metrix of the background mesh
3342 ;}
3343 }
3344 assert(NbVertexOnBThVertex == NbVerticesOnGeomVertex);
3345 // new stuff FH with curve
3346 // find the begin of the curve in BTh
3347 {
3348 Gh.UnMarkEdges();
3349 int bfind=0;
3350 /*
3351 cout << " nb curves = " << Gh.NbOfCurves << endl;
3352 for(int i=0;i<Gh.NbOfCurves ;i++)
3353 {
3354 cout << " Curve " << i << " begin e=" << Gh.Number(Gh.curves[i].be) << " k=" << Gh.curves[i].kb
3355 << " end e= " << Gh.Number(Gh.curves[i].ee) << " k=" << Gh.curves[i].ke << endl;
3356 }*/
3357 for (int i=0;i<Gh.NbOfCurves;i++)
3358 {
3359 bcurve[i]=-1;
3360 }
3361
3362 for (int iedge=0;iedge<BTh.nbe;iedge++)
3363 {
3364 Edge & ei = BTh.edges[iedge];
3365 for(int je=0;je<2;je++) // for the 2 extremites
3366 if (!ei.on->Mark() && ei[je].on->IsRequiredVertex() )
3367 {
3368 // a begin of curve
3369 int nc = ei.on->CurveNumber;
3370
3371 //cout << "curve " << nc << " v " << Gh.Number((GeometricalVertex *) *ei[je].on) << " "
3372 // << " e " << " " << Gh.Number(ei.on) << " vc " << Gh.Number((*Gh.curves[nc].be)[Gh.curves[nc].kb]) << endl;
3373 if(
3374 ei.on==Gh.curves[nc].be &&
3375 (GeometricalVertex *) *ei[je].on == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] // same extremity
3376 )
3377 {
3378 // cout << " find " << endl;
3379 bcurve[nc]=iedge*2+je;
3380 bfind++;
3381 }
3382 }
3383 }
3384 assert( bfind==Gh.NbOfCurves);
3385 }
3386 // method in 2 + 1 step
3387 // 0.0) compute the length and the number of vertex to do allocation
3388 // 1.0) recompute the length
3389 // 1.1) compute the vertex
3390 Int4 nbex=0,NbVerticesOnGeomEdgex=0;
3391 for (int step=0; step <2;step++)
3392 {
3393 Int4 NbOfNewPoints=0;
3394 Int4 NbOfNewEdge=0;
3395 Int4 iedge;
3396 Gh.UnMarkEdges();
3397 /* add Curve loop FH
3398 // find a starting back groud edges to walk
3399 for (iedge=0;iedge<BTh.nbe;iedge++) {
3400 Edge & ei = BTh.edges[iedge];
3401 for(int jedge=0;jedge<2;jedge++) // for the 2 extremites
3402 if (!ei.on->Mark() && ei[jedge].on->IsRequiredVertex() )
3403 {
3404 */
3405 // new code FH 2004
3406 Real8 L=0;
3407 for (int icurve=0;icurve<Gh.NbOfCurves;icurve++)
3408 {
3409 iedge=bcurve[icurve]/2;
3410 int jedge=bcurve[icurve]%2;
3411 if( ! Gh.curves[icurve].master) continue; // we skip all equi curve
3412 Edge & ei = BTh.edges[iedge];
3413 // warning: ei.on->Mark() can be change in
3414 // loop for(jedge=0;jedge<2;jedge++)
3415 // new curve
3416 // good the find a starting edge
3417 Real8 Lstep=0,Lcurve=0;// step between two points (phase==1)
3418 Int4 NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
3419
3420
3421 // cout.precision(16);
3422 for(int phase=0;phase<=step;phase++)
3423 {
3424
3425 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)
3426 {
3427
3428 int icurveequi= Gh.Number(curve);
3429
3430 if( phase == 0 && icurveequi != icurve) continue;
3431
3432 int k0=jedge,k1;
3433 Edge * pe= BTh.edges+iedge;
3434 //GeometricalEdge *ong = ei.on;
3435 int iedgeequi=bcurve[icurveequi]/2;
3436 int jedgeequi=bcurve[icurveequi]%2;
3437
3438 int k0equi=jedgeequi,k1equi;
3439 Edge * peequi= BTh.edges+iedgeequi;
3440 GeometricalEdge *ongequi = peequi->on;
3441
3442 Real8 sNew=Lstep;// abcisse of the new points (phase==1)
3443 L=0;// length of the curve
3444 Int4 i=0;// index of new points on the curve
3445 register GeometricalVertex * GA0 = *(*peequi)[k0equi].on;
3446 Vertex *A0;
3447 A0 = GA0->to; // the vertex in new mesh
3448 Vertex *A1;
3449 VertexOnGeom *GA1;
3450 Edge * PreviousNewEdge = 0;
3451 // cout << " --------------New Curve phase " << phase
3452 // << "---------- A0=" << *A0 << ei[k0] <<endl;
3453 assert (A0-vertices>=0 && A0-vertices <nbv);
3454 if(ongequi->Required() )
3455 {
3456 GeometricalVertex *GA1 = *(*peequi)[1-k0equi].on;
3457 A1 = GA1->to; //
3458 }
3459 else
3460 for(;;)
3461 {
3462 // assert(pe && BTh.Number(pe)>=0 && BTh.Number(pe)<=BTh.nbe);
3463 Edge &ee=*pe;
3464 Edge &eeequi=*peequi;
3465 k1 = 1-k0; // next vertex of the edge
3466 k1equi= 1 - k0equi;
3467
3468 assert(pe && ee.on);
3469 ee.on->SetMark();
3470 Vertex & v0=ee[0], & v1=ee[1];
3471 R2 AB= (R2) v1 - (R2) v0;
3472 Real8 L0=L,LAB;
3473 LAB = LengthInterpole(v0.m,v1.m,AB);
3474 L+= LAB;
3475 if (phase) {// computation of the new points
3476 while ((i!=NbCreatePointOnCurve) && sNew <= L) {
3477 // cout << " L0= " << L0 << " L " << L << " sN="
3478 // << sNew << " LAB=" << LAB << " NBPC =" <<NbCreatePointOnCurve<< " i " << i << endl;
3479 assert (sNew >= L0);
3480 assert(LAB);
3481
3482
3483 assert(vertices && nbv<nbvx);
3484 assert(edges && nbe < nbex);
3485 assert(VerticesOnGeomEdge && NbVerticesOnGeomEdge < NbVerticesOnGeomEdgex);
3486 // new vertex on edge
3487 A1=vertices+nbv++;
3488 GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
3489 Edge *e = edges + nbe++;
3490 Real8 se= (sNew-L0)/LAB;
3491 assert(se>=0 && se < 1.000000001);
3492 #ifdef DEBUG
3493 se = abscisseInterpole(v0.m,v1.m,AB,se); // because code \ref(xxx)
3494 #else
3495 se = abscisseInterpole(v0.m,v1.m,AB,se,1);
3496 #endif
3497 assert(se>=0 && se <= 1);
3498 //((k1==1) != (k1==k1equi))
3499 se = k1 ? se : 1. - se;
3500 se = k1==k1equi ? se : 1. - se;
3501 VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
3502 ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
3503 A1->ReferenceNumber = eeequi.ref;
3504 A1->DirOfSearch =NoDirOfSearch;
3505 //cout << icurveequi << " " << i << " " << *A1 << endl;
3506 e->on = ongequi;
3507 e->v[0]= A0;
3508 e->v[1]= A1;
3509 if(verbosity>99)
3510 cout << i << "+ New P "<< nbv-1 << " " <<sNew<< " L0=" << L0
3511 << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= "
3512 << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<" " << A1->r <<endl;
3513
3514 #ifdef DEBUG
3515 // code \label(xxx)
3516 R2 A1A0 = A1->r - A0->r;
3517 Real8 dp = LengthInterpole(A0->m,A1->m,A1A0);
3518 if (dp > 1.4) {
3519 cerr << " PB new Points "<< nbv-1 ;
3520 cerr << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= " ;
3521 cerr << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<endl;
3522 cerr << " PB calcul new on cuver points trop loin l=" << dp
3523 << " v=" << nbv-1 << " " << nbv-2 << " Lcurve = " << Lcurve << AB <<v0.m<< v1.m <<endl;
3524 }
3525
3526 #endif
3527 e->ref = eeequi.ref;
3528 e->adj[0]=PreviousNewEdge;
3529
3530 if (PreviousNewEdge)
3531 PreviousNewEdge->adj[1] = e;
3532 PreviousNewEdge = e;
3533 #ifdef DRAWING1
3534 e->Draw();
3535 // A0->Draw();
3536 A1->m.Draw(*A1);
3537 A1->Draw(Number(A1));
3538 #endif
3539 A0=A1;
3540 sNew += Lstep;
3541 // cout << " sNew = " << sNew << " L = " << L
3542 // << " ------" <<NbCreatePointOnCurve << " " << i << endl;
3543 if (++i== NbCreatePointOnCurve) break;
3544 }
3545
3546 }
3547 assert(ee.on->CurveNumber==ei.on->CurveNumber);
3548 if(verbosity>98) cout << BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() << endl;
3549 if ( ee[k1].on->IsRequiredVertex()) {
3550 assert(eeequi[k1equi].on->IsRequiredVertex());
3551 register GeometricalVertex * GA1 = *eeequi[k1equi].on;
3552 A1=GA1->to;// the vertex in new mesh
3553 assert (A1-vertices>=0 && A1-vertices <nbv);
3554 break;}
3555 if (!ee.adj[k1])
3556 {cerr << "Error adj edge " << BTh.Number(ee) << ", nbe = " << nbe
3557 << " Gh.vertices " << Gh.vertices
3558 << " k1 = " << k1 << " on=" << *ee[k1].on << endl;
3559 cerr << ee[k1].on->gv-Gh.vertices << endl;
3560 }
3561 pe = ee.adj[k1]; // next edge
3562 k0 = pe->Intersection(ee);
3563 peequi= eeequi.adj[k1equi]; // next edge
3564 k0equi=peequi->Intersection(eeequi);
3565 }// for(;;) end of the curve
3566
3567
3568 if (phase) // construction of the last edge
3569 {
3570 Edge *e = edges + nbe++;
3571 if (verbosity>10)
3572 cout << " Fin curve A1" << *A1 << " " << icurve << " <=> " << icurveequi <<"-----" <<
3573 NbCreatePointOnCurve << " == " <<i<<endl;
3574 e->on = ongequi;
3575 e->v[0]= A0;
3576 e->v[1]= A1;
3577 e->ref = peequi->ref;
3578 e->adj[0]=PreviousNewEdge;
3579 e->adj[1]=0;
3580 if (PreviousNewEdge)
3581 PreviousNewEdge->adj[1] = e;
3582 PreviousNewEdge = e;
3583 // cout << "Last new edge " << nbe << " " << " on " << Gh.Number(pe->on)
3584 // << " of curve =" <<pe->on->CurveNumber <<endl;
3585
3586
3587 #ifdef DRAWING1
3588 e->Draw();
3589 A1->Draw();
3590 A0->Draw();
3591 // inquire();
3592 #endif
3593 assert(i==NbCreatePointOnCurve);
3594
3595 }
3596 } // end loop on equi curve
3597
3598 if (!phase) { //
3599 Int4 NbSegOnCurve = Max((Int4)(L+0.5),(Int4) 1);// nb of seg
3600 Lstep = L/NbSegOnCurve;
3601 Lcurve = L;
3602 NbCreatePointOnCurve = NbSegOnCurve-1;
3603
3604 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)
3605 {
3606 NbOfNewEdge += NbSegOnCurve;
3607 NbOfNewPoints += NbCreatePointOnCurve;
3608 }
3609 if(verbosity>5)
3610 cout << icurve << " NbSegOnCurve = " << NbSegOnCurve << " Lstep="
3611 << Lstep <<" " << NbOfNewPoints<< " NBPC= " << NbCreatePointOnCurve <<endl;
3612 // do'nt
3613 // if(NbCreatePointOnCurve<1) break;
3614 }
3615 }//for(phase;;)
3616 /* modif FH add Curve class
3617 }}//for (iedge=0;iedge<BTh.nbe;iedge++)
3618 */
3619 // new code Curve class
3620 } // end of curve loop
3621 // end new code
3622 // do the allocation
3623 if(step==0)
3624 {
3625 //if(!NbOfNewPoints) break;// nothing ????? bug
3626 if(nbv+NbOfNewPoints > nbvx)
3627 {
3628 cerr << " Too much vertices on geometry " << nbv+NbOfNewPoints << " >= " << nbvx << endl;
3629 MeshError(3,this);
3630 }
3631 //cout << " NbOfNewEdge" << NbOfNewEdge << " NbOfNewPoints " << NbOfNewPoints << endl;
3632 edges = new Edge[NbOfNewEdge];
3633 nbex = NbOfNewEdge;
3634 if(NbOfNewPoints) { //
3635 VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints];
3636 NbVertexOnBThEdge =NbOfNewPoints;
3637 VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints];
3638 NbVerticesOnGeomEdgex = NbOfNewPoints; }
3639 NbOfNewPoints =0;
3640 NbOfNewEdge = 0;
3641 }
3642 } // for(step;;)
3643 assert(nbe);
3644
3645 delete [] bcurve;
3646
3647
3648 #ifdef DRAWING1
3649 reffecran();
3650 InitDraw();
3651 Draw();
3652 inquire();
3653 #endif
3654
3655 Insert();
3656 ForceBoundary();
3657 FindSubDomain();
3658
3659 #ifdef DRAWING1
3660 reffecran();
3661 Draw();
3662 inquire();
3663 #endif
3664 // NewPointsOld(*this) ;
3665 // BTh.ReMakeTriangleContainingTheVertex(); // FH change => put in NewPoints
3666 // for (Int4 iv=0;iv<BTh.nbv;iv++)
3667 // BTh[iv].i = toI2(BTh[iv].r);
3668 NewPoints(BTh,KeepBackVertices) ;
3669 CurrentTh = 0;
3670 //#undef DRAWING1
3671 }
3672
GeomToTriangles0(Int4 inbvx)3673 void Triangles::GeomToTriangles0(Int4 inbvx)
3674 {
3675 Gh.NbRef++;// add a ref to GH
3676
3677
3678 Int4 i,NbOfCurves=0,NbNewPoints,NbEdgeCurve;
3679 Real8 lcurve, lstep,s;
3680 #ifdef DRAWING
3681 if (withrgraphique)
3682 {
3683 Gh.InitDraw() ;
3684 CurrentTh = this; }
3685 #endif
3686
3687 R2 AB;
3688 GeometricalVertex *a,*b;
3689 Vertex *va,*vb;
3690 GeometricalEdge * e;
3691 PreInit(inbvx);
3692 int background = &BTh != this;
3693 // int SameGeom = background && (&BTh.Gh == &Gh);
3694 nbv = 0;
3695 NbVerticesOnGeomVertex=0;
3696 NbVerticesOnGeomEdge=0;
3697 for (i=0;i<Gh.nbv;i++)
3698 if (Gh[i].Required() && Gh[i].IsThe() ) NbVerticesOnGeomVertex++;
3699 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
3700 //
3701 if( NbVerticesOnGeomVertex >= nbvx)
3702 {
3703 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl;
3704 MeshError(1,this);
3705 }
3706 for (i=0;i<Gh.nbv;i++)
3707 if (Gh[i].Required()&& Gh[i].IsThe() ) {//Gh vertices Required
3708 if (nbv < nbvx)
3709 vertices[nbv] = Gh[i];
3710 Gh[i].to = vertices + nbv;// save Geom -> Th
3711 VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]);
3712 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;
3713 nbv++;
3714 }
3715 // assert( Gh.nbv < nbvx);
3716
3717 // Method in 2 step: 0 and 1
3718 // 1) compute de nb of edge
3719 // 2) construct the edge
3720 // generation of the curves
3721 assert(! edges);
3722 #ifdef DRAWING1
3723 reffecran();
3724 #endif
3725 // 2 step
3726 // --step=0 to compute the number of edges + alloc at end
3727 // --step=1 to construct the edges
3728 for (int step=0;step<2;step++)
3729 {// for (int step=0;step<2;step++)
3730 Int4 nbex = 0;
3731 nbe = 0;
3732 Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
3733 // cout << " -------------- step =" << step << endl;
3734 Gh.UnMarkEdges();
3735 NbOfCurves = 0;
3736 for (i=0;i<Gh.nbe;i++) {
3737 GeometricalEdge & ei = Gh.edges[i];
3738 if (!ei.Dup()) // a good curve (not dup )
3739 for(int j=0;j<2;j++)
3740 if (!ei.Mark() && ei[j].Required()) {
3741 // warning ei.Mark() can be change in loop for(j=0;j<2;j++)
3742 // cout << " New curve = " << NbOfCurves << endl;
3743 Int4 nbvend =0;
3744
3745 Edge * PreviousNewEdge=0;
3746
3747 lstep = -1;//to do not create points
3748 if(ei.Required())
3749 {
3750 if (j==0)
3751 if(step==0)
3752 nbe++;
3753 else
3754 {
3755 e = & ei;
3756 a=ei(0)->The();
3757 b=ei(1)->The();
3758 assert(edges);
3759 edges[nbe].v[0]=a->to;
3760 edges[nbe].v[1]=b->to;;
3761 edges[nbe].ref = e->ref;
3762 edges[nbe].on = e;
3763 edges[nbe].adj[0] = 0;
3764 edges[nbe].adj[1] = 0;
3765 #ifdef DRAWING1
3766 edges[nbe].Draw();
3767 #endif
3768 nbe++;}
3769 }
3770 else
3771 { // on curve ------
3772 for ( int kstep=0;kstep<= step;kstep++)
3773 { // begin for ( int kstep=0;kstep<= step;kstep++)
3774 // if 2nd step where 2 step
3775 // -- 1 compute le length of the curve
3776 // -- create the points and edge
3777 PreviousNewEdge=0;
3778 NbNewPoints=0;
3779 NbEdgeCurve=0;
3780 assert(nbvend < nbvx);
3781 lcurve =0;
3782 s = lstep;
3783 int k=j;
3784 e = & ei;
3785 a=ei(k)->The();
3786 va = a->to;
3787 e->SetMark();
3788 // cout << " New curve " ;
3789
3790 // if SameGeo We have go in the background geometry
3791 // to find the discretisation of the curve
3792
3793 for(;;)
3794 {
3795 k = 1-k;
3796 b= (*e)(k)->The();
3797 AB = b->r - a->r;
3798 Metric MA = background ? BTh.MetricAt(a->r) :a->m ;
3799 Metric MB = background ? BTh.MetricAt(b->r) :b->m ;
3800 Real8 ledge = (MA(AB) + MB(AB))/2;
3801 //
3802 const int MaxSubEdge = 10;
3803 int NbSubEdge = 1;
3804 Real8 lSubEdge[MaxSubEdge];
3805 R2 A,B;
3806 if (ledge < 1.5)
3807 lSubEdge[0] = ledge;
3808 else {
3809 NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
3810 A= a->r;
3811 Metric MAs =MA,MBs;
3812 // cout << " lSubEdge old=" << ledge
3813 // << " new " << A << MA << endl;
3814 ledge = 0;
3815 Real8 x =0, xstep= 1. / NbSubEdge;
3816 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
3817 x += xstep;
3818 B = e->F(k ? x : 1-x);
3819 MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB);
3820 AB = A-B;
3821 lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2);
3822 // cout << " " << lSubEdge[kk] << " x " << x
3823 // << " " << A << B << MA << MB<< endl ;
3824 }
3825 // cout << endl;
3826 }
3827
3828 Real8 lcurveb = lcurve+ ledge ;
3829 while (lcurve<=s && s <= lcurveb && nbv < nbvend)
3830 {
3831 // New points
3832
3833 // Real8 aa=(lcurveb-s)/ledge;
3834 // Real8 bb=(s-lcurve)/ledge;
3835
3836 Real8 ss = s-lcurve;
3837 // 1) find the SubEdge containing ss by dichotomie
3838 int kk0=-1,kk1=NbSubEdge-1,kkk;
3839 Real8 ll0=0,ll1=ledge,llk;
3840 while (kk1-kk0>1)
3841 {
3842 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
3843 kk1=kkk,ll1=llk;
3844 else
3845 kk0=kkk,ll0=llk;}
3846 assert(kk1 != kk0);
3847
3848 Real8 sbb = (ss-ll0 )/(ll1-ll0);
3849 Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb;
3850
3851 // new vertex on edge
3852 vb = &vertices[nbv++];
3853 vb->m = Metric(aa,a->m,bb,b->m);
3854 vb->ReferenceNumber = e->ref;
3855 vb->DirOfSearch =NoDirOfSearch;
3856 Real8 abcisse = k ? bb : aa;
3857 vb->r = e->F( abcisse );
3858 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);
3859
3860 // to take in account the sens of the edge
3861
3862 s += lstep;
3863 edges[nbe].v[0]=va;
3864 edges[nbe].v[1]=vb;
3865 edges[nbe].ref = e->ref;
3866 edges[nbe].on = e;
3867 edges[nbe].adj[0] = PreviousNewEdge;
3868 if(PreviousNewEdge)
3869 PreviousNewEdge->adj[1] = &edges[nbe];
3870 #ifdef DRAWING1
3871 vb->Draw();
3872 edges[nbe].Draw();
3873 #endif
3874 PreviousNewEdge = edges + nbe;
3875 nbe++;
3876 #ifdef DEBUG1
3877 cout << " new points " << nbv-1 << " " << vb->r ;
3878 cout << " new edge " << nbe-1 << " " ;
3879 cout << va << vb << " kk0 = " << kk0
3880 << " " << kk1 << " ss=" << ss ;
3881 cout << " " << sbb << endl;
3882 cout << " " << aa << va->r << bb << vb->r
3883 <<" length=" << Norme(va->r-vb->r) << endl;
3884 cout << " s " << s << " lstep= " << lstep
3885 << " ledge= " << ledge
3886 << " lcurve= " << lcurve << endl;
3887 #endif
3888 va = vb;
3889 }
3890 lcurve = lcurveb;
3891 e->SetMark();
3892 // cout << e-Gh.edges << ", " << k << " "
3893 // <<(*e)[k].r <<" " <<(*e)[1-k].r <<" "
3894 // << lcurve<< ";; " ;
3895 a=b;
3896 if (b->Required() ) break;
3897 int kprev=k;
3898 k = e->SensAdj[kprev];// next vertices
3899 e = e->Adj[kprev];
3900 assert(e);
3901 }// for(;;)
3902 vb = b->to;
3903 // cout << endl;
3904 NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1);
3905 NbNewPoints = NbEdgeCurve-1;
3906 if(!kstep)
3907 { NbVerticesOnGeomEdge0 += NbNewPoints;
3908 NbOfCurves++;}
3909
3910 nbvend=nbv+NbNewPoints;
3911
3912 lstep = lcurve / NbEdgeCurve;
3913 // cout <<"lstep " << lstep << " lcurve "
3914 // << lcurve << " NbEdgeCurve " << NbEdgeCurve << " " <<NbVerticesOnGeomEdge0<<" " <<NbVerticesOnGeomEdge<<" step =" <<step<< endl;
3915 }
3916 // end of curve --
3917 if (edges) { // last edges of the curves
3918 edges[nbe].v[0]=va;
3919 edges[nbe].v[1]=vb;
3920 edges[nbe].ref = e->ref;
3921 edges[nbe].on = e;
3922 edges[nbe].adj[0] = PreviousNewEdge;
3923 edges[nbe].adj[1] = 0;
3924 if(PreviousNewEdge)
3925 PreviousNewEdge->adj[1] = & edges[nbe];
3926
3927
3928 #ifdef DRAWING1
3929 edges[nbe].Draw();
3930 #endif
3931 nbe++;}
3932 else
3933 nbe += NbEdgeCurve;
3934 } // end on curve ---
3935 } // if (edges[i][j].Corner())
3936 } // for (i=0;i<nbe;i++)
3937 if(!step) {
3938 // cout << "edges " << edges << " VerticesOnGeomEdge " <<VerticesOnGeomEdge << endl;
3939 assert(!edges);
3940 assert(!VerticesOnGeomEdge);
3941 edges = new Edge[nbex=nbe];
3942 if(NbVerticesOnGeomEdge0)
3943 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
3944 assert(edges);
3945 assert(VerticesOnGeomEdge || NbVerticesOnGeomEdge0 ==0);
3946 // do the vertex on a geometrical vertex
3947 NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;
3948 }
3949 else
3950 assert(NbVerticesOnGeomEdge == NbVerticesOnGeomEdge0);
3951 // cout << " Nb of Curves = " << NbOfCurves << "nbe = " << nbe
3952 // << "== " << nbex << " nbv = " << nbv << endl;
3953 assert(nbex=nbe);
3954 } // for (step=0;step<2;step++)
3955
3956 #ifdef DRAWING1
3957 reffecran();
3958 InitDraw();
3959 Draw();
3960 inquire();
3961 #endif
3962
3963 Insert();
3964 ForceBoundary();
3965 FindSubDomain();
3966
3967 #ifdef DRAWING1
3968 reffecran();
3969 Draw();
3970 inquire();
3971 #endif
3972 // NewPointsOld(*this) ;
3973 NewPoints(*this,0) ;
3974 CurrentTh = 0;
3975 }
3976
MakeGeometricalEdgeToEdge()3977 Edge** Triangles::MakeGeometricalEdgeToEdge()
3978 {
3979 assert(Gh.nbe);
3980 Edge **e= new Edge* [this->Gh.nbe];
3981
3982 Int4 i;
3983 for ( i=0;i<Gh.nbe ; i++)
3984 e[i]=NULL;
3985 for ( i=0;i<nbe ; i++)
3986 {
3987 Edge * ei = edges+i;
3988 GeometricalEdge *on = ei->on;
3989 e[Gh.Number(on)] = ei;
3990 }
3991 for ( i=0;i<nbe ; i++)
3992 for (int ii=0;ii<2;ii++) {
3993 Edge * ei = edges+i;
3994 GeometricalEdge *on = ei->on;
3995 int j= ii;
3996 while (!(*on)[j].Required()) {
3997 // cout << i << " " << ii << " j= " << j << " curve = "
3998 // << on->CurveNumber << " " << Gh.Number(on) << " on " << j
3999 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1])
4000 // << " -> " ;
4001 Adj(on,j); // next geom edge
4002 j=1-j;
4003 // cout << Gh.Number(on) << " " << j << " e[ON] = " << e[Gh.Number(on)]
4004 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) << endl;
4005 if (e[Gh.Number(on)]) break; // optimisation
4006 e[Gh.Number(on)] = ei;
4007 }
4008 }
4009
4010 int kk=0;
4011 for ( i=0;i<Gh.nbe ; i++)
4012 if (!e[i])
4013 if(kk++<10) {
4014 cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber
4015 << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1 " << Gh.Number( Gh.edges[i][1]) << endl;
4016 // assert( e[i]);
4017 }
4018 if(kk) MeshError(997,this);
4019
4020 return e;
4021 }
4022
~Triangles()4023 Triangles::~Triangles()
4024 {
4025 assert(NbRef<=0);
4026 if (CurrentTh == this) CurrentTh=0;
4027 if(verbosity>10)
4028 cout << " ~Triangles "<< this <<" "<< identity << endl;
4029 if(vertices) delete [] vertices;
4030 if(edges) delete [] edges;
4031 if(triangles) delete [] triangles;
4032 if(quadtree) delete quadtree;
4033 if(ordre) delete [] ordre;
4034 if( subdomains) delete [] subdomains;
4035 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
4036 if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex;
4037 if (name) delete [] name;
4038 if (identity) delete [] identity;
4039 if (VertexOnBThVertex) delete [] VertexOnBThVertex;
4040 if (VertexOnBThEdge) delete [] VertexOnBThEdge;
4041
4042 if (&Gh)
4043 {
4044 if (Gh.NbRef>0) Gh.NbRef--;
4045 else if (Gh.NbRef==0) delete &Gh;
4046 }
4047 if (&BTh && (&BTh != this))
4048 {
4049 if (BTh.NbRef>0) BTh.NbRef--;
4050 else if (BTh.NbRef==0) delete &BTh;
4051 }
4052 PreInit(0); // set all to zero
4053
4054 }
4055
SetIntCoor(const char * strfrom)4056 void Triangles::SetIntCoor(const char * strfrom)
4057 {
4058 pmin = vertices[0].r;
4059 pmax = vertices[0].r;
4060
4061 // recherche des extrema des vertices pmin,pmax
4062 Int4 i;
4063 for (i=0;i<nbv;i++) {
4064 pmin.x = Min(pmin.x,vertices[i].r.x);
4065 pmin.y = Min(pmin.y,vertices[i].r.y);
4066 pmax.x = Max(pmax.x,vertices[i].r.x);
4067 pmax.y = Max(pmax.y,vertices[i].r.y);
4068 }
4069 R2 DD = (pmax-pmin)*0.05;
4070 pmin = pmin-DD;
4071 pmax = pmax+DD;
4072 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
4073 assert(coefIcoor >0);
4074
4075 // generation of integer coord
4076 for (i=0;i<nbv;i++) {
4077 vertices[i].i = toI2(vertices[i].r);
4078 }
4079 #ifdef DRAWING
4080 xGrafCoef = coefIcoor;
4081 yGrafCoef = coefIcoor;
4082 xGrafOffSet = pmin.x;
4083 yGrafOffSet = pmin.y;
4084 #ifdef DRAWING1
4085 rattente(1);
4086 #endif
4087 #endif
4088
4089 // computation of the det
4090 int Nberr=0;
4091 for (i=0;i<nbt;i++)
4092 {
4093 Vertex & v0 = triangles[i][0];
4094 Vertex & v1 = triangles[i][1];
4095 Vertex & v2 = triangles[i][2];
4096 if ( &v0 && &v1 && &v2 ) // a good triangles;
4097 {
4098 triangles[i].det= det(v0,v1,v2);
4099 if (triangles[i].det <=0 && Nberr++ <10)
4100 {
4101 if(Nberr==1)
4102 if (strfrom)
4103 cerr << "+++ Fatal Error " << strfrom << "(SetInCoor) Error : area of Triangle < 0 " << endl;
4104 else
4105 cerr << "+++ Fatal Error Triangle (in SetInCoor) area of Triangle < 0" << endl;
4106 cerr << " Triangle " << i << " det (I2) = " << triangles[i].det ;
4107 cerr << " (R2) " << Det(v1.r-v0.r,v2.r-v0.r);
4108 cerr << "; The 3 vertices " << endl;
4109 cerr << Number(v0) << " " << Number(v1) << " "
4110 << Number(v2) << " : " ;
4111 cerr << v0.r << v1.r << v2.r << " ; ";
4112 cerr << v0.i << v1.i << v2.i << endl;
4113 }
4114 }
4115 else
4116 triangles[i].det= -1; // Null triangle;
4117 }
4118 if (Nberr) MeshError(899,this);
4119
4120 }
4121
FillHoleInMesh()4122 void Triangles::FillHoleInMesh()
4123 {
4124 Triangles * OldCurrentTh =CurrentTh;
4125 CurrentTh=this;
4126 // Int4 NbTold = nbt;
4127 // generation of the integer coor
4128 {
4129
4130 // double coef = coefIcoor;
4131 // recherche des extrema des vertices pmin,pmax
4132 Int4 i;
4133 if(verbosity>2)
4134 cout << " -- FillHoleInMesh: Nb of vertices =" << nbv
4135 << " Pmin = "<< pmin << " Pmax = "<< pmax << endl;
4136
4137 assert(ordre);
4138 for (i=0;i<nbv;i++)
4139 ordre[i]= 0 ;
4140
4141
4142 NbSubDomains =0;
4143
4144 // generation of the adjacence of the triangles
4145 SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv);
4146 Int4 * st = new Int4[nbt*3];
4147 for (i=0;i<nbt*3;i++)
4148 st[i]=-1;
4149 Int4 kk =0;
4150 for (i=0;i<nbe;i++)
4151 kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
4152 if (kk != nbe)
4153 {
4154 cerr << " Some Double edge in the mesh, the number is " << kk-nbe << endl;
4155 MeshError(1002,this);
4156 }
4157 for (i=0;i<nbt;i++)
4158 for (int j=0;j<3;j++)
4159 {
4160 // Int4 i0,i1;
4161 Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
4162 Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
4163 Int4 invisible = triangles[i].Hidden(j);
4164 if(st[k]==-1)
4165 st[k]=3*i+j;
4166 else if(st[k]>=0) {
4167 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)));
4168
4169 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
4170 if (invisible) triangles[i].SetHidden(j);
4171 if (k<nbe) {
4172 triangles[i].SetLocked(j);
4173 }
4174 st[k]=-2-st[k]; }
4175 else {
4176 cerr << " The edge ("
4177 << Number(triangles[i][VerticesOfTriangularEdge[j][0]])
4178 << " , "
4179 << Number(triangles[i][VerticesOfTriangularEdge[j][1]])
4180 << " ) is in more than 2 triangles " <<k <<endl;
4181 cerr << " Edge " << j << " Of Triangle " << i << endl;
4182 cerr << " Edge " << (-st[k]+2)%3 << " Of Triangle " << (-st[k]+2)/3 << endl;
4183 cerr << " Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3))
4184 << " Of Triangle " << Number(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))) << endl;
4185 MeshError(9999,this);}
4186
4187
4188 }
4189 if(verbosity>5) {
4190 cout << " On Mesh " << name << endl;
4191 cout << " - The number of Vertices = " << nbv << endl;
4192 cout << " - The number of Triangles = " << nbt << endl;
4193 cout << " - The number of given edge = " << nbe << endl;
4194 cout << " - The number of all edges = " << edge4->nb() << endl;
4195 cout << " - The Euler number = 1-Nb Of Hole = " << nbt-edge4->nb()+nbv << endl; }
4196
4197
4198 // check the consistant of edge[].adj and the geometrical required vertex
4199 Int4 k=0;
4200 for (i=0;i<edge4->nb();i++)
4201 if (st[i] >=0) // edge alone
4202 if (i < nbe)
4203 {
4204 Int4 i0=edge4->i(i);ordre[i0] = vertices+i0;
4205 Int4 i1=edge4->j(i);ordre[i1] = vertices+i1;
4206 }
4207 else {
4208 k++;
4209 if (verbosity>20 && k <20)
4210 {
4211 Int4 i0=edge4->i(i);
4212 Int4 i1=edge4->j(i);
4213 cerr << " Lose boundary edges " << i << " : " << i0 << " " << i1 << endl;
4214 }
4215 }
4216
4217 if(k != 0) {
4218 if (verbosity>20)
4219 {
4220 cout << " The given edge are " << endl;
4221 for (int i=0;i< nbe;i++)
4222 cout << " Edge " << i << " : " << Number(edges[i][0]) << " " << Number(edges[i][1])
4223 << " " << edges[i].ref << endl;
4224 }
4225 cerr << k << " boundary edges are not defined as edges " << endl;
4226 MeshError(9998,this);
4227 }
4228 // generation of the mesh with boundary points
4229 Int4 nbvb = 0;
4230 for (i=0;i<nbv;i++)
4231 {
4232 vertices[i].t=0;
4233 vertices[i].vint=0;
4234 if (ordre[i])
4235 ordre[nbvb++] = ordre[i];
4236 }
4237
4238 Triangle *savetriangles= triangles;
4239 Int4 savenbt=nbt;
4240 Int4 savenbtx=nbtx;
4241 SubDomain * savesubdomains = subdomains;
4242 subdomains = 0;
4243
4244 Int4 Nbtriafillhole = 2*nbvb;
4245 Triangle * triafillhole =new Triangle[Nbtriafillhole];
4246 if (verbosity>9)
4247 cout << " Nbtriafillhole triafillhole*" << triafillhole << endl;
4248 triangles = triafillhole;
4249
4250 nbt=2;
4251 nbtx= Nbtriafillhole;
4252
4253 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
4254 if ( ++i >= nbvb) {
4255 cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;
4256 MeshError(998,this); }
4257 Exchange( ordre[2], ordre[i]);
4258
4259 Vertex * v0=ordre[0], *v1=ordre[1];
4260
4261
4262 triangles[0](0) = 0; // sommet pour infini
4263 triangles[0](1) = v0;
4264 triangles[0](2) = v1;
4265
4266 triangles[1](0) = 0;// sommet pour infini
4267 triangles[1](2) = v0;
4268 triangles[1](1) = v1;
4269 const int e0 = OppositeEdge[0];
4270 const int e1 = NextEdge[e0];
4271 const int e2 = PreviousEdge[e0];
4272 triangles[0].SetAdj2(e0, &triangles[1] ,e0);
4273 triangles[0].SetAdj2(e1, &triangles[1] ,e2);
4274 triangles[0].SetAdj2(e2, &triangles[1] ,e1);
4275
4276 triangles[0].det = -1; // faux triangles
4277 triangles[1].det = -1; // faux triangles
4278
4279 triangles[0].SetTriangleContainingTheVertex();
4280 triangles[1].SetTriangleContainingTheVertex();
4281
4282 triangles[0].link=&triangles[1];
4283 triangles[1].link=&triangles[0];
4284
4285 #ifdef DEBUG
4286 triangles[0].check();
4287 triangles[1].check();
4288 #endif
4289 // nbtf = 2;
4290 if ( !quadtree )
4291 delete quadtree; // ->ReInitialise();
4292
4293 quadtree = new QuadTree(this,0);
4294 quadtree->Add(*v0);
4295 quadtree->Add(*v1);
4296
4297 // on ajoute les sommets un a un
4298 Int4 NbSwap=0;
4299 for (Int4 icount=2; icount<nbvb; icount++) {
4300
4301 Vertex *vi = ordre[icount];
4302 // cout << " Add vertex " << Number(vi) << endl;
4303 Icoor2 dete[3];
4304 Triangle *tcvi = FindTriangleContening(vi->i,dete);
4305 quadtree->Add(*vi);
4306 Add(*vi,tcvi,dete);
4307 NbSwap += vi->Optim(1,1);
4308
4309 #ifdef DRAWING2
4310 cout << Number(vi) << " " << NbSwap << endl;
4311 reffecran();
4312 Draw();
4313 vi->Draw();
4314 inquire();
4315 #endif
4316 }// end loop on icount
4317 #ifdef DRAWING1
4318 inquire();
4319 #endif
4320
4321 //Int4 nbtfillhole = nbt;
4322 // inforce the boundary
4323 TriangleAdjacent ta(0,0);
4324 Int4 nbloss = 0,knbe=0;
4325 for ( i = 0; i < nbe; i++)
4326 if (st[i] >=0) // edge alone => on border ... FH oct 2009
4327 {
4328 Vertex & a=edges[i][0], & b = edges[i][1];
4329 if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
4330 {
4331 knbe++;
4332 if (ForceEdge(a,b,ta)<0)
4333 nbloss++;
4334 }
4335 }
4336 if(nbloss)
4337 {
4338 cerr << " we loss some " << nbloss << " " << " edges other " << knbe << endl;
4339 MeshError(1100,this);
4340 }
4341 FindSubDomain(1);
4342 // remove all the hole
4343 // remove all the good sub domain
4344 Int4 krm =0;
4345 for (i=0;i<nbt;i++)
4346 if (triangles[i].link) // remove triangles
4347 {
4348 krm++;
4349 for (int j=0;j<3;j++)
4350 {
4351 TriangleAdjacent ta = triangles[i].Adj(j);
4352 Triangle & tta = * (Triangle *) ta;
4353 if(! tta.link) // edge between remove and not remove
4354 { // change the link of ta;
4355 int ja = ta;
4356 Vertex *v0= ta.EdgeVertex(0);
4357 Vertex *v1= ta.EdgeVertex(1);
4358 Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
4359 assert(st[k] >=0);
4360 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
4361 ta.SetLock();
4362 st[k]=-2-st[k];
4363 }
4364 }
4365 }
4366 Int4 NbTfillHoll =0;
4367 for (i=0;i<nbt;i++)
4368 if (triangles[i].link) {
4369 triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
4370 triangles[i].color=-1;
4371 }
4372 else
4373 {
4374 triangles[i].color= savenbt+ NbTfillHoll++;
4375 #ifdef DEBUG
4376 triangles[i].check();
4377 #endif
4378 }
4379 // cout << savenbt+NbTfillHoll << " " << savenbtx << endl;
4380 assert(savenbt+NbTfillHoll <= savenbtx );
4381 // copy of the outside triangles in saveTriangles
4382 for (i=0;i<nbt;i++)
4383 if(triangles[i].color>=0)
4384 {
4385 savetriangles[savenbt]=triangles[i];
4386 savetriangles[savenbt].link=0;
4387 savenbt++;
4388 }
4389 // gestion of the adj
4390 k =0;
4391 Triangle * tmax = triangles + nbt;
4392 for (i=0;i<savenbt;i++)
4393 {
4394 Triangle & ti = savetriangles[i];
4395 for (int j=0;j<3;j++)
4396 {
4397 Triangle * ta = ti.TriangleAdj(j);
4398 int aa = ti.NuEdgeTriangleAdj(j);
4399 int lck = ti.Locked(j);
4400 if (!ta) k++; // bug
4401 else if ( ta >= triangles && ta < tmax)
4402 {
4403 ta= savetriangles + ta->color;
4404 ti.SetAdj2(j,ta,aa);
4405 if(lck) ti.SetLocked(j);
4406 }
4407 }
4408 }
4409 // OutSidesTriangles = triangles;
4410 // Int4 NbOutSidesTriangles = nbt;
4411
4412 // restore triangles;
4413 nbt=savenbt;
4414 nbtx=savenbtx;
4415 delete [] triangles;
4416 delete [] subdomains;
4417 triangles = savetriangles;
4418 subdomains = savesubdomains;
4419 // cout << triangles << " <> " << OutSidesTriangles << endl;
4420 /* k=0;
4421 for (i=0;i<nbt;i++)
4422 for (int j=0;j<3;j++)
4423 if (!triangles[i].TriangleAdj(j))
4424 k++;
4425 */
4426 if (k) {
4427 cerr << "Error Nb of triangles edge alone = " << k << endl;
4428 MeshError(9997,this);
4429 }
4430 FindSubDomain();
4431 // cout << " NbTOld = " << NbTold << " == " << nbt - NbOutT << " " << nbt << endl;
4432
4433 //
4434
4435 delete edge4;
4436 delete [] st;
4437 for (i=0;i<nbv;i++)
4438 quadtree->Add(vertices[i]);
4439
4440 SetVertexFieldOn();
4441
4442 for (i=0;i<nbe;i++)
4443 if(edges[i].on)
4444 for(int j=0;j<2;j++)
4445 if (!edges[i].adj[j])
4446 if(!edges[i][j].on->IsRequiredVertex()) {
4447 cerr << " Erreur adj et sommet requis edges [" << i << "][ " << j << "]= "
4448 << Number(edges[i][j]) << " : " << " on = " << Gh.Number(edges[i].on) ;
4449 if (edges[i][j].on->OnGeomVertex())
4450 cerr << " vertex " << Gh.Number(edges[i][j].on->gv);
4451 else if (edges[i][j].on->OnGeomEdge())
4452 cerr << "Edges " << Gh.Number(edges[i][j].on->ge);
4453 else
4454 cerr << " = " << edges[i][j].on ;
4455 cerr << endl;
4456 }
4457
4458 #ifdef DRAWING1
4459 InitDraw();
4460 #endif
4461
4462 }
4463 CurrentTh=OldCurrentTh;
4464 }
4465
Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx)4466 Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) // COPY OPERATOR
4467 : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this))
4468 {
4469 Gh.NbRef++;
4470 nbvxx = Max(nbvxx,Th.nbv);
4471 Int4 i;
4472 // do all the allocation to be sure all the pointer existe
4473
4474 char * cname = 0;
4475 if (Th.name)
4476 {
4477 cname = new char[strlen(Th.name)+1];
4478 strcpy(cname,Th.name);
4479 }
4480 PreInit(nbvxx,cname);// to make the allocation
4481 // copy of triangles
4482 nt=Th.nt;
4483 nbv = Th.nbv;
4484 nbt = Th.nbt;
4485 nbiv = Th.nbiv;
4486 nbe = Th.nbe;
4487 NbSubDomains = Th.NbSubDomains;
4488 NbOutT = Th.NbOutT;
4489 NbOfQuad = Th.NbOfQuad ;
4490 NbOfSwapTriangle =0;
4491 NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex;
4492 if(NbVerticesOnGeomVertex)
4493 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
4494 NbVerticesOnGeomEdge = Th.NbVerticesOnGeomEdge;
4495 if (NbVerticesOnGeomEdge)
4496 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge] ;
4497 if (& BTh == & Th.BTh) // same back ground
4498 {
4499 BTh.NbRef++;
4500 NbVertexOnBThVertex = Th.NbVertexOnBThVertex;
4501 if(NbVertexOnBThVertex)
4502 VertexOnBThVertex = new VertexOnVertex[NbVertexOnBThVertex];
4503 NbVertexOnBThEdge = Th.NbVertexOnBThEdge;
4504 if(NbVertexOnBThEdge)
4505 VertexOnBThEdge = new VertexOnEdge[NbVertexOnBThEdge];
4506 }
4507 else
4508 { // no add on back ground mesh
4509 BTh.NbRef++;
4510 NbVertexOnBThVertex=0;
4511 VertexOnBThVertex=0;
4512 NbVertexOnBThEdge=0;
4513 VertexOnBThEdge=0;
4514 // assert (& BTh == this); // --- a voir
4515
4516 }
4517
4518
4519 if(nbe)
4520 edges = new Edge[nbe];
4521 if(NbSubDomains)
4522 subdomains = new SubDomain[NbSubDomains];
4523 pmin = Th.pmin;
4524 pmax = Th.pmax;
4525 coefIcoor = Th.coefIcoor;
4526 for(i=0;i<nbt;i++)
4527 triangles[i].Set(Th.triangles[i],Th,*this);
4528 for(i=0;i<nbe;i++)
4529 edges[i].Set(Th,i,*this);
4530 for(i=0;i<nbv;i++)
4531 vertices[i].Set(Th.vertices[i],Th,*this);
4532 for(i=0;i<NbSubDomains;i++)
4533 subdomains[i].Set(Th,i,*this);
4534 for (i=0;i<NbVerticesOnGeomVertex;i++)
4535 VerticesOnGeomVertex[i].Set(Th.VerticesOnGeomVertex[i],Th,*this);
4536 for (i=0;i<NbVerticesOnGeomEdge;i++)
4537 VerticesOnGeomEdge[i].Set(Th.VerticesOnGeomEdge[i],Th,*this);
4538 quadtree=0;
4539
4540
4541 // assert(!OutSidesTriangles);
4542 }
4543
4544 /** -- old with a bug we loss some time last swap
4545
4546 Int4 Triangle::Optim(Int2 i,int koption)
4547 {
4548 // turn in the positif sens around vertex s
4549 register Int4 NbSwap =0;
4550 register Vertex * s = ns[i];
4551 register Triangle * tbegin=0 , *t = this , *ttc;
4552 register int k=0,j = EdgesVertexTriangle[i][0],jc;
4553 tbegin=t;
4554 do {
4555 k++;
4556 #ifdef DEBUG
4557 assert( s == & (*t)[VerticesOfTriangularEdge[j][1]] );
4558 #endif
4559 #ifdef DRAWING1
4560 t->Draw();
4561 DrawMark( s->r);
4562 #endif
4563 ttc = t->at[j];
4564 jc = NextEdge[t->aa[j]&3];
4565 cout << *t << " " << VerticesOfTriangularEdge[j][1] << "\n\t try swap " << * ttc << " " << jc ;
4566 while ( ttc->swap(jc,koption)) {
4567 NbSwap++,assert(k++<20000);
4568 ttc = t->at[j];
4569 jc = NextEdge[t->aa[j]&3];
4570 cout << "\n\t s " << *ttc << " " << jc << endl;
4571 }
4572 cout << endl;
4573 t = ttc;
4574 j = NextEdge[jc];
4575 assert(k<20000);
4576 } while ( (tbegin != t));
4577
4578 return NbSwap;
4579 }
4580 */
Optim(Int2 i,int koption)4581 Int4 Triangle::Optim(Int2 i,int koption)
4582 {
4583 // turne around in positif sens
4584 Int4 NbSwap =0;
4585 #ifdef DEBUG
4586 Vertex * s = ns[i];
4587 #endif
4588 Triangle *t = this;
4589 int k=0,j =OppositeEdge[i];
4590 int jp = PreviousEdge[j];
4591 // initialise tp, jp the previous triangle & edge
4592 Triangle *tp= at[jp];
4593 jp = aa[jp]&3;
4594 #ifdef DEBUG
4595 assert(tp->at[jp] == this);
4596 #endif
4597 do {
4598 #ifdef DEBUG
4599 assert(k++<20000);
4600 assert( s == & (*t)[OppositeVertex[j]] );
4601 #endif
4602 // cout << *t << " " << j << "\n\t try swap " ;
4603 while (t->swap(j,koption))
4604 {
4605 NbSwap++;
4606 assert(k++<20000);
4607 t= tp->at[jp]; // set unchange t qnd j for previous triangles
4608 j= NextEdge[tp->aa[jp]&3];
4609 // cout << "\n\t s " << *t << " " << j << endl;
4610 #ifdef DEBUG
4611 assert( s == & (*t)[OppositeVertex[j]] );
4612 #endif
4613 }
4614 // end on this Triangle
4615 tp = t;
4616 jp = NextEdge[j];
4617
4618 t= tp->at[jp]; // set unchange t qnd j for previous triangles
4619 j= NextEdge[tp->aa[jp]&3];
4620
4621 } while( t != this);
4622 return NbSwap;
4623 }
4624
SmoothingVertex(int nbiter,Real8 omega)4625 void Triangles::SmoothingVertex(int nbiter,Real8 omega )
4626 {
4627 // if quatree exist remove it end reconstruct
4628 if (quadtree) delete quadtree;
4629 quadtree=0;
4630 ReMakeTriangleContainingTheVertex();
4631 Triangle vide; // a triangle to mark the boundary vertex
4632 Triangle ** tstart= new Triangle* [nbv];
4633 Int4 i,j,k;
4634 // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide
4635 if ( this == & BTh)
4636 for ( i=0;i<nbv;i++)
4637 tstart[i]=vertices[i].t;
4638 else
4639 for ( i=0;i<nbv;i++)
4640 tstart[i]=0;
4641 for ( j=0;j<NbVerticesOnGeomVertex;j++ )
4642 tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;
4643 for ( k=0;k<NbVerticesOnGeomEdge;k++ )
4644 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;
4645 if(verbosity>2)
4646 cout << " -- SmoothingVertex: nb Iteration = " << nbiter << " Omega = " << omega << endl;
4647 for (k=0;k<nbiter;k++)
4648 {
4649 Int4 i,NbSwap =0;
4650 Real8 delta =0;
4651 for ( i=0;i<nbv;i++)
4652 if (tstart[i] != &vide) // not a boundary vertex
4653 delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
4654 if (!NbOfQuad)
4655 for ( i=0;i<nbv;i++)
4656 if (tstart[i] != &vide) // not a boundary vertex
4657 NbSwap += vertices[i].Optim(1);
4658 if (verbosity>3)
4659 cout << " Move max = " << sqrt(delta) << " iteration = "
4660 << k << " Nb of Swap = " << NbSwap << endl;
4661 }
4662
4663 delete [] tstart;
4664 if (quadtree) quadtree= new QuadTree(this);
4665 }
MakeQuadTree()4666 void Triangles::MakeQuadTree()
4667 {
4668 if(verbosity>8)
4669 cout << " MakeQuadTree" << endl;
4670 if ( !quadtree ) quadtree = new QuadTree(this);
4671
4672
4673 #ifdef DRAWING1
4674 quadtree->Draw();
4675 rattente(1);
4676 reffecran();
4677 quadtree->Draw();
4678 rattente(1);
4679 #endif
4680
4681 }
ShowRegulaty() const4682 void Triangles::ShowRegulaty() const// Add FH avril 2007
4683 {
4684 const Real8 sqrt32=sqrt(3.)*0.5;
4685 const Real8 aireKh=sqrt32*0.5;
4686 D2 Beq(1,0),Heq(0.5,sqrt32);
4687 D2xD2 Br(D2xD2(Beq,Heq).t());
4688 D2xD2 B1r(Br.inv());
4689 /* D2xD2 BB = Br.t()*Br;
4690 cout << " BB = " << BB << " " << Br*B1r << endl;
4691 MetricAnIso MMM(BB.x.x,BB.x.y,BB.y.y);
4692 MatVVP2x2 VMM(MMM);
4693 cout << " " << VMM.lambda1 << " " << VMM.lambda2 << endl;
4694 */
4695 double gammamn=1e100,hmin=1e100;
4696 double gammamx=0,hmax=0;
4697 double beta=1e100;
4698 double beta0=0;
4699 double alpha2=0;
4700 double area=0,Marea=0;
4701 // Real8 cf= Real8(coefIcoor);
4702 // Real8 cf2= 6.*cf*cf;
4703 int nt=0;
4704 for (int it=0;it<nbt;it++)
4705 if ( triangles[it].link)
4706 {
4707 nt++;
4708 Triangle &K=triangles[it];
4709 Real8 area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
4710 area+= area3;
4711 D2xD2 B_Kt(K[0],K[1],K[2]);
4712 D2xD2 B_K(B_Kt.t());
4713 D2xD2 B1K = Br*B_K.inv();
4714 D2xD2 BK = B_K*B1r;
4715 D2xD2 B1B1 = B1K.t()*B1K;
4716 MetricAnIso MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
4717 MatVVP2x2 VMK(MK);
4718 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
4719 // cout << B_K << " * " << B1r << " == " << BK << " " << B_K*B_K.inv() << endl;
4720 Real8 betaK=0;
4721
4722 for (int j=0;j<3;j++)
4723 {
4724 Real8 he= Norme2(R2(K[j],K[(j+1)%3]));
4725 hmin=Min(hmin,he);
4726 hmax=Max(hmax,he);
4727 Vertex & v=K[j];
4728 D2xD2 M((MetricAnIso)v);
4729 betaK += sqrt(M.det());
4730
4731 D2xD2 BMB = BK.t()*M*BK;
4732 MetricAnIso M1(BMB.x.x,BMB.x.y,BMB.y.y);
4733 MatVVP2x2 VM1(M1);
4734 //cout << B_K <<" " << M << " " << he << " " << BMB << " " << VM1.lambda1 << " " << VM1.lambda2<< endl;
4735 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
4736 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);
4737 }
4738 betaK *= area3;// 1/2 (somme sqrt(det))* area(K)
4739 Marea+= betaK;
4740 // cout << betaK << " " << area3 << " " << beta << " " << beta0 << " " << area3*3*3*3 <<endl;
4741 beta=min(beta,betaK);
4742 beta0=max(beta0,betaK);
4743 }
4744 area*=3;
4745 gammamn=sqrt(gammamn);
4746 gammamx=sqrt(gammamx);
4747 cout << " -- adaptmesh Regulary: Nb triangles " << nt << " , h min " << hmin << " , h max " << hmax << endl;
4748 cout << " area = " << area << " , M area = " << Marea << " , M area/( |Khat| nt) " << Marea/(aireKh*nt) << endl;
4749 cout << " infiny-regulaty: min " << gammamn << " max " << gammamx << endl;
4750 cout << " anisomax "<< sqrt(alpha2) << ", beta max = " << 1./sqrt(beta/aireKh)
4751 << " min "<< 1./sqrt(beta0/aireKh) << endl;
4752 }
ShowHistogram() const4753 void Triangles::ShowHistogram() const
4754 {
4755
4756 const Int4 kmax=10;
4757 const Real8 llmin = 0.5,llmax=2;
4758 const Real8 lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
4759 Int4 histo[kmax+1];
4760 Int4 i,it,k, nbedges =0;
4761 for (i=0;i<=kmax;i++) histo[i]=0;
4762 for (it=0;it<nbt;it++)
4763 if ( triangles[it].link)
4764 {
4765
4766 for (int j=0;j<3;j++)
4767 {
4768 Triangle *ta = triangles[it].TriangleAdj(j);
4769 if ( !ta || !ta->link || Number(ta) >= it)
4770 {
4771 Vertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
4772 Vertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
4773 if ( !& vP || !&vQ) continue;
4774 R2 PQ = vQ.r - vP.r;
4775 Real8 l = log(LengthInterpole(vP,vQ,PQ));
4776 #ifdef DRAWING2
4777 if (l>1.4) {
4778 penthickness(3);
4779 vP.MoveTo(),vQ.LineTo();
4780 penthickness(1);
4781 cout << " l = " << l << Number(vP) << " edge = " << Number(vQ) << endl;
4782 }
4783 #endif
4784 nbedges++;
4785 k = (int) ((l - lmin)*delta);
4786 k = Min(Max(k,0L),kmax);
4787 histo[k]++;
4788 }
4789 }
4790 }
4791 cout << " -- Histogram of the unit mesh, nb of edges" << nbedges << endl <<endl;
4792
4793 cout << " length of edge in | % of edge | Nb of edges " << endl;
4794 cout << " ------------------- | ---------- | ----------- " << endl;
4795 for (i=0;i<=kmax;i++)
4796 {
4797 cout << " " ;
4798 cout.width(10);
4799 if (i==0) cout << " 0 " ;
4800 else cout << exp(lmin+i/delta) ;
4801 cout.width(); cout << "," ;
4802 cout.width(10);
4803 if (i==kmax) cout << " +infty " ;
4804 else cout << exp(lmin+(i+1)/delta) ;
4805 cout.width();cout << " | " ;
4806
4807 cout.precision(4);
4808 cout.width(6);
4809 cout << ((long) ((10000.0 * histo[i])/ nbedges))/100.0 ;
4810 cout.width();
4811 cout.precision();
4812 cout << " | " << histo[i] <<endl;
4813 }
4814 cout << " ------------------- | ---------- | ----------- " << endl <<endl;
4815
4816 }
4817
Crack()4818 int Triangles::Crack()
4819 {
4820 assert(NbCrackedEdges ==0 || NbCrackedVertices >0);
4821 for (int i=0;i<NbCrackedEdges;i++)
4822 CrackedEdges[i].Crack();
4823 return NbCrackedEdges;
4824 }
4825
UnCrack()4826 int Triangles::UnCrack()
4827 {
4828 assert(NbCrackedEdges ==0 || NbCrackedVertices >0);
4829 for (int i=0;i<NbCrackedEdges;i++)
4830 CrackedEdges[i].UnCrack();
4831 return NbCrackedEdges;
4832 }
4833
CrackMesh()4834 int Triangles::CrackMesh()
4835 {
4836 Triangles *CurrentThOld = CurrentTh;
4837 // computed the number of cracked edge
4838 int i,k;
4839 for (k=i=0;i<nbe;i++)
4840 if(edges[i].on->Cracked()) k++;
4841 if( k==0) return 0;
4842 CurrentTh = this;
4843 cout << " Nb of Cracked Edges = " << k << endl;
4844 NbCrackedEdges =k;
4845 CrackedEdges = new CrackedEdge[k];
4846 // new edge
4847 Edge * e = new Edge[ nbe + k];
4848
4849 // copy
4850 for (i=0;i<nbe;i++)
4851 e[i] = edges[i];
4852 delete edges;
4853 edges = e;
4854
4855 const int nbe0 = nbe;
4856 for (k=i=0;i<nbe0;i++) // on double les arete cracked
4857 if(edges[i].on->Cracked())
4858 {
4859 e[nbe] = e[i];
4860 // return the edge
4861 e[nbe].v[0] = e[i].v[1];
4862 e[nbe].v[1] = e[i].v[0];
4863 e[nbe].on = e[i].on->link ; // fqux
4864 CrackedEdges[k++]=CrackedEdge(edges,i,nbe);
4865 nbe++;
4866 }
4867 ReMakeTriangleContainingTheVertex() ;
4868 //
4869 int nbcrakev =0;
4870 Vertex *vlast = vertices + nbv;
4871 Vertex *vend = vertices + nbvx; // end of array
4872 for (int iv=0;iv<nbv;iv++) // vertex
4873 {
4874 Vertex & v= vertices[iv];
4875 Vertex * vv = & v;
4876 int kk=0; // nb cracked
4877 int kc=0;
4878 int kkk =0; // nb triangle with same number
4879 Triangle * tbegin = v.t;
4880 int i = v.vint;
4881 assert(tbegin && (i >= 0 ) && (i <3));
4882 // turn around the vertex v
4883 TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge
4884 int k=0;
4885 do {
4886 int kv = VerticesOfTriangularEdge[ta][1];
4887 k++;
4888 Triangle * tt (ta);
4889 if ( ta.Cracked() )
4890 {
4891 TriangleAdjacent tta=(ta.Adj());
4892 assert(tta.Cracked());
4893 if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart
4894 if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use
4895 kk++;// number of cracked edge view
4896 }
4897 if ( tt->link ) { // if good triangles store the value
4898 int it = Number(tt);
4899 assert(it < nt);
4900 (*tt)(kv)= vv; // Change the vertex of triangle
4901 if(vv<vend) {*vv= v;vv->ReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref
4902 // tt->SetTriangleContainingTheVertex();
4903 kkk++;
4904 } else if (kk) { // crack + boundary
4905 if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use
4906 }
4907
4908 ta = Next(ta).Adj();
4909 } while ( (tbegin != ta));
4910 assert(k);
4911 if (kc) nbcrakev++;
4912 }
4913
4914 if ( nbcrakev )
4915 for (int iec =0;iec < NbCrackedEdges; iec ++)
4916 CrackedEdges[iec].Set();
4917
4918 // set the ref
4919 cout << " set the ref " << endl ;
4920 NbCrackedVertices = nbcrakev;
4921 // int nbvo = nbv;
4922 nbv = vlast - vertices;
4923 int nbnewv = nbv - nbv; // nb of new vrtices
4924 if (nbcrakev && verbosity > 1 )
4925 cout << " Nb of craked vertices = " << nbcrakev << " Nb of created vertices " << nbnewv<< endl;
4926 // all the new vertices are on geometry
4927 // BOFBO-- A VOIR
4928 if (nbnewv)
4929 { //
4930 Int4 n = nbnewv+NbVerticesOnGeomVertex;
4931 Int4 i,j,k;
4932 VertexOnGeom * vog = new VertexOnGeom[n];
4933 for ( i =0; i<NbVerticesOnGeomVertex;i++)
4934 vog[i]=VerticesOnGeomVertex[i];
4935 delete [] VerticesOnGeomVertex;
4936 VerticesOnGeomVertex = vog;
4937 // loop on cracked edge
4938 Vertex * LastOld = vertices + nbv - nbnewv;
4939 for (int iec =0;iec < NbCrackedEdges; iec ++)
4940 for (k=0;k<2;k++)
4941 {
4942 Edge & e = *( k ? CrackedEdges[iec].a.edge : CrackedEdges[iec].b.edge);
4943 for (j=0;j<2;j++)
4944 {
4945 Vertex * v = e(j);
4946 if ( v >= LastOld)
4947 { // a new vertex
4948 Int4 old = v->ReferenceNumber ; // the old same vertex
4949 Int4 i = ( v - LastOld);
4950 // if the old is on vertex => warning
4951 // else the old is on edge => ok
4952 vog[i] = vog[old];
4953 // vog[i].mv = v;
4954 //g[i].ge = ;
4955 //og[i].abcisse = ;
4956 }
4957
4958 }
4959 }
4960
4961 NbVerticesOnGeomVertex = n;
4962 }
4963 SetVertexFieldOn();
4964
4965
4966 if (vlast >= vend)
4967 {
4968 cerr << " Not enougth vertices to crack the mesh we need " << nbv << " vertices " << endl;
4969 MeshError(555,this);
4970 }
4971 cout << " NbCrackedVertices " << NbCrackedVertices << endl;
4972 CurrentTh = CurrentThOld;
4973 return NbCrackedVertices;
4974 }
4975
Triangles(const Triangles & Tho,const int * flag,const int * bb)4976 Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
4977 : Gh(*(new Geometry())), BTh(*this)
4978 { // truncature
4979 //
4980
4981 char cname[] = "trunc";
4982
4983 int i,k,itadj;
4984 int kt=0;
4985 int * kk = new int [Tho.nbv];
4986 Int4 * reft = new Int4[Tho.nbt];
4987 Int4 nbInT = Tho.ConsRefTriangle(reft);
4988 Int4 * refv = new Int4[Tho.nbv];
4989
4990 for (i=0;i<Tho.nbv;i++)
4991 kk[i]=-1;
4992 for (i=0;i<Tho.nbv;i++)
4993 refv[i]=0;
4994 int nbNewBedge =0;
4995 // int nbOldBedge =0;
4996 for (i=0;i<Tho.nbt;i++)
4997 if( reft[i] >=0 && flag[i])
4998 {
4999 const Triangle & t = Tho.triangles[i];
5000 kt++;
5001 kk[Tho.Number(t[0])]=1;
5002 kk[Tho.Number(t[1])]=1;
5003 kk[Tho.Number(t[2])]=1;
5004 itadj=Tho.Number(t.TriangleAdj(0));
5005 if ( reft[itadj] >=0 && !flag[itadj])
5006 { nbNewBedge++;
5007 refv[Tho.Number(t[VerticesOfTriangularEdge[0][0]])]=bb[i];
5008 refv[Tho.Number(t[VerticesOfTriangularEdge[0][1]])]=bb[i];
5009 }
5010 itadj=Tho.Number(t.TriangleAdj(1));
5011 if ( reft[itadj] >=0 && !flag[itadj])
5012 { nbNewBedge++;
5013 refv[Tho.Number(t[VerticesOfTriangularEdge[1][0]])]=bb[i];
5014 refv[Tho.Number(t[VerticesOfTriangularEdge[1][1]])]=bb[i];}
5015 itadj=Tho.Number(t.TriangleAdj(2));
5016 if ( reft[itadj] >=0 && !flag[itadj])
5017 { nbNewBedge++;
5018 refv[Tho.Number(t[VerticesOfTriangularEdge[2][0]])]=bb[i];
5019 refv[Tho.Number(t[VerticesOfTriangularEdge[2][1]])]=bb[i];}
5020 }
5021 k=0;
5022 for (i=0;i<Tho.nbv;i++)
5023 if (kk[i]>=0)
5024 kk[i]=k++;
5025 cout << " number of vertices " << k << " remove = " << Tho.nbv - k << endl;
5026 cout << " number of triangles " << kt << " remove = " << nbInT-kt << endl;
5027 cout << " number of New boundary edge " << nbNewBedge << endl;
5028 Int4 inbvx =k;
5029 PreInit(inbvx,cname);
5030 for (i=0;i<Tho.nbv;i++)
5031 if (kk[i]>=0)
5032 {
5033 vertices[nbv] = Tho.vertices[i];
5034 if (!vertices[nbv].ref())
5035 vertices[nbv].ReferenceNumber = refv[i];
5036 nbv++;
5037 }
5038 assert(inbvx == nbv);
5039 for (i=0;i<Tho.nbt;i++)
5040 if( reft[i] >=0 && flag[i])
5041 {
5042 const Triangle & t = Tho.triangles[i];
5043 int i0 = Tho.Number(t[0]);
5044 int i1 = Tho.Number(t[1]);
5045 int i2 = Tho.Number(t[2]);
5046 assert(i0>=0 && i1 >= 0 && i2 >= 0);
5047 assert(i0<Tho.nbv && i1 <Tho.nbv && i2 <Tho.nbv);
5048 // cout <<i<< " F" << flag[i] << " T " << nbt << " = " << kk[i0] << " " << kk[i1] << " " << kk[i2] ;
5049 // cout << " OT " << i0 << " " << i1 << " " << i2 << " " << reft[i] << endl;
5050 triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]);
5051 triangles[nbt].color = Tho.subdomains[reft[i]].ref;
5052 nbt++;
5053 }
5054 assert(kt==nbt);
5055 if (nbt ==0 && nbv ==0) {
5056 cout << "Error all triangles was remove " << endl;
5057 MeshError(999,this);
5058 }
5059 delete [] kk;
5060 delete [] reft;
5061 delete [] refv;
5062 double cutoffradian = 10.0/180.0*Pi;
5063 ConsGeometry(cutoffradian);
5064 Gh.AfterRead();
5065 SetIntCoor();
5066 FillHoleInMesh();
5067
5068 assert(NbSubDomains);
5069 assert(subdomains[0].head && subdomains[0].head->link);
5070
5071 }
5072
FindTriangleContening(const I2 & B,Icoor2 dete[3],Triangle * tstart) const5073 Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const
5074 { // in: B
5075 // out: t
5076 // out : dete[3]
5077 // t the triangle and s0,s1,s2 the 3 vertices of t
5078 // in dete[3] = { det(B,s1,s2) , det(s0,B,s2), det(s0,s1,B)}
5079 // with det(a,b,c ) = -1 if one of 3 vertices a,b,c is NULL
5080 Triangle * t=0;
5081 int j,jp,jn,jj;
5082 if (tstart)
5083 t=tstart;
5084 else
5085 {
5086 assert(quadtree);
5087 Vertex *a = quadtree->NearestVertex(B.x,B.y) ;
5088
5089 if (! a || !a->t ) {
5090 if (a)
5091 {cerr << " Attention PB TriangleConteningTheVertex vertex number=" << Number(a) << endl;
5092 cerr << "We forget a call to ReMakeTriangleContainingTheVertex" << endl;}
5093 cerr << " Pb with " << B << toR2(B) << endl;
5094 MeshError(7777);
5095 }
5096 assert(a>= vertices && a < vertices+nbv);
5097 #ifdef DRAWING1
5098 a->Draw();
5099 #endif
5100 // int k=0;
5101 t = a->t;
5102 assert(t>= triangles && t < triangles+nbt);
5103
5104 }
5105 Icoor2 detop ;
5106 int kkkk =0; // number of test triangle
5107
5108 while ( t->det < 0)
5109 { // the initial triangles is outside
5110 int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
5111 assert(k0>=0); // k0 the NULL vertex
5112 int k1=NextVertex[k0],k2=PreviousVertex[k0];
5113 dete[k0]=det(B,(*t)[k1],(*t)[k2]);
5114 dete[k1]=dete[k2]=-1;
5115 if (dete[k0] > 0) // outside B
5116 return t;
5117 t = t->TriangleAdj(OppositeEdge[k0]);
5118 assert(kkkk++ < 2);
5119 }
5120
5121 jj=0;
5122 detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
5123
5124 while(t->det > 0 )
5125 {
5126 assert( kkkk++ < 2000 );
5127 j= OppositeVertex[jj];
5128
5129 #ifdef DRAWING1
5130 t->Draw();
5131 #endif
5132 dete[j] = detop; //det(*b,*s1,*s2);
5133 jn = NextVertex[j];
5134 jp = PreviousVertex[j];
5135 dete[jp]= det(*(*t)(j),*(*t)(jn),B);
5136 dete[jn] = t->det-dete[j] -dete[jp];
5137
5138 #ifdef DEBUG
5139 const Vertex * s0 = (*t)(0);
5140 const Vertex * s1 = (*t)(1);
5141 const Vertex * s2 = (*t)(2);
5142 assert(dete[0] == det(B ,*s1,*s2));
5143 assert(dete[1] == det(*s0,B ,*s2));
5144 assert(dete[2] == det(*s0,*s1,B ));
5145 assert(t->det== (dete[0] + dete[1] +dete[2]));
5146 #endif
5147 // count the number k of dete <0
5148 int k=0,ii[3];
5149 if (dete[0] < 0 ) ii[k++]=0;
5150 if (dete[1] < 0 ) ii[k++]=1;
5151 if (dete[2] < 0 ) ii[k++]=2;
5152 // 0 => ok
5153 // 1 => go in way 1
5154 // 2 => two way go in way 1 or 2 randomly
5155
5156 if (k==0)
5157 break;
5158 if (k == 2 && BinaryRand())
5159 Exchange(ii[0],ii[1]);
5160 assert ( k < 3);
5161 TriangleAdjacent t1 = t->Adj(jj=ii[0]);
5162 if ((t1.det() < 0 ) && (k == 2))
5163 t1 = t->Adj(jj=ii[1]);
5164 t=t1;
5165 j=t1;// for optimisation we now the -det[OppositeVertex[j]];
5166 detop = -dete[OppositeVertex[jj]];
5167 jj = j;
5168 }
5169
5170 if (t->det<0) // outside triangle
5171 dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
5172 // NbOfTriangleSearchFind += kkkk;
5173 return t;
5174 }
5175
5176 }
5177
5178