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

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

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