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

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

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