Changeset 2870


Ignore:
Timestamp:
01/20/10 09:07:36 (15 years ago)
Author:
Mathieu Morlighem
Message:

cosmetics

Location:
issm/trunk/src/c/Bamgx
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Bamgx/Mesh2.h

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

  • issm/trunk/src/c/Bamgx/objects/Triangles.cpp

    r2865 r2870  
    2020        long NbUnSwap =0;
    2121        int  Triangles::counter = 0;
    22         Triangles* CurrentTh=0; //To delete
     22        Triangles* CurrentTh=NULL; //To delete
    2323
    2424        /*Constructors/Destructors*/
Note: See TracChangeset for help on using the changeset viewer.