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

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

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