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

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

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