Changeset 2865


Ignore:
Timestamp:
01/20/10 08:08:23 (15 years ago)
Author:
Mathieu Morlighem
Message:

reverted back to 2863 (cleaner version)

Location:
issm/trunk/src/c
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

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

    r2864 r2865  
    1717#include <stdio.h>
    1818#include <string.h>
     19#include "Mesh2.h"
    1920#include "QuadTree.h"
    2021

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

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

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

    r2864 r2865  
    44#include <time.h>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

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

    r2864 r2865  
    44#include <ctime>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

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

    r2864 r2865  
    44#include <ctime>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

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

    r2864 r2865  
    44#include <ctime>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

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

    r2864 r2865  
    44#include <time.h>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

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

    r2864 r2865  
    33#include <stdlib.h>
    44
    5 #include "../Bamgx.h"
     5#include "../Mesh2.h"
    66#include "../QuadTree.h"
    77

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

    r2864 r2865  
    44#include <ctime>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

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

    r2864 r2865  
    44#include <ctime>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

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

    r2864 r2865  
    44#include <ctime>
    55
    6 #include "../Bamgx.h"
     6#include "../Mesh2.h"
    77#include "../QuadTree.h"
    88#include "../SetOfE4.h"

  • issm/trunk/src/c/Makefile.am

    r2864 r2865  
    320320                                        ./Bamgx/Bamgx.h\
    321321                                        ./Bamgx/Bamgx.cpp\
     322                                        ./Bamgx/Mesh2.h \
    322323                                        ./Bamgx/meshtype.h \
    323324                                        ./Bamgx/objects/MatVVP2x2.cpp \
    … …  
    660661                                        ./Bamgx/Bamgx.h\
    661662                                        ./Bamgx/Bamgx.cpp\
     663                                        ./Bamgx/Mesh2.h \
    662664                                        ./Bamgx/meshtype.h \
    663665                                        ./Bamgx/objects/MatVVP2x2.cpp \
Note: See TracChangeset for help on using the changeset viewer.