source: issm/trunk/src/c/Bamgx/Mesh2.h@ 3230

Last change on this file since 3230 was 3230, checked in by Mathieu Morlighem, 15 years ago

cosmetics
File size: 41.5 KB
Line 
1#ifndef _MESH2_H_
2#define _MESH2_H_
3
4#include "../objects/objects.h"
5#include "../shared/shared.h"
6#include "../include/macros.h"
7#include "../toolkits/toolkits.h"
8
9#include "meshtype.h"
10#include "Metric.h"
11
12namespace bamg {
13
14 //classes
15 class Geometry;
16 class Triangles;
17 class Triangle;
18 class QuadTree;
19 class GeometricalEdge;
20 class VertexOnGeom;
21 class VertexOnEdge;
22
23 /*CLASS: DoubleAndInt4 {{{1*/
24 class DoubleAndInt4 {
25 //class used by Triangles::MakeQuadrangles
26
27 public:
28 double q;
29 Int4 i3j;
30
31 //Operators
32 int operator<(DoubleAndInt4 a){return q > a.q;}
33 };
34 /*}}}1*/
35 /*CLASS: Direction{{{1*/
36 class Direction { //
37 private:
38 Icoor1 dir;
39
40 public:
41 //Methods
42 Direction(): dir(MaxICoor){}; // no direction set
43 Direction(Icoor1 i,Icoor1 j) {
44 Icoor2 n2 = 2*(Abs(i)+Abs(j));
45 Icoor2 r = MaxICoor* (Icoor2) i;
46 Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number
47 dir = (j>0) ? r1 : r1+1; // odd-> j>0 even-> j<0
48 }
49 int sens(Icoor1 i,Icoor1 j) {
50 int r =1;
51 if (dir!= MaxICoor) {
52 Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1);
53 r = (x*i + y*j) >=0;
54 }
55 return r;
56 }
57 };
58 /*}}}1*/
59 /*CLASS: Vertex{{{1*/
60 class Vertex {
61
62 public:
63 I2 i; // integer coordinates
64 R2 r; // real coordinates
65 Metric m;
66 Int4 ReferenceNumber;
67 Direction DirOfSearch;
68 Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t
69 union {
70 Triangle* t; // one triangle which is containing the vertex
71 Int4 color;
72 Vertex* to; // use in geometry Vertex to now the Mesh Vertex associed
73 VertexOnGeom* onGeometry; // if vint == 8; // set with Triangles::SetVertexFieldOn()
74 Vertex* onBackgroundVertex; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
75 VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge
76 };
77
78 //Operators
79 operator I2() const {return i;} // Cast operator
80 operator const R2 & () const {return r;} // Cast operator
81 operator Metric () const {return m;} // Cast operator
82 Real8 operator()(R2 x) const { return m(x);} // Get x in the metric m
83
84 //methods (No constructor and no destructors...)
85 Int4 Optim(int =1,int =0);
86 Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1);
87 void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
88 void Echo();
89 int ref() const { return ReferenceNumber;}
90
91 //inline functions
92 inline void Set(const Vertex & rec,const Triangles &,Triangles &);
93 };
94 /*}}}1*/
95 inline Vertex* TheVertex(Vertex * a); // for remove crak in mesh
96 double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &);
97 /*CLASS: TriangleAdjacent{{{1*/
98 class TriangleAdjacent {
99
100 public:
101 Triangle* t; //pointer toward the triangle
102 int a; //Edge number
103
104 //Constructors
105 TriangleAdjacent() {};
106 TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {};
107
108 //Operators
109 operator Triangle * () const {return t;}
110 operator Triangle & () const {return *t;}
111 operator int() const {return a;}
112 TriangleAdjacent & operator++(){ a= NextEdge[a]; return *this; }
113 TriangleAdjacent operator--(){ a= PreviousEdge[a]; return *this; }
114
115 //Methods
116 int swap();
117
118 //Inline methods
119 inline TriangleAdjacent Adj() const ;
120 inline void SetAdj2(const TriangleAdjacent& , int =0);
121 inline Vertex * EdgeVertex(const int &) const ;
122 inline Vertex * OppositeVertex() const ;
123 inline Icoor2 & det() const;
124 inline int Locked() const ;
125 inline int GetAllFlag_UnSwap() const ;
126 inline void SetLock();
127 inline int MarkUnSwap() const;
128 inline void SetMarkUnSwap();
129 inline void SetCracked();
130 inline int Cracked() const ;
131 };
132 /*}}}1*/
133 /*CLASS: Edge{{{1*/
134 class Edge {
135
136 public:
137 Vertex* v[2];
138 Int4 ref;
139 GeometricalEdge* onGeometry;
140 Edge* adj[2]; // the 2 adj edges if on the same curve
141
142 //Operators
143 Vertex & operator[](int i){return *v[i];};
144 Vertex * operator()(int i){return v[i];};
145 R2 operator()(double t) const; // return the point
146 const Vertex & operator[](int i) const { return *v[i];};
147
148 //Methods
149 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){
150 if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
151 if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
152 }
153 int Intersection(const Edge & e) const {
154 if (!(adj[0]==&e || adj[1]==&e)){
155 throw ErrorException(__FUNCT__,exprintf("Intersection bug"));
156 }
157 if (adj[0]!=&e && adj[1]!=&e){
158 throw ErrorException(__FUNCT__,exprintf("adj[0]!=&e && adj[1]!=&e"));
159 }
160 return adj[0]==&e ? 0 : 1;
161 }
162
163 //Inline methods
164 inline void Set(const Triangles &,Int4,Triangles &);
165
166 };
167 /*}}}1*/
168 inline TriangleAdjacent FindTriangleAdjacent(Edge &E);
169 /*CLASS: GeometricalVertex{{{1*/
170 class GeometricalVertex : public Vertex {
171
172 public:
173 friend class Geometry;
174
175 int cas;
176 GeometricalVertex* link; // link all the same GeometricalVertex circular (Crack)
177
178 //Methods
179 int Corner() const {return cas&4;}
180 int Required()const {return cas&6;}// a corner is required
181 int IsThe() const { return link == this;}
182 void SetCorner(){ cas |= 4;}
183 void SetRequired(){ cas |= 2;}
184 void Set(){cas=0;}
185 GeometricalVertex() :cas(0), link(this) {};
186 GeometricalVertex* The() {
187 if (!link){ throw ErrorException(__FUNCT__,exprintf("!link"));}
188 return link;
189 }// return a unique vertex
190
191 //Inline methods
192 inline void Set(const GeometricalVertex & rec,const Geometry & Gh ,const Geometry & GhNew);
193 };
194 /*}}}1*/
195 /*CLASS: GeometricalEdge{{{1*/
196 class GeometricalEdge {
197
198 public:
199 GeometricalVertex* v[2];
200 Int4 ref;
201 Int4 CurveNumber;
202 R2 tg[2]; // the 2 tangentes (tg[0] =0 => no continuity)
203 GeometricalEdge* Adj[2];
204 int DirAdj[2];
205 int flag ;
206 GeometricalEdge* link; // if Cracked() or Equi()
207
208 //Operators
209 GeometricalVertex & operator[](int i){return *v[i];};
210 const GeometricalVertex & operator[](int i) const { return *v[i];};
211 GeometricalVertex * operator()(int i){return v[i];};
212
213 //Methods
214 R2 F(Real8 theta) const ; // parametrization of the curve edge
215 Real8 R1tg(Real8 theta,R2 &t) const ; // 1/radius of curvature + tangente
216 int Cracked() const {return flag & 1;}
217 int Dup() const { return flag & 32;}
218 int Equi()const {return flag & 2;}
219 int ReverseEqui()const {return flag & 128;}
220 int TgA()const {return flag &4;}
221 int TgB()const {return flag &8;}
222 int Tg(int i) const{return i==0 ? TgA() : TgB();}
223 int Mark()const {return flag &16;}
224 int Required() { return flag & 64;}
225 void SetCracked() { flag |= 1;}
226 void SetDup() { flag |= 32;} // not a real edge
227 void SetEqui() { flag |= 2;}
228 void SetTgA() { flag|=4;}
229 void SetTgB() { flag|=8;}
230 void SetMark() { flag|=16;}
231 void SetUnMark() { flag &= 1007 /* 1023-16*/;}
232 void SetRequired() { flag |= 64;}
233 void SetReverseEqui() {flag |= 128;}
234
235 //Inline methods
236 inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew);
237 };
238 /*}}}1*/
239 /*CLASS: Curve{{{1*/
240 class Curve {
241 public:
242 GeometricalEdge * be,*ee; // begin et end edge
243 int kb,ke; // begin vetex and end vertex
244 Curve *next; // next curve equi to this
245 bool master; // true => of equi curve point on this curve
246
247 //Methods
248 Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {}
249 void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse
250
251 //Inline methods
252 inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
253 };
254 /*}}}1*/
255 /*CLASS: Triangle{{{1*/
256 class Triangle {
257
258 friend class TriangleAdjacent;
259
260 private:
261 Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
262 Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles
263 Int1 TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1
264
265 public:
266 Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
267 union {
268 Triangle * link ;
269 Int4 color;
270 };
271
272 //Constructors/Destructors
273 Triangle() {}
274 Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k);
275 Triangle(Vertex *v0,Vertex *v1,Vertex *v2);
276
277 //Operators
278 const Vertex & operator[](int i) const {return *TriaVertices[i];};
279 Vertex & operator[](int i) {return *TriaVertices[i];};
280 const Vertex * operator()(int i) const {return TriaVertices[i];};
281 Vertex * & operator()(int i) {return TriaVertices[i];};
282
283 //Methods
284 void Echo();
285 int swap(Int2 a1,int=0);
286 Int4 Optim(Int2 a,int =0);
287 int Locked(int a)const { return TriaAdjSharedEdge[a]&4;}
288 int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;}
289 int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;}
290 int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;}
291 void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);}
292 double QualityQuad(int a,int option=1) const;
293 Int1 NuEdgeTriangleAdj(int i) const {return TriaAdjSharedEdge[i&3]&3;} // Number of the adjacent edge in adj tria
294 TriangleAdjacent FindBoundaryEdge(int i) const;
295 TriangleAdjacent Adj(int i) const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};
296 Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
297 Triangle* Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
298 void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){
299 if (link >=tb && link <te) link = tb + renu[link -tb];
300 if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];
301 if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb];
302 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];
303 }
304 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){
305 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
306 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
307 if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb];
308 }
309 void SetAdjAdj(Int1 a){
310 a &= 3;
311 register Triangle *tt=TriaAdjTriangles[a];
312 TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap
313 register Int1 aatt = TriaAdjSharedEdge[a] & 3;
314 if(tt){
315 tt->TriaAdjTriangles[aatt]=this;
316 tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark
317 }
318 void SetAdj2(Int1 a,Triangle *t,Int1 aat){
319 TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t
320 TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle
321 if(t) { //if t!=NULL add adjacent triangle to t (this)
322 t->TriaAdjTriangles[aat]=this;
323 t->TriaAdjSharedEdge[aat]=a;
324 }
325 }
326 void SetTriangleContainingTheVertex() {
327 if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);
328 if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);
329 if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2);
330 }
331 void SetHidden(int a){
332 //Get Adjacent Triangle number a
333 register Triangle* t = TriaAdjTriangles[a];
334 //if it exist
335 //C|=D -> C=(C|D) bitwise inclusive OR
336 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16;
337 TriaAdjSharedEdge[a] |= 16;
338 }
339 void SetCracked(int a){
340 register Triangle * t = TriaAdjTriangles[a];
341 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32;
342 TriaAdjSharedEdge[a] |= 32;
343 }
344
345 void SetLocked(int a){
346 //mark the edge as on Boundary
347 register Triangle * t = TriaAdjTriangles[a];
348 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4;
349 TriaAdjSharedEdge[a] |= 4;
350 }
351 void SetMarkUnSwap(int a){
352 register Triangle * t = TriaAdjTriangles[a];
353 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8;
354 TriaAdjSharedEdge[a] |=8 ;
355 }
356 void SetUnMarkUnSwap(int a){
357 register Triangle * t = TriaAdjTriangles[a];
358 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32
359 TriaAdjSharedEdge[a] &=55 ;
360 }
361 void SetDet() {
362 if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
363 else det = -1; }
364
365 //Inline methods
366 inline Real4 qualite() ;
367 inline void Set(const Triangle &,const Triangles &,Triangles &);
368 inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
369
370
371 }; // end of Triangle class
372 /*}}}1*/
373 /*CLASS: ListofIntersectionTriangles{{{1*/
374 class ListofIntersectionTriangles {
375
376 class IntersectionTriangles {
377 public:
378 Triangle* t;
379 Real8 bary[3]; // use if t != 0
380 R2 x;
381 Metric m;
382 Real8 s; // curvilinear coordinate
383 Real8 sp;// length of the previous segment in m
384 Real8 sn;// length of the next segment in m
385 };
386
387 class SegInterpolation {
388 public:
389 GeometricalEdge * e;
390 Real8 sBegin,sEnd; // abscisse of the seg on edge parameter
391 Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
392 int last;// last index in ListofIntersectionTriangles for this Sub seg of edge
393
394 //Methods
395 R2 F(Real8 s){
396 Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01;
397 if (lBegin>s || s>lEnd){
398 throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd"));
399 }
400 return e->F(sBegin*c0+sEnd*c1);
401 }
402 };
403
404 public:
405
406 int MaxSize;
407 int Size;
408 Real8 len;
409 int state;
410 IntersectionTriangles * lIntTria;
411 int NbSeg;
412 int MaxNbSeg;
413 SegInterpolation * lSegsI;
414
415 //Constructors/Destructors
416 ListofIntersectionTriangles(int n=256,int m=16)
417 : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
418 NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
419 long int verbosity=0;
420 if (verbosity>9) printf(" construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg);
421 };
422 ~ListofIntersectionTriangles(){
423 if (lIntTria) delete [] lIntTria,lIntTria=0;
424 if (lSegsI) delete [] lSegsI,lSegsI=0;
425 }
426
427 //Operators
428 IntersectionTriangles & operator[](int i) {return lIntTria[i];}
429 operator int&() {return Size;}
430
431 //Methods
432 void init(){state=0;len=0;Size=0;}
433 int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2);
434 int NewItem(R2,const Metric & );
435 void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0);
436 Real8 Length();
437 Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx);
438 void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1){
439 long int verbosity=0;
440 if (NbSeg>=MaxNbSeg) {
441 int mneo= MaxNbSeg;
442 MaxNbSeg *= 2;
443 if (verbosity>3){
444 printf(" reshape lSegsI from %i to %i\n",mneo,MaxNbSeg);
445 }
446 SegInterpolation * lEn = new SegInterpolation[MaxNbSeg];
447 if (!lSegsI || NbSeg>=MaxNbSeg){
448 throw ErrorException(__FUNCT__,exprintf("!lSegsI || NbSeg>=MaxNbSeg"));
449 }
450 for (int i=0;i< NbSeg;i++)
451 lEn[i] = lSegsI[MaxNbSeg]; // copy old to new
452 delete [] lSegsI; // remove old
453 lSegsI = lEn;
454 }
455 if (NbSeg) lSegsI[NbSeg-1].last=Size;
456 lSegsI[NbSeg].e=e;
457 lSegsI[NbSeg].sBegin=s0;
458 lSegsI[NbSeg].sEnd=s1;
459 NbSeg++;
460 }
461 void ReShape(){
462 register int newsize = MaxSize*2;
463 IntersectionTriangles* nw = new IntersectionTriangles[newsize];
464 if (!nw){ throw ErrorException(__FUNCT__,exprintf("!nw"));}
465 // recopy
466 for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i];
467 long int verbosity=0;
468 if(verbosity>3) printf(" ListofIntersectionTriangles ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg);
469 MaxSize = newsize;
470 delete [] lIntTria;// remove old
471 lIntTria = nw; // copy pointer
472 }
473
474 };
475 /*}}}1*/
476 /*CLASS: GeometricalSubDomain{{{1*/
477 class GeometricalSubDomain {
478 public:
479 GeometricalEdge *edge;
480 int sens; // -1 or 1
481 Int4 ref;
482
483 //Inline methods
484 inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &);
485
486 };
487 /*}}}1*/
488 /*CLASS: SubDomain{{{1*/
489 class SubDomain {
490 public:
491 Triangle * head;
492 Int4 ref;
493 int sens; // -1 or 1
494 Edge * edge; // to geometrical
495
496 //Inline methods
497 inline void Set(const Triangles &,Int4,Triangles &);
498 };
499 /*}}}1*/
500 /*CLASS: VertexOnGeom{{{1*/
501 class VertexOnGeom{
502
503 public:
504
505 Vertex* mv;
506 Real8 abscisse;
507 union{
508 GeometricalVertex * gv; // if abscisse <0;
509 GeometricalEdge * ge; // if abscisse in [0..1]
510 };
511
512 //Constructors/Destructors
513 VertexOnGeom(): mv(0),abscisse(0){gv=0;}
514 VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
515 VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;}
516
517 //Operators
518 operator Vertex * () const {return mv;}
519 operator GeometricalVertex * () const {return gv;}
520 operator GeometricalEdge * () const {return ge;}
521 operator const Real8 & () const {return abscisse;}
522
523 //Methods
524 int OnGeomVertex()const {return this? abscisse <0 :0;}
525 int OnGeomEdge() const {return this? abscisse >=0 :0;}
526 int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;}
527 void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;}
528
529 //Inline methods
530 inline void Set(const Triangles &,Int4,Triangles &);
531 inline void Set(const VertexOnGeom&,const Triangles &,Triangles &);
532
533 };
534 /*}}}1*/
535 /*CLASS: VertexOnVertex{{{1*/
536 class VertexOnVertex {
537
538 public:
539 Vertex* v;
540 Vertex* bv;
541
542 //Constructors
543 VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){}
544 VertexOnVertex() {};
545
546 //Methods
547 void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}
548
549 //Inline methods
550 inline void Set(const Triangles &,Int4,Triangles &);
551 };
552 /*}}}1*/
553 /*CLASS: VertexOnEdge{{{1*/
554 class VertexOnEdge {
555
556 public:
557 Vertex* v;
558 Edge* be;
559 Real8 abcisse;
560
561 //Constructors
562 VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {}
563 VertexOnEdge(){}
564
565 //Operators
566 operator Real8 () const { return abcisse;}
567 operator Vertex* () const { return v;}
568 operator Edge* () const { return be;}
569 Vertex & operator[](int i) const { return (*be)[i];}
570
571 //Methods
572 void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;}
573
574 //Inline methods
575 inline void Set(const Triangles &,Int4,Triangles &);
576 };
577 /*}}}1*/
578 /*CLASS: CrackedEdge{{{1*/
579 class CrackedEdge {
580
581 friend class Triangles;
582
583 class CrackedTriangle {
584 friend class Triangles;
585 friend class CrackedEdge;
586 Triangle* t; // edge of triangle t
587 int i; // edge number of in triangle
588 Edge *edge; // the 2 edge
589 Vertex *New[2]; // new vertex number
590
591 //Constructors
592 CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;}
593 CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;}
594
595 //Methods
596 void Crack(){
597 Triangle & T(*t);
598 int i0=VerticesOfTriangularEdge[i][0];
599 int i1=VerticesOfTriangularEdge[i][0];
600 if (!New[0] && !New[1]){
601 throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
602 }
603 T(i0) = New[0];
604 T(i1) = New[1];
605 }
606 void UnCrack(){
607 Triangle & T(*t);
608 int i0=VerticesOfTriangularEdge[i][0];
609 int i1=VerticesOfTriangularEdge[i][0];
610 if (!New[0] && !New[1]){
611 throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
612 }
613 T(i0) = TheVertex(T(i0));
614 T(i1) = TheVertex(T(i1));
615 }
616 void Set() {
617 TriangleAdjacent ta ( FindTriangleAdjacent(*edge));
618 t = ta;
619 i = ta;
620
621 New[0] = ta.EdgeVertex(0);
622 New[1] = ta.EdgeVertex(1);
623 // warning the ref
624 }
625 };
626
627 public:
628 CrackedTriangle a,b;
629
630 //Constructors
631 CrackedEdge() :a(),b() {}
632 CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {};
633 CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {};
634
635 //Methods
636 void Crack() { a.Crack(); b.Crack();}
637 void UnCrack() { a.UnCrack(); b.UnCrack();}
638 void Set() { a.Set(), b.Set();}
639 };
640 /*}}}1*/
641 /*CLASS: Triangles{{{1*/
642 class Triangles {
643 public:
644
645 Geometry & Gh; // Geometry
646 Triangles & BTh; // Background Mesh Bth==*this =>no background
647 Int4 NbRef; // counter of ref on the this class if 0 we can delete
648 Int4 nbvx,nbtx; // nombre max de sommets , de triangles
649 Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference,
650 Int4 NbOfQuad; // nb of quadrangle
651 Int4 NbSubDomains;
652 Int4 NbOutT; // Nb of oudeside triangle
653 Int4 NbOfTriangleSearchFind;
654 Int4 NbOfSwapTriangle;
655 Vertex* vertices;
656 Int4 NbVerticesOnGeomVertex;
657 VertexOnGeom * VerticesOnGeomVertex;
658 Int4 NbVerticesOnGeomEdge;
659 VertexOnGeom * VerticesOnGeomEdge;
660 Int4 NbVertexOnBThVertex;
661 VertexOnVertex *VertexOnBThVertex;
662 Int4 NbVertexOnBThEdge;
663 VertexOnEdge *VertexOnBThEdge;
664 Int4 NbCrackedVertices;
665 Int4 NbCrackedEdges;
666 CrackedEdge *CrackedEdges;
667 R2 pmin,pmax; // extrema
668 Real8 coefIcoor; // coef to integer Icoor1;
669 Triangle* triangles;
670 Edge* edges;
671 QuadTree *quadtree;
672 Vertex** ordre;
673 SubDomain* subdomains;
674 ListofIntersectionTriangles lIntTria;
675
676 //Constructors/Destructors
677 Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts);
678 Triangles(double* index,double* x,double* y,int nods,int nels);
679 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); //copy operator
680 Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
681 Triangles(Int4 nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
682 try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}
683 catch(...) { this->~Triangles(); throw; }
684 }
685 Triangles(Int4 nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
686 try { GeomToTriangles0(nbvx,bamgopts);}
687 catch(...) { this->~Triangles(); throw; }
688 }
689 ~Triangles();
690
691 //Operators
692 const Vertex & operator[] (Int4 i) const { return vertices[i];};
693 Vertex & operator[](Int4 i) { return vertices[i];};
694 const Triangle & operator() (Int4 i) const { return triangles[i];};
695 Triangle & operator()(Int4 i) { return triangles[i];};
696
697 //Methods
698 void SetIntCoor(const char * from =0);
699 Real8 MinimalHmin() {return 2.0/coefIcoor;}
700 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
701 I2 toI2(const R2 & P) const {
702 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
703 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
704 R2 toR2(const I2 & P) const {
705 return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
706 }
707 void Add( Vertex & s,Triangle * t,Icoor2 * =0) ;
708 void Insert();
709 void ForceBoundary();
710 void Renumber(BamgOpts* bamgopts);
711 void FindSubDomain(int OutSide=0);
712 Int4 ConsRefTriangle(Int4 *) const;
713 void ShowHistogram() const;
714 void ShowRegulaty() const;
715 void ReMakeTriangleContainingTheVertex();
716 void UnMarkUnSwapTriangle();
717 void SmoothMetric(Real8 raisonmax) ;
718 void BoundAnisotropy(Real8 anisomax,double hminaniso= 1e-100) ;
719 void MaxSubDivision(Real8 maxsubdiv);
720 Edge** MakeGeometricalEdgeToEdge();
721 void SetVertexFieldOn();
722 void SetVertexFieldOnBTh();
723 Int4 SplitInternalEdgeWithBorderVertices();
724 void MakeQuadrangles(double costheta);
725 int SplitElement(int choice);
726 void MakeQuadTree();
727 void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1);
728 Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ;
729 void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
730 void ReNumberingVertex(Int4 * renu);
731 void SmoothingVertex(int =3,Real8=0.3);
732 Metric MetricAt (const R2 &) const;
733 GeometricalEdge* ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta, Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
734 Int4 Number(const Triangle & t) const { return &t - triangles;}
735 Int4 Number(const Triangle * t) const { return t - triangles;}
736 Int4 Number(const Vertex & t) const { return &t - vertices;}
737 Int4 Number(const Vertex * t) const { return t - vertices;}
738 Int4 Number(const Edge & t) const { return &t - edges;}
739 Int4 Number(const Edge * t) const { return t - edges;}
740 Int4 Number2(const Triangle * t) const {
741 return t - triangles;
742 }
743 Vertex* NearestVertex(Icoor1 i,Icoor1 j) ;
744 Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
745 void ReadMesh(double* index,double* x,double* y,int nods,int nels);
746 void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
747 void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
748 void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);
749 void WriteMetric(BamgOpts* bamgopts);
750 void IntersectConsMetric(BamgOpts* bamgopts);
751 void BuildMetric0(BamgOpts* bamgopts);
752 void BuildMetric1(BamgOpts* bamgopts);
753 void IntersectGeomMetric(BamgOpts* bamgopts);
754 int isCracked() const {return NbCrackedVertices != 0;}
755 int Crack();
756 int UnCrack();
757 void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
758 void FillHoleInMesh() ;
759 int CrackMesh();
760
761 private:
762 void GeomToTriangles1(Int4 nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
763 void GeomToTriangles0(Int4 nbvx,BamgOpts* bamgopts);// the real constructor mesh generator
764 void PreInit(Int4);
765
766 };
767 /*}}}1*/
768 /*CLASS: Geometry{{{1*/
769 class Geometry {
770
771 public:
772 Int4 NbRef; // counter of ref on the this class if 0 we can delete
773 Int4 nbvx,nbtx; // maximum number of vertices
774 Int4 nbv,nbt,nbiv,nbe; // number of vertices
775 Int4 NbSubDomains; //
776 Int4 NbOfCurves;
777 GeometricalVertex* vertices;
778 Triangle* triangles;
779 GeometricalEdge* edges;
780 QuadTree* quadtree;
781 GeometricalSubDomain* subdomains;
782 Curve* curves;
783 R2 pmin,pmax; // extrema
784 Real8 coefIcoor; // coef to integer Icoor1;
785 Real8 MaximalAngleOfCorner;
786
787 //Constructor/Destructors
788 ~Geometry();
789 Geometry(const Geometry & Gh); //Copy Operator
790 Geometry(int nbg,const Geometry** ag); // intersection operator
791
792 //Operators
793 const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];};
794 GeometricalVertex & operator[](Int4 i) { return vertices[i];};
795 const GeometricalEdge & operator() (Int4 i) const { return edges[i];};
796 GeometricalEdge & operator()(Int4 i) { return edges[i];};
797
798 //Methods
799 int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }
800 void Echo();
801 I2 toI2(const R2 & P) const {
802 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
803 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
804 Real8 MinimalHmin() {return 2.0/coefIcoor;}
805 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
806 void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
807 void EmptyGeometry();
808 Geometry() {EmptyGeometry();}// empty Geometry
809 void AfterRead();
810 Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();}
811 Int4 Number(const GeometricalVertex & t) const { return &t - vertices;}
812 Int4 Number(const GeometricalVertex * t) const { return t - vertices;}
813 Int4 Number(const GeometricalEdge & t) const { return &t - edges;}
814 Int4 Number(const GeometricalEdge * t) const { return t - edges;}
815 Int4 Number(const Curve * c) const { return c - curves;}
816 void UnMarkEdges() {for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();}
817 GeometricalEdge * ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ;
818 GeometricalEdge * Contening(const R2 P, GeometricalEdge * start) const;
819 void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
820 };
821 /*}}}1*/
822
823 /*INLINE functions{{{1*/
824 // to sort in descending order
825 template<class T> inline void HeapSort(T *c,long n){
826 /*Intermediary*/
827 int l,j,r,i;
828 T crit;
829 c--; //the array must starts at 1 and not 0
830 if(n<=1) return; //return if size <=1
831 l=n/2+1; //initialize l and r
832 r=n;
833 for(;;){
834 if(l<=1){
835 crit =c[r];
836 c[r--]=c[1];
837 if (r==1){c[1]=crit; return;}
838 }
839 else crit = c[--l];
840 j=l;
841 for(;;){
842 i=j;
843 j=2*j;
844 if (j>r) {c[i]=crit;break;}
845 if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
846 if (crit < c[j]) c[i]=c[j]; //c[j] > crit -> stock this large value in i(<j)
847 else{c[i]=crit;break;} //c[j] < crit -> stock crit in i (<j)
848 }
849 }
850 }
851 // to sort in descending order and return ordering
852 template<class T> inline void HeapSort(int** porder,T* c,int n){
853 /*Intermediary*/
854 int l,j,r,i;
855 T crit;
856 int pos;
857 int* order=NULL;
858 order=(int*)xmalloc(n*sizeof(int));
859 for(i=0;i<n;i++) order[i]=i+1;
860 c--; //the array must starts at 1 and not 0
861 order--;
862 if(n<=1) return; //return if size <=1
863 l=n/2+1; //initialize l and r
864 r=n;
865 for(;;){
866 if(l<=1){
867 crit =c[r]; pos=order[r];
868 c[r--]=c[1]; order[r+1]=order[1];
869 if (r==1){
870 c[1]=crit; order[1]=pos;
871 order++;
872 *porder=order;
873 return;
874 }
875 }
876 else {crit=c[--l]; pos=order[l];}
877 j=l;
878 for(;;){
879 i=j;
880 j=2*j;
881 if (j>r) {c[i]=crit;order[i]=pos;break;}
882 if ((j<r) && (c[j] < c[j+1]))j++;
883 if (crit < c[j]) {c[i]=c[j];order[i]=order[j];}
884 else{c[i]=crit;order[i]=pos;break;}
885 }
886 }
887 }
888 inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){
889 return A[0]*( B[1]*C[2]-B[2]*C[1])
890 - A[1]*( B[0]*C[2]-B[2]*C[0])
891 + A[2]*( B[0]*C[1]-B[1]*C[0]);
892 }
893 inline TriangleAdjacent Adj(const TriangleAdjacent & a)
894 { return a.Adj();}
895
896 inline TriangleAdjacent Next(const TriangleAdjacent & ta)
897 { return TriangleAdjacent(ta.t,NextEdge[ta.a]);}
898
899 inline TriangleAdjacent Previous(const TriangleAdjacent & ta)
900 { return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);}
901 inline void Adj(GeometricalEdge * & on,int &i)
902 {int j=i;i=on->DirAdj[i];on=on->Adj[j];}
903 inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc)
904 {
905 Real4 ret;
906 I2 ia=va,ib=vb,ic=vc;
907 I2 ab=ib-ia,bc=ic-ib,ac=ic-ia;
908 Icoor2 deta=Det(ab,ac);
909 if (deta <=0) ret = -1;
910 else {
911 Real8 a = sqrt((Real8) (ac,ac)),
912 b = sqrt((Real8) (bc,bc)),
913 c = sqrt((Real8) (ab,ab)),
914 p = a+b+c;
915 Real8 h= Max(Max(a,b),c),ro=deta/p;
916 ret = ro/h;}
917 return ret;
918 }
919 Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c){
920 register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y;
921 register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y;
922 return bax*cay - bay*cax;
923 }
924 inline TriangleAdjacent FindTriangleAdjacent(Edge &E){
925 Vertex * a = E.v[0];
926 Vertex * b = E.v[1];
927
928 Triangle * t = a->t;
929 int i = a->vint;
930 TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
931 if (!t || i<0 || i>=3){
932 throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
933 }
934 if ( a!=(*t)(i)){
935 throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
936 }
937 int k=0;
938 do { // turn around vertex in direct sens (trigo)
939 k++;
940 if (k>=20000){
941 throw ErrorException(__FUNCT__,exprintf("k>=20000"));
942 }
943 // in no crack => ta.EdgeVertex(1) == a otherwise ???
944 if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find
945 ta = ta.Adj();
946 if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find
947 --ta;
948 } while (t != (Triangle *)ta);
949 throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found"));
950 return TriangleAdjacent(0,0);//for compiler
951 }
952 inline Vertex* TheVertex(Vertex * a){
953 // give a unique vertex with smallest number
954 // in case on crack in mesh
955 Vertex * r(a), *rr;
956 Triangle * t = a->t;
957 int i = a->vint;
958 TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
959 if (!t || i<0 || i>=3){
960 throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
961 }
962 if ( a!=(*t)(i)){
963 throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
964 }
965 int k=0;
966 do { // turn around vertex in direct sens (trigo)
967 k++;
968 if (k>=20000){
969 throw ErrorException(__FUNCT__,exprintf("k>=20000"));
970 }
971 // in no crack => ta.EdgeVertex(1) == a
972 if ((rr=ta.EdgeVertex(0)) < r) r = rr;
973 ta = ta.Adj();
974 if ((rr=ta.EdgeVertex(1)) < r) r =rr;
975 --ta;
976 } while (t != (Triangle*) ta);
977 return r;
978 }
979
980 /*}}}1*/
981
982 /*INLINE functions of CLASS Vertex{{{1*/
983 inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){
984 *this = rec;
985 }
986 Int4 inline Vertex::Optim(int i,int koption){
987 Int4 ret=0;
988 if ( t && (vint >= 0 ) && (vint <3) )
989 {
990 ret = t->Optim(vint,koption);
991 if(!i)
992 {
993 t =0; // for no future optime
994 vint= 0; }
995 }
996 return ret;
997 }
998 /*}}}1*/
999 /*INLINE functions of CLASS GeometricalVertex{{{1*/
1000 inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){
1001 *this = rec;
1002 }
1003 /*}}}1*/
1004 /*INLINE functions of CLASS VertexOnVertex{{{1*/
1005 inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) {
1006 *this = Th.VertexOnBThVertex[i];
1007 v = ThNew.vertices + Th.Number(v);
1008
1009 }
1010/*}}}1*/
1011 /*INLINE functions of CLASS Triangles{{{1*/
1012 inline void Triangles::ReMakeTriangleContainingTheVertex(){
1013 register Int4 i;
1014 for ( i=0;i<nbv;i++){
1015 vertices[i].vint=0;
1016 vertices[i].t=0;
1017 }
1018 for ( i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
1019 }
1020
1021 inline void Triangles::UnMarkUnSwapTriangle()
1022 {
1023 register Int4 i;
1024 for ( i=0;i<nbt;i++)
1025 for(int j=0;j<3;j++)
1026 triangles[i].SetUnMarkUnSwap(j);
1027 }
1028
1029 inline void Triangles::SetVertexFieldOn(){
1030 for (Int4 i=0;i<nbv;i++)
1031 vertices[i].onGeometry=0;
1032 for (Int4 j=0;j<NbVerticesOnGeomVertex;j++ )
1033 VerticesOnGeomVertex[j].SetOn();
1034 for (Int4 k=0;k<NbVerticesOnGeomEdge;k++ )
1035 VerticesOnGeomEdge[k].SetOn();
1036 }
1037 inline void Triangles::SetVertexFieldOnBTh(){
1038 for (Int4 i=0;i<nbv;i++)
1039 vertices[i].onGeometry=0;
1040 for (Int4 j=0;j<NbVertexOnBThVertex;j++ )
1041 VertexOnBThVertex[j].SetOnBTh();
1042 for (Int4 k=0;k<NbVertexOnBThEdge;k++ )
1043 VertexOnBThEdge[k].SetOnBTh();
1044
1045 }
1046 /*}}}1*/
1047 /*INLINE functions of CLASS Triangle{{{1*/
1048 inline Triangle* Triangle::Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const
1049 {
1050 // return the other triangle of the quad if a quad or 0 if not a quat
1051 Triangle * t =0;
1052 if (link) {
1053 int a=-1;
1054 if (TriaAdjSharedEdge[0] & 16 ) a=0;
1055 if (TriaAdjSharedEdge[1] & 16 ) a=1;
1056 if (TriaAdjSharedEdge[2] & 16 ) a=2;
1057 if (a>=0) {
1058 t = TriaAdjTriangles[a];
1059 // if (t-this<0) return 0;
1060 v2 = TriaVertices[VerticesOfTriangularEdge[a][0]];
1061 v0 = TriaVertices[VerticesOfTriangularEdge[a][1]];
1062 v1 = TriaVertices[OppositeEdge[a]];
1063 v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
1064 }
1065 }
1066 return t;
1067 }
1068
1069 inline double Triangle::QualityQuad(int a,int option) const
1070 { // first do the logique part
1071 double q;
1072 if (!link || TriaAdjSharedEdge[a] &4)
1073 q= -1;
1074 else {
1075 Triangle * t = TriaAdjTriangles[a];
1076 if (t-this<0) q= -1;// because we do 2 times
1077 else if (!t->link ) q= -1;
1078 else if (TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16 || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 )
1079 q= -1;
1080 else if(option)
1081 {
1082 const Vertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
1083 const Vertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
1084 const Vertex & v1 = *TriaVertices[OppositeEdge[a]];
1085 const Vertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
1086 q = QuadQuality(v0,v1,v2,v3); // do the float part
1087 }
1088 else q= 1;
1089 }
1090 return q;
1091 }
1092 inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew)
1093 {
1094 *this = rec;
1095 if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices + Th.Number(TriaVertices[0]);
1096 if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices + Th.Number(TriaVertices[1]);
1097 if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices + Th.Number(TriaVertices[2]);
1098 if(TriaAdjTriangles[0]) TriaAdjTriangles[0] = ThNew.triangles + Th.Number(TriaAdjTriangles[0]);
1099 if(TriaAdjTriangles[1]) TriaAdjTriangles[1] = ThNew.triangles + Th.Number(TriaAdjTriangles[1]);
1100 if(TriaAdjTriangles[2]) TriaAdjTriangles[2] = ThNew.triangles + Th.Number(TriaAdjTriangles[2]);
1101 if (link >= Th.triangles && link < Th.triangles + Th.nbt)
1102 link = ThNew.triangles + Th.Number(link);
1103 }
1104 inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) {
1105 Vertex *v=Th->vertices;
1106 Int4 nbv = Th->nbv;
1107 if (i<0 || j<0 || k<0){
1108 throw ErrorException(__FUNCT__,exprintf("i<0 || j<0 || k<0"));
1109 }
1110 if (i>=nbv || j>=nbv || k>=nbv){
1111 throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv"));
1112 }
1113 TriaVertices[0]=v+i;
1114 TriaVertices[1]=v+j;
1115 TriaVertices[2]=v+k;
1116 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
1117 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
1118 det=0;
1119 }
1120 inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){
1121 TriaVertices[0]=v0;
1122 TriaVertices[1]=v1;
1123 TriaVertices[2]=v2;
1124 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
1125 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
1126 if (v0) det=0;
1127 else {
1128 det=-1;
1129 link=NULL;};
1130 }
1131 inline Real4 Triangle::qualite()
1132 {
1133 return det < 0 ? -1 : bamg::qualite(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
1134 }
1135
1136 /*}}}1*/
1137 /*INLINE functions of CLASS Edge{{{1*/
1138 inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
1139 {
1140 *this = Th.edges[i];
1141 v[0] = ThNew.vertices + Th.Number(v[0]);
1142 v[1] = ThNew.vertices + Th.Number(v[1]);
1143 if (onGeometry)
1144 onGeometry = ThNew.Gh.edges+Th.Gh.Number(onGeometry);
1145 if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]);
1146 if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]);
1147
1148 }
1149
1150 /*}}}1*/
1151 /*INLINE functions of CLASS GeometricalEdge{{{1*/
1152 inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew)
1153 {
1154 *this = rec;
1155 v[0] = GhNew.vertices + Gh.Number(v[0]);
1156 v[1] = GhNew.vertices + Gh.Number(v[1]);
1157 if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]);
1158 if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]);
1159 }
1160 /*}}}1*/
1161 /*INLINE functions of CLASS Curve{{{1*/
1162 inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew)
1163 {
1164 *this = rec;
1165 be = GhNew.edges + Gh.Number(be);
1166 ee = GhNew.edges + Gh.Number(ee);
1167 if(next) next= GhNew.curves + Gh.Number(next);
1168 }
1169 /*}}}1*/
1170 /*INLINE functions of CLASS SubDomain{{{1*/
1171 inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
1172 {
1173 *this = Th.subdomains[i];
1174 if ( head-Th.triangles<0 || head-Th.triangles>=Th.nbt){
1175 throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt"));
1176 }
1177 head = ThNew.triangles + Th.Number(head) ;
1178 if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{
1179 throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe"));
1180 }
1181 edge = ThNew.edges+ Th.Number(edge);
1182 }
1183 /*}}}1*/
1184 /*INLINE functions of CLASS GeometricalSubDomain{{{1*/
1185 inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew)
1186 {
1187 *this = rec;
1188 edge = Gh.Number(edge) + GhNew.edges;
1189 }
1190 /*}}}1*/
1191 /*INLINE functions of CLASS VertexOnEdge{{{1*/
1192 inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
1193 {
1194 *this = Th.VertexOnBThEdge[i];
1195 v = ThNew.vertices + Th.Number(v);
1196 }
1197 /*}}}1*/
1198 /*INLINE functions of CLASS VertexOnGeom{{{1*/
1199 inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew)
1200 {
1201 *this = rec;
1202 mv = ThNew.vertices + Th.Number(mv);
1203 if (gv)
1204 if (abscisse < 0 )
1205 gv = ThNew.Gh.vertices + Th.Gh.Number(gv);
1206 else
1207 ge = ThNew.Gh.edges + Th.Gh.Number(ge);
1208
1209 }
1210 /*}}}1*/
1211 /*INLINE functions of CLASS TriangleAdjacent{{{1*/
1212 inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l ){
1213 //set Adjacent Triangle of a triangle
1214 if(t) {
1215 t->TriaAdjTriangles[a]=ta.t;
1216 t->TriaAdjSharedEdge[a]=ta.a|l;
1217 }
1218 if(ta.t) {
1219 ta.t->TriaAdjTriangles[ta.a] = t ;
1220 ta.t->TriaAdjSharedEdge[ta.a] = a| l ;
1221 }
1222 }
1223 inline int TriangleAdjacent::Locked() const
1224 { return t->TriaAdjSharedEdge[a] &4;}
1225 inline int TriangleAdjacent::Cracked() const
1226 { return t->TriaAdjSharedEdge[a] &32;}
1227 inline int TriangleAdjacent::GetAllFlag_UnSwap() const
1228 { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap
1229 inline int TriangleAdjacent::MarkUnSwap() const
1230 { return t->TriaAdjSharedEdge[a] &8;}
1231 inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);}
1232 inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);}
1233 inline TriangleAdjacent TriangleAdjacent::Adj() const
1234 { return t->Adj(a);}
1235 inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const
1236 {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; }
1237 inline Vertex * TriangleAdjacent::OppositeVertex() const
1238 {return t->TriaVertices[bamg::OppositeVertex[a]]; }
1239 inline Icoor2 & TriangleAdjacent::det() const
1240 { return t->det;}
1241 int inline TriangleAdjacent::swap()
1242 { return t->swap(a);}
1243 /*}}}1*/
1244
1245 /*Other prototypes {{{1*/
1246 TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ;
1247 TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b);
1248 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside);
1249 void swap(Triangle *t1,Int1 a1,
1250 Triangle *t2,Int1 a2,
1251 Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2);
1252 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
1253 TriangleAdjacent & tt1,Icoor2 & dets1,
1254 Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
1255 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ;
1256 /*}}}1*/
1257
1258}
1259#endif
Note: See TracBrowser for help on using the repository browser.