Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 3229) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 3230) @@ -23,6 +23,12 @@ /*CLASS: DoubleAndInt4 {{{1*/ class DoubleAndInt4 { - public: double q; Int4 i3j; - int operator<(DoubleAndInt4 a){return q > a.q;} + //class used by Triangles::MakeQuadrangles + + public: + double q; + Int4 i3j; + + //Operators + int operator<(DoubleAndInt4 a){return q > a.q;} }; /*}}}1*/ @@ -31,49 +37,58 @@ private: Icoor1 dir; - public: + + public: + //Methods Direction(): dir(MaxICoor){}; // no direction set - Direction(Icoor1 i,Icoor1 j) { Icoor2 n2 = 2*(Abs(i)+Abs(j)); - Icoor2 r = MaxICoor* (Icoor2) i; + Direction(Icoor1 i,Icoor1 j) { + Icoor2 n2 = 2*(Abs(i)+Abs(j)); + Icoor2 r = MaxICoor* (Icoor2) i; Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number - dir = (j>0) ? r1 : r1+1; // odd -> j>0 even -> j<0 - } - int sens( Icoor1 i,Icoor1 j) { int r =1; + dir = (j>0) ? r1 : r1+1; // odd-> j>0 even-> j<0 + } + int sens(Icoor1 i,Icoor1 j) { + int r =1; if (dir!= MaxICoor) { Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1); - r = (x*i + y*j) >=0;} - return r;} + r = (x*i + y*j) >=0; + } + return r; + } }; /*}}}1*/ /*CLASS: Vertex{{{1*/ class Vertex { - public: - I2 i; // allow to use i.x, and i.y in long int (beware of scale and centering) - R2 r; // allow to use r.x, and r.y in double + + public: + I2 i; // integer coordinates + R2 r; // real coordinates Metric m; Int4 ReferenceNumber; Direction DirOfSearch; + Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t union { - Triangle* t; // one triangle which contained the vertex + Triangle* t; // one triangle which is containing the vertex Int4 color; - Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed - VertexOnGeom * onGeometry; // if vint 8; // set with Triangles::SetVertexFieldOn() - Vertex * onBackgroundVertex; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() - VertexOnEdge * onBackgroundEdge; // if vint == 32 on Background edge + Vertex* to; // use in geometry Vertex to now the Mesh Vertex associed + VertexOnGeom* onGeometry; // if vint == 8; // set with Triangles::SetVertexFieldOn() + Vertex* onBackgroundVertex; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() + VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge }; - Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t - operator I2 () const {return i;} // operator de cast - operator const R2 & () const {return r;}// operator de cast - // operator R2 & () {return r;}// operator de cast - Real8 operator()(R2 x) const { return m(x);} - operator Metric () const {return m;}// operator de cast - Int4 Optim(int = 1,int =0); - // Vertex(){} - // ~Vertex(){} - Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1); - void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts); - void Echo(); - int ref() const { return ReferenceNumber;} - - inline void Set(const Vertex & rec,const Triangles &,Triangles &); + + //Operators + operator I2() const {return i;} // Cast operator + operator const R2 & () const {return r;} // Cast operator + operator Metric () const {return m;} // Cast operator + Real8 operator()(R2 x) const { return m(x);} // Get x in the metric m + + //methods (No constructor and no destructors...) + Int4 Optim(int =1,int =0); + Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1); + void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts); + void Echo(); + int ref() const { return ReferenceNumber;} + + //inline functions + inline void Set(const Vertex & rec,const Triangles &,Triangles &); }; /*}}}1*/ @@ -87,20 +102,20 @@ int a; //Edge number + //Constructors + TriangleAdjacent() {}; TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {}; - TriangleAdjacent() {}; - + + //Operators operator Triangle * () const {return t;} operator Triangle & () const {return *t;} operator int() const {return a;} - TriangleAdjacent & operator++(){ - a= NextEdge[a]; - return *this; - } - TriangleAdjacent operator--(){ - a= PreviousEdge[a]; - return *this; - } + TriangleAdjacent & operator++(){ a= NextEdge[a]; return *this; } + TriangleAdjacent operator--(){ a= PreviousEdge[a]; return *this; } + + //Methods + int swap(); + + //Inline methods inline TriangleAdjacent Adj() const ; - int swap(); inline void SetAdj2(const TriangleAdjacent& , int =0); inline Vertex * EdgeVertex(const int &) const ; @@ -118,33 +133,34 @@ /*CLASS: Edge{{{1*/ class Edge { - public: - Vertex * v[2]; - Int4 ref; - GeometricalEdge* onGeometry; - Vertex & operator[](int i){return *v[i];}; - Vertex * operator()(int i){return v[i];}; - - void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu) - { - if (v[0] >=vb && v[0] =vb && v[1] =vb && v[0] =vb && v[1] no continuity + R2 tg[2]; // the 2 tangentes (tg[0] =0 => no continuity) GeometricalEdge* Adj[2]; int DirAdj[2]; int flag ; GeometricalEdge* link; // if Cracked() or Equi() - // end of data - - GeometricalVertex & operator[](int i){return *v[i];}; + + //Operators + GeometricalVertex & operator[](int i){return *v[i];}; const GeometricalVertex & operator[](int i) const { return *v[i];}; - GeometricalVertex * operator()(int i){return v[i];}; - // inline void Set(const Geometry &,Int4,Geometry &); - + GeometricalVertex * operator()(int i){return v[i];}; + + //Methods R2 F(Real8 theta) const ; // parametrization of the curve edge Real8 R1tg(Real8 theta,R2 &t) const ; // 1/radius of curvature + tangente - int Cracked() const {return flag & 1;} - int Dup() const { return flag & 32;} - int Equi()const {return flag & 2;} - int ReverseEqui()const {return flag & 128;} - int TgA()const {return flag &4;} - int TgB()const {return flag &8;} - int Tg(int i) const{return i==0 ? TgA() : TgB();} - int Mark()const {return flag &16;} - int Required() { return flag & 64;} + int Cracked() const {return flag & 1;} + int Dup() const { return flag & 32;} + int Equi()const {return flag & 2;} + int ReverseEqui()const {return flag & 128;} + int TgA()const {return flag &4;} + int TgB()const {return flag &8;} + int Tg(int i) const{return i==0 ? TgA() : TgB();} + int Mark()const {return flag &16;} + int Required() { return flag & 64;} void SetCracked() { flag |= 1;} void SetDup() { flag |= 32;} // not a real edge @@ -216,147 +233,139 @@ void SetReverseEqui() {flag |= 128;} + //Inline methods inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); - }; /*}}}1*/ /*CLASS: Curve{{{1*/ - class Curve {public: - GeometricalEdge * be,*ee; // begin et end edge - int kb,ke; // begin vetex and end vertex - Curve *next; // next curve equi to this - bool master; // true => of equi curve point on this curve - inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew); - Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {} - void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse + class Curve { + public: + GeometricalEdge * be,*ee; // begin et end edge + int kb,ke; // begin vetex and end vertex + Curve *next; // next curve equi to this + bool master; // true => of equi curve point on this curve + + //Methods + Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {} + void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse + + //Inline methods + inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew); }; /*}}}1*/ /*CLASS: Triangle{{{1*/ class Triangle { + friend class TriangleAdjacent; - private: // les arete sont opposes a un sommet - Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer - Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles - Int1 TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1 + private: + Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer + Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles + Int1 TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1 + public: - Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) - union { - Triangle * link ; - Int4 color; - }; - void Echo(); - void SetDet() { - if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); - else det = -1; } - Triangle() {} - Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k); - Triangle(Vertex *v0,Vertex *v1,Vertex *v2); - inline void Set(const Triangle &,const Triangles &,Triangles &); - inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;} - TriangleAdjacent FindBoundaryEdge(int i) const; - - void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){ - if (link >=tb && link =tb && TriaAdjTriangles[0] =tb && TriaAdjTriangles[1] =tb && TriaAdjTriangles[2] =vb && TriaVertices[0] =vb && TriaVertices[1] =vb && TriaVertices[2] TriaAdjTriangles[aatt]=this; - tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark - } - - void SetAdj2(Int1 a,Triangle *t,Int1 aat){ - //TriaAdjTriangles = pointer toward the adjacent triangles of this - //TriaAdjSharedEdge = position of the edge in the triangle (mod 3) - TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t - TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle - //if t!=NULL add adjacent triangle to t (this) - if(t) { - t->TriaAdjTriangles[aat]=this; - t->TriaAdjSharedEdge[aat]=a; - } - } - - void SetTriangleContainingTheVertex() { - if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0); - if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1); - if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2); - } - - int swap(Int2 a1,int=0); - Int4 Optim(Int2 a,int =0); - - int Locked(int a)const { return TriaAdjSharedEdge[a]&4;} - int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} - int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;} - // for optimisation - int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;} - void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);} - - void SetHidden(int a){ - - //Get Adjacent Triangle number a - register Triangle* t = TriaAdjTriangles[a]; - - //if it exist - //C|=D -> C=(C|D) bitwise inclusive OR - if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16; - TriaAdjSharedEdge[a] |= 16; - } - void SetCracked(int a){ - register Triangle * t = TriaAdjTriangles[a]; - if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32; - TriaAdjSharedEdge[a] |= 32; - } - - double QualityQuad(int a,int option=1) const; - Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ; - - void SetLocked(int a){ - //mark the edge as on Boundary - register Triangle * t = TriaAdjTriangles[a]; - t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4; - TriaAdjSharedEdge[a] |= 4; - } - - void SetMarkUnSwap(int a){ - register Triangle * t = TriaAdjTriangles[a]; - t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8; - TriaAdjSharedEdge[a] |=8 ; - } - - void SetUnMarkUnSwap(int a){ - register Triangle * t = TriaAdjTriangles[a]; - t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 - TriaAdjSharedEdge[a] &=55 ; - } + Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) + union { + Triangle * link ; + Int4 color; + }; + + //Constructors/Destructors + Triangle() {} + Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k); + Triangle(Vertex *v0,Vertex *v1,Vertex *v2); + + //Operators + const Vertex & operator[](int i) const {return *TriaVertices[i];}; + Vertex & operator[](int i) {return *TriaVertices[i];}; + const Vertex * operator()(int i) const {return TriaVertices[i];}; + Vertex * & operator()(int i) {return TriaVertices[i];}; + + //Methods + void Echo(); + int swap(Int2 a1,int=0); + Int4 Optim(Int2 a,int =0); + int Locked(int a)const { return TriaAdjSharedEdge[a]&4;} + int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} + int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;} + int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;} + void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);} + double QualityQuad(int a,int option=1) const; + Int1 NuEdgeTriangleAdj(int i) const {return TriaAdjSharedEdge[i&3]&3;} // Number of the adjacent edge in adj tria + TriangleAdjacent FindBoundaryEdge(int i) const; + TriangleAdjacent Adj(int i) const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);}; + Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];} + Triangle* Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ; + void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){ + if (link >=tb && link =tb && TriaAdjTriangles[0] =tb && TriaAdjTriangles[1] =tb && TriaAdjTriangles[2] =vb && TriaVertices[0] =vb && TriaVertices[1] =vb && TriaVertices[2] TriaAdjTriangles[aatt]=this; + tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark + } + void SetAdj2(Int1 a,Triangle *t,Int1 aat){ + TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t + TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle + if(t) { //if t!=NULL add adjacent triangle to t (this) + t->TriaAdjTriangles[aat]=this; + t->TriaAdjSharedEdge[aat]=a; + } + } + void SetTriangleContainingTheVertex() { + if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0); + if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1); + if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2); + } + void SetHidden(int a){ + //Get Adjacent Triangle number a + register Triangle* t = TriaAdjTriangles[a]; + //if it exist + //C|=D -> C=(C|D) bitwise inclusive OR + if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16; + TriaAdjSharedEdge[a] |= 16; + } + void SetCracked(int a){ + register Triangle * t = TriaAdjTriangles[a]; + if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32; + TriaAdjSharedEdge[a] |= 32; + } + + void SetLocked(int a){ + //mark the edge as on Boundary + register Triangle * t = TriaAdjTriangles[a]; + t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4; + TriaAdjSharedEdge[a] |= 4; + } + void SetMarkUnSwap(int a){ + register Triangle * t = TriaAdjTriangles[a]; + t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8; + TriaAdjSharedEdge[a] |=8 ; + } + void SetUnMarkUnSwap(int a){ + register Triangle * t = TriaAdjTriangles[a]; + t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 + TriaAdjSharedEdge[a] &=55 ; + } + void SetDet() { + if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); + else det = -1; } + + //Inline methods + inline Real4 qualite() ; + inline void Set(const Triangle &,const Triangles &,Triangles &); + inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;} + }; // end of Triangle class @@ -382,4 +391,6 @@ Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length int last;// last index in ListofIntersectionTriangles for this Sub seg of edge + + //Methods R2 F(Real8 s){ Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01; @@ -387,36 +398,44 @@ throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd")); } - return e->F(sBegin*c0+sEnd*c1);} + return e->F(sBegin*c0+sEnd*c1); + } }; - int MaxSize; - int Size; - Real8 len; - int state; - IntersectionTriangles * lIntTria; - int NbSeg; - int MaxNbSeg; - SegInterpolation * lSegsI; - public: - IntersectionTriangles & operator[](int i) {return lIntTria[i];} - operator int&() {return Size;} - - ListofIntersectionTriangles(int n=256,int m=16) - : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) , - NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){ - long int verbosity=0; - if (verbosity>9) printf(" construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg); - }; - - ~ListofIntersectionTriangles(){ - if (lIntTria) delete [] lIntTria,lIntTria=0; - if (lSegsI) delete [] lSegsI,lSegsI=0;} - void init(){state=0;len=0;Size=0;} - - int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2); - int NewItem(R2,const Metric & ); - void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1){ + public: + + int MaxSize; + int Size; + Real8 len; + int state; + IntersectionTriangles * lIntTria; + int NbSeg; + int MaxNbSeg; + SegInterpolation * lSegsI; + + //Constructors/Destructors + ListofIntersectionTriangles(int n=256,int m=16) + : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) , + NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){ + long int verbosity=0; + if (verbosity>9) printf(" construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg); + }; + ~ListofIntersectionTriangles(){ + if (lIntTria) delete [] lIntTria,lIntTria=0; + if (lSegsI) delete [] lSegsI,lSegsI=0; + } + + //Operators + IntersectionTriangles & operator[](int i) {return lIntTria[i];} + operator int&() {return Size;} + + //Methods + void init(){state=0;len=0;Size=0;} + int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2); + int NewItem(R2,const Metric & ); + void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0); + Real8 Length(); + Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx); + void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1){ long int verbosity=0; - if (NbSeg>=MaxNbSeg) { int mneo= MaxNbSeg; @@ -439,12 +458,9 @@ lSegsI[NbSeg].sEnd=s1; NbSeg++; - } - - void ReShape() { + } + void ReShape(){ register int newsize = MaxSize*2; IntersectionTriangles* nw = new IntersectionTriangles[newsize]; - if (!nw){ - throw ErrorException(__FUNCT__,exprintf("!nw")); - } + if (!nw){ throw ErrorException(__FUNCT__,exprintf("!nw"));} // recopy for (int i=0;i=0 :0;} - VertexOnGeom(): mv(0),abscisse(0){gv=0;} - VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} - VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;} - operator Vertex * () const {return mv;} - operator GeometricalVertex * () const {return gv;} - operator GeometricalEdge * () const {return ge;} - operator const Real8 & () const {return abscisse;} - int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} - void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;} - inline void Set(const Triangles &,Int4,Triangles &); + + public: + + Vertex* mv; + Real8 abscisse; + union{ + GeometricalVertex * gv; // if abscisse <0; + GeometricalEdge * ge; // if abscisse in [0..1] + }; + + //Constructors/Destructors + VertexOnGeom(): mv(0),abscisse(0){gv=0;} + VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} + VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;} + + //Operators + operator Vertex * () const {return mv;} + operator GeometricalVertex * () const {return gv;} + operator GeometricalEdge * () const {return ge;} + operator const Real8 & () const {return abscisse;} + + //Methods + int OnGeomVertex()const {return this? abscisse <0 :0;} + int OnGeomEdge() const {return this? abscisse >=0 :0;} + int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} + void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;} + + //Inline methods + inline void Set(const Triangles &,Int4,Triangles &); + inline void Set(const VertexOnGeom&,const Triangles &,Triangles &); }; @@ -510,30 +535,48 @@ /*CLASS: VertexOnVertex{{{1*/ class VertexOnVertex { - public: - Vertex * v, *bv; + + public: + Vertex* v; + Vertex* bv; + + //Constructors VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){} VertexOnVertex() {}; + + //Methods + void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;} + + //Inline methods inline void Set(const Triangles &,Int4,Triangles &); - void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;} }; /*}}}1*/ /*CLASS: VertexOnEdge{{{1*/ class VertexOnEdge { - public: - Vertex * v; - Edge * be; - Real8 abcisse; - VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {} - VertexOnEdge(){} - inline void Set(const Triangles &,Int4,Triangles &); - void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;} - Vertex & operator[](int i) const { return (*be)[i];} - operator Real8 () const { return abcisse;} - operator Vertex * () const { return v;} - operator Edge * () const { return be;} + + public: + Vertex* v; + Edge* be; + Real8 abcisse; + + //Constructors + VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {} + VertexOnEdge(){} + + //Operators + operator Real8 () const { return abcisse;} + operator Vertex* () const { return v;} + operator Edge* () const { return be;} + Vertex & operator[](int i) const { return (*be)[i];} + + //Methods + void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;} + + //Inline methods + inline void Set(const Triangles &,Int4,Triangles &); }; /*}}}1*/ /*CLASS: CrackedEdge{{{1*/ class CrackedEdge { + friend class Triangles; @@ -541,10 +584,14 @@ friend class Triangles; friend class CrackedEdge; - Triangle * t; // edge of triangle t - int i; // edge number of in triangle - Edge *edge; // the 2 edge + Triangle* t; // edge of triangle t + int i; // edge number of in triangle + Edge *edge; // the 2 edge Vertex *New[2]; // new vertex number + + //Constructors CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;} CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;} + + //Methods void Crack(){ Triangle & T(*t); @@ -555,34 +602,39 @@ } T(i0) = New[0]; - T(i1) = New[1];} - void UnCrack(){ - Triangle & T(*t); - int i0=VerticesOfTriangularEdge[i][0]; - int i1=VerticesOfTriangularEdge[i][0]; - if (!New[0] && !New[1]){ - throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]")); - } - T(i0) = TheVertex(T(i0)); - T(i1) = TheVertex(T(i1));} - void Set() { - TriangleAdjacent ta ( FindTriangleAdjacent(*edge)); - t = ta; - i = ta; - - New[0] = ta.EdgeVertex(0); - New[1] = ta.EdgeVertex(1); - // warning the ref - } + T(i1) = New[1]; + } + void UnCrack(){ + Triangle & T(*t); + int i0=VerticesOfTriangularEdge[i][0]; + int i1=VerticesOfTriangularEdge[i][0]; + if (!New[0] && !New[1]){ + throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]")); + } + T(i0) = TheVertex(T(i0)); + T(i1) = TheVertex(T(i1)); + } + void Set() { + TriangleAdjacent ta ( FindTriangleAdjacent(*edge)); + t = ta; + i = ta; + + New[0] = ta.EdgeVertex(0); + New[1] = ta.EdgeVertex(1); + // warning the ref + } }; public: - CrackedTriangle a,b; - CrackedEdge() :a(),b() {} - CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {}; - CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {}; - - void Crack() { a.Crack(); b.Crack();} - void UnCrack() { a.UnCrack(); b.UnCrack();} - void Set() { a.Set(), b.Set();} + CrackedTriangle a,b; + + //Constructors + CrackedEdge() :a(),b() {} + CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {}; + CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {}; + + //Methods + void Crack() { a.Crack(); b.Crack();} + void UnCrack() { a.UnCrack(); b.UnCrack();} + void Set() { a.Set(), b.Set();} }; /*}}}1*/ @@ -593,76 +645,58 @@ Geometry & Gh; // Geometry Triangles & BTh; // Background Mesh Bth==*this =>no background - - Int4 NbRef; // counter of ref on the this class if 0 we can delete + Int4 NbRef; // counter of ref on the this class if 0 we can delete Int4 nbvx,nbtx; // nombre max de sommets , de triangles - - Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, - // of initial vertices, of edges with reference, - Int4 NbOfQuad; // nb of quadrangle - - Int4 NbSubDomains; // - Int4 NbOutT; // Nb of oudeside triangle + Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference, + Int4 NbOfQuad; // nb of quadrangle + Int4 NbSubDomains; + Int4 NbOutT; // Nb of oudeside triangle Int4 NbOfTriangleSearchFind; Int4 NbOfSwapTriangle; - Vertex * vertices; // data of vertices des sommets - + Vertex* vertices; Int4 NbVerticesOnGeomVertex; VertexOnGeom * VerticesOnGeomVertex; - Int4 NbVerticesOnGeomEdge; VertexOnGeom * VerticesOnGeomEdge; - Int4 NbVertexOnBThVertex; VertexOnVertex *VertexOnBThVertex; - Int4 NbVertexOnBThEdge; VertexOnEdge *VertexOnBThEdge; - Int4 NbCrackedVertices; - Int4 NbCrackedEdges; CrackedEdge *CrackedEdges; - - R2 pmin,pmax; // extrema - Real8 coefIcoor; // coef to integer Icoor1; - - Triangle * triangles; - Edge * edges; - + R2 pmin,pmax; // extrema + Real8 coefIcoor; // coef to integer Icoor1; + Triangle* triangles; + Edge* edges; QuadTree *quadtree; - Vertex ** ordre; - SubDomain * subdomains; + Vertex** ordre; + SubDomain* subdomains; ListofIntersectionTriangles lIntTria; - // end of variable - - Triangles(Int4 i);//:BTh(*this),Gh(*new Geometry()){PreInit(i);} - - ~Triangles(); - Triangles(const char * ,Real8=-1) ; + + //Constructors/Destructors Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts); Triangles(double* index,double* x,double* y,int nods,int nels); - + Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); //copy operator + Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature Triangles(Int4 nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) { try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);} catch(...) { this->~Triangles(); throw; } } - - //Genetare mesh from geometry Triangles(Int4 nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){ try { GeomToTriangles0(nbvx,bamgopts);} catch(...) { this->~Triangles(); throw; } } - - Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR - Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature - - void SetIntCoor(const char * from =0); - - Real8 MinimalHmin() {return 2.0/coefIcoor;} - Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} + ~Triangles(); + + //Operators const Vertex & operator[] (Int4 i) const { return vertices[i];}; Vertex & operator[](Int4 i) { return vertices[i];}; const Triangle & operator() (Int4 i) const { return triangles[i];}; Triangle & operator()(Int4 i) { return triangles[i];}; + + //Methods + void SetIntCoor(const char * from =0); + Real8 MinimalHmin() {return 2.0/coefIcoor;} + Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} I2 toI2(const R2 & P) const { return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) @@ -674,11 +708,9 @@ void Insert(); void ForceBoundary(); - //void Heap(); void Renumber(BamgOpts* bamgopts); void FindSubDomain(int OutSide=0); Int4 ConsRefTriangle(Int4 *) const; void ShowHistogram() const; - void ShowRegulaty() const; // Add FH avril 2007 - + void ShowRegulaty() const; void ReMakeTriangleContainingTheVertex(); void UnMarkUnSwapTriangle(); @@ -691,5 +723,5 @@ Int4 SplitInternalEdgeWithBorderVertices(); void MakeQuadrangles(double costheta); - int SplitElement(int choice); + int SplitElement(int choice); void MakeQuadTree(); void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1); @@ -699,7 +731,5 @@ void SmoothingVertex(int =3,Real8=0.3); Metric MetricAt (const R2 &) const; - GeometricalEdge * ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta, - Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR); - + GeometricalEdge* ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta, Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR); Int4 Number(const Triangle & t) const { return &t - triangles;} Int4 Number(const Triangle * t) const { return t - triangles;} @@ -709,16 +739,11 @@ Int4 Number(const Edge * t) const { return t - edges;} Int4 Number2(const Triangle * t) const { - // if(t>= triangles && t < triangles + nbt ) return t - triangles; - // else return t - OutSidesTriangles; - } - - Vertex * NearestVertex(Icoor1 i,Icoor1 j) ; - Triangle * FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; - + } + Vertex* NearestVertex(Icoor1 i,Icoor1 j) ; + Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; void ReadMesh(double* index,double* x,double* y,int nods,int nels); void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts); void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts); - void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef); void WriteMetric(BamgOpts* bamgopts); @@ -727,12 +752,11 @@ void BuildMetric1(BamgOpts* bamgopts); void IntersectGeomMetric(BamgOpts* bamgopts); - - int isCracked() const {return NbCrackedVertices != 0;} - int Crack(); - int UnCrack(); - + int isCracked() const {return NbCrackedVertices != 0;} + int Crack(); + int UnCrack(); void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); void FillHoleInMesh() ; int CrackMesh(); + private: void GeomToTriangles1(Int4 nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption @@ -744,33 +768,38 @@ /*CLASS: Geometry{{{1*/ class Geometry { - public: - Int4 NbRef; // counter of ref on the this class if 0 we can delete - - Int4 nbvx,nbtx; // nombre max de sommets , de Geometry - Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux, + + public: + Int4 NbRef; // counter of ref on the this class if 0 we can delete + Int4 nbvx,nbtx; // maximum number of vertices + Int4 nbv,nbt,nbiv,nbe; // number of vertices Int4 NbSubDomains; // - // Int4 nbtf;// de triangle frontiere Int4 NbOfCurves; - int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); } GeometricalVertex* vertices; - Triangle * triangles; - GeometricalEdge * edges; - QuadTree *quadtree; - GeometricalSubDomain *subdomains; - Curve *curves; - ~Geometry(); - Geometry(const Geometry & Gh); //Copy Operator - Geometry(int nbg,const Geometry ** ag); // intersection operator - + Triangle* triangles; + GeometricalEdge* edges; + QuadTree* quadtree; + GeometricalSubDomain* subdomains; + Curve* curves; R2 pmin,pmax; // extrema Real8 coefIcoor; // coef to integer Icoor1; Real8 MaximalAngleOfCorner; + //Constructor/Destructors + ~Geometry(); + Geometry(const Geometry & Gh); //Copy Operator + Geometry(int nbg,const Geometry** ag); // intersection operator + + //Operators + const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];}; + GeometricalVertex & operator[](Int4 i) { return vertices[i];}; + const GeometricalEdge & operator() (Int4 i) const { return edges[i];}; + GeometricalEdge & operator()(Int4 i) { return edges[i];}; + + //Methods + int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); } void Echo(); - I2 toI2(const R2 & P) const { return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} - Real8 MinimalHmin() {return 2.0/coefIcoor;} Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} @@ -780,9 +809,4 @@ void AfterRead(); Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();} - - const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];}; - GeometricalVertex & operator[](Int4 i) { return vertices[i];}; - const GeometricalEdge & operator() (Int4 i) const { return edges[i];}; - GeometricalEdge & operator()(Int4 i) { return edges[i];}; Int4 Number(const GeometricalVertex & t) const { return &t - vertices;} Int4 Number(const GeometricalVertex * t) const { return t - vertices;} @@ -790,8 +814,5 @@ Int4 Number(const GeometricalEdge * t) const { return t - edges;} Int4 Number(const Curve * c) const { return c - curves;} - - void UnMarkEdges() { - for (Int4 i=0;im,v[1]->m,v[1]->r - v[0]->r) ; - } /*}}}1*/ Index: /issm/trunk/src/c/Bamgx/meshtype.h =================================================================== --- /issm/trunk/src/c/Bamgx/meshtype.h (revision 3229) +++ /issm/trunk/src/c/Bamgx/meshtype.h (revision 3230) @@ -38,5 +38,5 @@ static const Int2 PreviousVertex[3] = {2,0,1}; #if LONG_BIT > 63 - const Icoor1 MaxICoor = 1073741823; // 2^30-1 + const Icoor1 MaxICoor = 1073741823; // 2^30-1 =111...111 (29 times) #else const Icoor1 MaxICoor = 8388608; // 2^23 Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3229) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3230) @@ -141,29 +141,4 @@ } /*}}}1*/ - /*FUNCTION Triangles::~Triangles(){{{1*/ - Triangles::~Triangles() { - long int verbosity=2; - //if(vertices) delete [] vertices; //TEST crash if not commented - if(edges) delete [] edges; - if(triangles) delete [] triangles; - if(quadtree) delete quadtree; - //if(ordre) delete [] ordre; //TEST crash if not commented - if( subdomains) delete [] subdomains; - if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge; - if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex; - if (VertexOnBThVertex) delete [] VertexOnBThVertex; - if (VertexOnBThEdge) delete [] VertexOnBThEdge; - - if (&Gh) { - if (Gh.NbRef>0) Gh.NbRef--; - else if (Gh.NbRef==0) delete &Gh; - } - if (&BTh && (&BTh != this)) { - if (BTh.NbRef>0) BTh.NbRef--; - else if (BTh.NbRef==0) delete &BTh; - } - PreInit(0); // set all to zero - } - /*}}}1*/ /*FUNCTION Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) COPY{{{1*/ Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) @@ -231,4 +206,29 @@ } + /*}}}1*/ + /*FUNCTION Triangles::~Triangles(){{{1*/ + Triangles::~Triangles() { + long int verbosity=2; + //if(vertices) delete [] vertices; //TEST crash if not commented + if(edges) delete [] edges; + if(triangles) delete [] triangles; + if(quadtree) delete quadtree; + //if(ordre) delete [] ordre; //TEST crash if not commented + if( subdomains) delete [] subdomains; + if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge; + if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex; + if (VertexOnBThVertex) delete [] VertexOnBThVertex; + if (VertexOnBThEdge) delete [] VertexOnBThEdge; + + if (&Gh) { + if (Gh.NbRef>0) Gh.NbRef--; + else if (Gh.NbRef==0) delete &Gh; + } + if (&BTh && (&BTh != this)) { + if (BTh.NbRef>0) BTh.NbRef--; + else if (BTh.NbRef==0) delete &BTh; + } + PreInit(0); // set all to zero + } /*}}}1*/ @@ -2190,11 +2190,16 @@ TriangleAdjacent ta(0,0); for (Int4 i = 0; i < nbe; i++){ + + //Force edge i nbswp = ForceEdge(edges[i][0],edges[i][1],ta); - - if ( nbswp < 0) k++; + if (nbswp<0) k++; else Nbswap += nbswp; + if (nbswp) nbfe++; if ( nbswp < 0 && k < 5){ - throw ErrorException(__FUNCT__,exprintf("Missing Edge %i, v0=%i,v1=%i",i ,Number(edges[i][0]),Number(edges[i][1]))); + for (Int4 j = 0; j < nbe; j++){ + printf("Edge %i: %i %i\n",j,Number(edges[j][0]),Number(edges[j][1])); + } + throw ErrorException(__FUNCT__,exprintf("Missing Edge %i, v0=%i,v1=%i",i,Number(edges[i][0]),Number(edges[i][1]))); } if ( nbswp >=0 && edges[i].onGeometry->Cracked()) @@ -3154,5 +3159,4 @@ FindSubDomain(); - // NewPointsOld(*this) ; if (verbose>3) printf(" Inserting internal points\n"); NewPoints(*this,bamgopts,0) ; Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 3229) +++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 3230) @@ -20,4 +20,5 @@ /*FUNCTION Vertex::Smoothing{{{1*/ Real8 Vertex::Smoothing(Triangles &Th,const Triangles &BTh,Triangle* &tstart ,Real8 omega){ + register Vertex* s=this; Vertex &vP = *s,vPsave=vP; @@ -43,5 +44,5 @@ j = NextEdge[jc]; if (k>=2000){ - throw ErrorException(__FUNCT__,exprintf("k>=2000")); + throw ErrorException(__FUNCT__,exprintf("k>=2000 (Maximum number of iterations reached)")); } } while ( tbegin != tria);