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

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

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