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

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

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