Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2865) @@ -17,4 +17,5 @@ #include #include +#include "Mesh2.h" #include "QuadTree.h" Index: /issm/trunk/src/c/Bamgx/Bamgx.h =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.h (revision 2864) +++ /issm/trunk/src/c/Bamgx/Bamgx.h (revision 2865) @@ -7,1312 +7,7 @@ #include "../objects/objects.h" -#include "../shared/shared.h" -#include "../include/macros.h" -#include "../toolkits/toolkits.h" - -//From MeshIo -#include -#include -#include -#include -#include -#include -#include -#include - -#if (defined(unix) || defined(__unix)) && !defined(__AIX) -#define SYSTIMES -#include -#include -#endif - -#include "meshtype.h" -#include "R2.h" /* local prototypes: */ int Bamgx(BamgMesh* bamgmesh,BamgGeom* bamggeom,BamgOpts* bamgopts); -/*Others*/ -using namespace std; -namespace bamg { - - const double Pi = 3.14159265358979323846264338328; - const float fPi = 3.14159265358979323846264338328; - - extern int hinterpole; - - typedef P2 I2; - - inline int BinaryRand(){ -#ifdef RAND_MAX - const long HalfRandMax = RAND_MAX/2; - return rand() < HalfRandMax; -#else - return rand() & 16384; // 2^14 (for sun because RAND_MAX is not def in stdlib.h) -#endif - - } - typedef P2 R2; - typedef P2xP2 I2xI2; - typedef P2 R2xR2; -} - -#include "Metric.h" - -namespace bamg { - inline float OppositeAngle(float a) - {return a<0 ? fPi + a :a - fPi ;} - inline double OppositeAngle(double a) - {return a<0 ? Pi + a :a - Pi ;} - - Icoor2 inline det(const I2 &a,const I2 & b,const I2 &c) - { - register Icoor2 bax = b.x - a.x ,bay = b.y - a.y; - register Icoor2 cax = c.x - a.x ,cay = c.y - a.y; - return bax*cay - bay*cax;} - - - // def de numerotation dans un triangles - static const Int2 VerticesOfTriangularEdge[3][2] = {{1,2},{2,0},{0,1}}; - static const Int2 EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}}; - static const Int2 OppositeVertex[3] = {0,1,2}; - static const Int2 OppositeEdge[3] = {0,1,2}; - static const Int2 NextEdge[3] = {1,2,0}; - static const Int2 PreviousEdge[3] = {2,0,1}; - static const Int2 NextVertex[3] = {1,2,0}; - static const Int2 PreviousVertex[3] = {2,0,1}; - - Int4 AGoodNumberPrimeWith(Int4 n); - - // remark all the angle are in radian beetwen [-Pi,Pi] - - - class Geometry; - //static Geometry *NULLGeometry=0; - class Triangles; - class Triangle; - class QuadTree; - class GeometricalEdge; - class VertexOnGeom; - class VertexOnEdge; - ///////////////////////////////////////////////////////////////////////////////////// - const int IsVertexOnGeom = 8; - const int IsVertexOnVertex = 16; - const int IsVertexOnEdge = 32; - ///////////////////////////////////////////////////////////////////////////////////// - - //class from MeshGeom - class DoubleAndInt4 { - public: double q; Int4 i3j; - int operator<(DoubleAndInt4 a){return q > a.q;} - };// to sort by decroissant - template inline void HeapSort(T *c,long n){ - long l,j,r,i; - T crit; - c--; // on decale de 1 pour que le tableau commence a 1 - if( n <= 1) return; - l = n/2 + 1; - r = n; - while (1) { // label 2 - if(l <= 1 ) { // label 20 - crit = c[r]; - c[r--] = c[1]; - if ( r == 1 ) { c[1]=crit; return;} - } else crit = c[--l]; - j=l; - while (1) {// label 4 - i=j; - j=2*j; - if (j>r) {c[i]=crit;break;} // L8 -> G2 - if ((j G2 - } - } - } - //End from MeshGeom - - - class Direction { // - private: - Icoor1 dir; - public: - Direction(): dir(MaxICoor){}; // no direction set - 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; - 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;} - }; - ///////////////////////////////////////////////////////////////////////////////////// - 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 - Metric m; - Int4 ReferenceNumber; - Direction DirOfSearch; - union { - Triangle * t; // one triangle which contained the vertex - Int4 color; - Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed - VertexOnGeom * on; // if vint 8; // set with Triangles::SetVertexFieldOn() - Vertex * onbv; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() - VertexOnEdge * onbe; // 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); - int ref() const { return ReferenceNumber;} - - inline void Set(const Vertex & rec,const Triangles &,Triangles &); - }; - - double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &); - - class TriangleAdjacent { - - public: - Triangle * t; // le triangle - int a; // le numero de l arete - - TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {}; - TriangleAdjacent() {}; - - 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;} - inline TriangleAdjacent Adj() const ; - int swap(); - inline void SetAdj2(const TriangleAdjacent& , int =0); - inline Vertex * EdgeVertex(const int &) const ; - inline Vertex * OppositeVertex() const ; - inline Icoor2 & det() const; - inline int Locked() const ; - inline int GetAllFlag_UnSwap() const ; - inline void SetLock(); - inline int MarkUnSwap() const; - inline void SetMarkUnSwap(); - inline void SetCracked(); - inline int Cracked() const ; - };// end of Vertex class - - class Edge { public: - Vertex * v[2]; - Int4 ref; - GeometricalEdge * on; - 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] no continuite - GeometricalEdge * Adj [2]; - int SensAdj[2]; - // private: - int flag ; - public: - GeometricalEdge * link; // if Cracked() or Equi() - - // end of data - - 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 &); - - 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;} - void SetCracked() { flag |= 1;} - void SetDup() { flag |= 32;} // not a real edge - void SetEqui() { flag |= 2;} - void SetTgA() { flag|=4;} - void SetTgB() { flag|=8;} - void SetMark() { flag|=16;} - void SetUnMark() { flag &= 1007 /* 1023-16*/;} - void SetRequired() { flag |= 64;} - void SetReverseEqui() {flag |= 128;} - - inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); - - }; - - 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 Triangle { - friend class TriangleAdjacent; - - - private: // les arete sont opposes a un sommet - Vertex * ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer - Triangle * at [3]; // nu triangle adjacent - Int1 aa[3]; // les nu des arete dans le triangles (mod 4) - public: - Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) - union { - Triangle * link ; - Int4 color; - }; - void SetDet() { - if(ns[0] && ns[1] && ns[2]) det = bamg::det(*ns[0],*ns[1],*ns[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 ns[0]==v || ns[1]==v || ns[2]==v ;} - TriangleAdjacent FindBoundaryEdge(int ) const; - - void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu) - { - if (link >=tb && link =tb && at[0] =tb && at[1] =tb && at[2] =vb && ns[0] =vb && ns[1] =vb && ns[2] at[aatt]=this; - tt->aa[aatt]=a + (aa[a] & 60 ) ;}// Copy all the mark - } - - void SetAdj2(Int1 a,Triangle *t,Int1 aat) - { at[a]=t;aa[a]=aat; - if(t) {t->at[aat]=this;t->aa[aat]=a;} - } - - void SetTriangleContainingTheVertex() - { - if (ns[0]) (ns[0]->t=this,ns[0]->vint=0); - if (ns[1]) (ns[1]->t=this,ns[1]->vint=1); - if (ns[2]) (ns[2]->t=this,ns[2]->vint=2); - } - - int swap(Int2 a1,int=0); - Int4 Optim(Int2 a,int =0); - - int Locked(int a)const { return aa[a]&4;} - int Hidden(int a)const { return aa[a]&16;} - int Cracked(int a) const { return aa[a] & 32;} - // for optimisation - int GetAllflag(int a){return aa[a] & 1020;} - void SetAllFlag(int a,int f){aa[a] = (aa[a] &3) + (1020 & f);} - - void SetHidden(int a){ - register Triangle * t = at[a]; - if(t) t->aa[aa[a] & 3] |=16; - aa[a] |= 16;} - void SetCracked(int a){ - register Triangle * t = at[a]; - if(t) t->aa[aa[a] & 3] |=32; - aa[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){ - register Triangle * t = at[a]; - t->aa[aa[a] & 3] |=4; - aa[a] |= 4;} - - void SetMarkUnSwap(int a){ - register Triangle * t = at[a]; - t->aa[aa[a] & 3] |=8; - aa[a] |=8 ;} - - - void SetUnMarkUnSwap(int a){ - register Triangle * t = at[a]; - t->aa[aa[a] & 3] &=55; // 23 + 32 - aa[a] &=55 ;} - - }; // end of Triangle class - - class ListofIntersectionTriangles { - ///////////////////////////////////////////////////////////////////////////////////// - class IntersectionTriangles { - public: - Triangle *t; - Real8 bary[3]; // use if t != 0 - R2 x; - Metric m; - Real8 s;// abscisse curviline - Real8 sp; // len of the previous seg in m - Real8 sn;// len of the next seg in m - }; - ///////////////////////////////////////////////////////////////////////////////////// - class SegInterpolation { - public: - GeometricalEdge * e; - Real8 sBegin,sEnd; // abscisse of the seg on edge parameter - Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length - int last;// last index in ListofIntersectionTriangles for this Sub seg of edge - R2 F(Real8 s){ - Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01; - if (lBegin>s || s>lEnd){ - throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd")); - } - 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) - { - long int verbosity=0; - - if (NbSeg>=MaxNbSeg) { - int mneo= MaxNbSeg; - MaxNbSeg *= 2; - if (verbosity>3){ - printf(" reshape lSegsI from %i to %i\n",mneo,MaxNbSeg); - } - SegInterpolation * lEn = new SegInterpolation[MaxNbSeg]; - if (!lSegsI || NbSeg>=MaxNbSeg){ - throw ErrorException(__FUNCT__,exprintf("!lSegsI || NbSeg>=MaxNbSeg")); - } - for (int i=0;i< NbSeg;i++) - lEn[i] = lSegsI[MaxNbSeg]; // copy old to new - delete [] lSegsI; // remove old - lSegsI = lEn; - } - if (NbSeg) - lSegsI[NbSeg-1].last=Size; - lSegsI[NbSeg].e=e; - lSegsI[NbSeg].sBegin=s0; - lSegsI[NbSeg].sEnd=s1; - NbSeg++; - } - - // void CopyMetric(int i,int j){ lIntTria[j].m=lIntTria[i].m;} - // void CopyMetric(const Metric & mm,int j){ lIntTria[j].m=mm;} - - void ReShape() { - register int newsize = MaxSize*2; - IntersectionTriangles* nw = new IntersectionTriangles[newsize]; - if (!nw){ - throw ErrorException(__FUNCT__,exprintf("!nw")); - } - for (int i=0;i3) printf(" ListofIntersectionTriangles ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg); - MaxSize = newsize; - delete [] lIntTria;// remove old - lIntTria = nw; // copy pointer - } - - void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0); - Real8 Length(); - Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx); - }; - - - ///////////////////////////////////////////////////////////////////////////////////// - class GeometricalSubDomain { - public: - GeometricalEdge *edge; - int sens; // -1 or 1 - Int4 ref; - inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &); - - }; - ///////////////////////////////////////////////////////////////////////////////////// - class SubDomain { - public: - Triangle * head; - Int4 ref; - int sens; // -1 or 1 - Edge * edge; // to geometrical - inline void Set(const Triangles &,Int4,Triangles &); - - }; - ///////////////////////////////////////////////////////////////////////////////////// - class VertexOnGeom { public: - - Vertex * mv; - Real8 abscisse; - union{ - GeometricalVertex * gv; // if abscisse <0; - GeometricalEdge * ge; // if abscisse in [0..1] - }; - inline void Set(const VertexOnGeom&,const Triangles &,Triangles &); - int OnGeomVertex()const {return this? abscisse <0 :0;} - int OnGeomEdge() const {return this? abscisse >=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->on=this;mv->vint=IsVertexOnGeom;} - inline void Set(const Triangles &,Int4,Triangles &); - - }; - ///////////////////////////////////////////////////////////////////////////////////// - class VertexOnVertex {public: - Vertex * v, *bv; - VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){} - VertexOnVertex() {}; - inline void Set(const Triangles &,Int4,Triangles &); - void SetOnBTh(){v->onbv=bv;v->vint=IsVertexOnVertex;} - }; - ///////////////////////////////////////////////////////////////////////////////////// - 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->onbe=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;} - }; - - inline TriangleAdjacent FindTriangleAdjacent(Edge &E); - inline Vertex * TheVertex(Vertex * a); // for remove crak in mesh - ///////////////////////////////////////////////////////////////////////////////////// - - class CrackedEdge { // a small class to store on crack an uncrack information - friend class Triangles; - class CrackedTriangle { - friend class Triangles; - friend class CrackedEdge; - Triangle * t; // edge of triangle t - int i; // edge number of in triangle - Edge *edge; // the 2 edge - Vertex *New[2]; // new vertex number - 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;} - void Crack(){ - 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) = 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 - - } - - }; // end of class CrackedTriangle - 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();} - }; - - ///////////////////////////////////////////////////////////////////////////////////// - class Triangles { - public: - - enum TypeFileMesh { - AutoMesh=0,BDMesh=1,NOPOMesh=2,amMesh=3,am_fmtMesh=4,amdbaMesh=5, - ftqMesh=6,mshMesh=7}; - - int static counter; // to kown the number of mesh in memory - int OnDisk; // true if on disk - 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 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 NbOfTriangleSearchFind; - Int4 NbOfSwapTriangle; - char * name, *identity; - Vertex * vertices; // data of vertices des sommets - - 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; - - QuadTree *quadtree; - 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) ; - Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts); - - Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) - :Gh(BT.Gh),BTh(BT) { - try {GeomToTriangles1(nbvx,keepBackVertices);} - catch(...) { this->~Triangles(); throw; } } - - Triangles(Int4 nbvx,Geometry & G) - :Gh(G),BTh(*this){ - try { GeomToTriangles0(nbvx);} - 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); // truncature - - void SetIntCoor(const char * from =0); - - // void RandomInit(); - // void CubeInit(int ,int); - - Real8 MinimalHmin() {return 2.0/coefIcoor;} - Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} - 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];}; - I2 toI2(const R2 & P) const { - return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) - ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} - R2 toR2(const I2 & P) const { - return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);} - void Add( Vertex & s,Triangle * t,Icoor2 * =0) ; - void Insert(); - // void InsertOld(); - void ForceBoundary(); - void Heap(); - void FindSubDomain(int ); - Int4 ConsRefTriangle(Int4 *) const; - void ShowHistogram() const; - void ShowRegulaty() const; // Add FH avril 2007 - // void ConsLinkTriangle(); - - void ReMakeTriangleContainingTheVertex(); - void UnMarkUnSwapTriangle(); - void SmoothMetric(Real8 raisonmax) ; - void BoundAnisotropy(Real8 anisomax,double hminaniso= 1e-100) ; - void MaxSubDivision(Real8 maxsubdiv); - Edge** MakeGeometricalEdgeToEdge(); - void SetVertexFieldOn(); - void SetVertexFieldOnBTh(); - Int4 SplitInternalEdgeWithBorderVertices(); - void MakeQuadrangles(double costheta); - int SplitElement(int choice); - void MakeQuadTree(); - void NewPoints( Triangles &,int KeepBackVertex =1 ); - Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; - void NewPoints(int KeepBackVertex=1){ NewPoints(*this,KeepBackVertex);} - void ReNumberingTheTriangleBySubDomain(bool justcompress=false); - void ReNumberingVertex(Int4 * renu); - 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); - - Int4 Number(const Triangle & t) const { return &t - triangles;} - Int4 Number(const Triangle * t) const { return t - triangles;} - Int4 Number(const Vertex & t) const { return &t - vertices;} - Int4 Number(const Vertex * t) const { return t - vertices;} - Int4 Number(const Edge & t) const { return &t - edges;} - 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 * FindTriangleContening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; - - 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); - void IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols, - const Real8 hmin,const Real8 hmax, const Real8 coef, - const Real8 anisomax,const Real8 CutOff=1.e-4,const int NbJacobi=1, - const int DoNormalisation=1, - const double power=1.0, - const int choise=0); - void IntersectGeomMetric(const Real8 err,const int iso); - - - int isCracked() const {return NbCrackedVertices != 0;} - int Crack(); - int UnCrack(); - - void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo - void FillHoleInMesh() ; - int CrackMesh(); - private: - void GeomToTriangles1(Int4 nbvx,int KeepBackVertices=1);// the real constructor mesh adaption - void GeomToTriangles0(Int4 nbvx);// the real constructor mesh generator - void PreInit(Int4,char * =0 ); - - }; // End Class Triangles - - - ///////////////////////////////////////////////////////////////////////////////////// - class Geometry { - public: - int OnDisk; - Int4 NbRef; // counter of ref on the this class if 0 we can delete - - char *name; - Int4 nbvx,nbtx; // nombre max de sommets , de Geometry - Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux, - Int4 NbSubDomains; // - Int4 NbEquiEdges; - Int4 NbCrackedEdges; - // Int4 nbtf;// de triangle frontiere - Int4 NbOfCurves; - int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); } - GeometricalVertex * vertices; // data of vertices des sommets - 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 - - R2 pmin,pmax; // extrema - Real8 coefIcoor; // coef to integer Icoor1; - Real8 MaximalAngleOfCorner; - - // end of data - - - 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);} - void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); - void EmptyGeometry(); - Geometry() {EmptyGeometry();}// empty Geometry - void AfterRead(); - Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();OnDisk=1;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;} - Int4 Number(const GeometricalEdge & t) const { return &t - edges;} - Int4 Number(const GeometricalEdge * t) const { return t - edges;} - Int4 Number(const Curve * c) const { return c - curves;} - - void UnMarkEdges() { - for (Int4 i=0;i=0) { - t = at[a]; - // if (t-this<0) return 0; - v2 = ns[VerticesOfTriangularEdge[a][0]]; - v0 = ns[VerticesOfTriangularEdge[a][1]]; - v1 = ns[OppositeEdge[a]]; - v3 = t->ns[OppositeEdge[aa[a]&3]]; - } - } - return t; - } - - inline double Triangle::QualityQuad(int a,int option) const - { // first do the logique part - double q; - if (!link || aa[a] &4) - q= -1; - else { - Triangle * t = at[a]; - if (t-this<0) q= -1;// because we do 2 times - else if (!t->link ) q= -1; - else if (aa[0] & 16 || aa[1] & 16 || aa[2] & 16 || t->aa[0] & 16 || t->aa[1] & 16 || t->aa[2] & 16 ) - q= -1; - else if(option) - { - const Vertex & v2 = *ns[VerticesOfTriangularEdge[a][0]]; - const Vertex & v0 = *ns[VerticesOfTriangularEdge[a][1]]; - const Vertex & v1 = *ns[OppositeEdge[a]]; - const Vertex & v3 = * t->ns[OppositeEdge[aa[a]&3]]; - q = QuadQuality(v0,v1,v2,v3); // do the float part - } - else q= 1; - } - return q; - } - - - inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ) - { - *this = rec; - } - inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ) - { - *this = rec; - } - inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) - { - *this = Th.edges[i]; - v[0] = ThNew.vertices + Th.Number(v[0]); - v[1] = ThNew.vertices + Th.Number(v[1]); - if (on) - on = ThNew.Gh.edges+Th.Gh.Number(on); - if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]); - if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]); - - } - inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew) - { - *this = rec; - v[0] = GhNew.vertices + Gh.Number(v[0]); - v[1] = GhNew.vertices + Gh.Number(v[1]); - if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]); - if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]); - } - - inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew) - { - *this = rec; - be = GhNew.edges + Gh.Number(be); - ee = GhNew.edges + Gh.Number(ee); - if(next) next= GhNew.curves + Gh.Number(next); - } - - inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew) - { - *this = rec; - if ( ns[0] ) ns[0] = ThNew.vertices + Th.Number(ns[0]); - if ( ns[1] ) ns[1] = ThNew.vertices + Th.Number(ns[1]); - if ( ns[2] ) ns[2] = ThNew.vertices + Th.Number(ns[2]); - if(at[0]) at[0] = ThNew.triangles + Th.Number(at[0]); - if(at[1]) at[1] = ThNew.triangles + Th.Number(at[1]); - if(at[2]) at[2] = ThNew.triangles + Th.Number(at[2]); - if (link >= Th.triangles && link < Th.triangles + Th.nbt) - link = ThNew.triangles + Th.Number(link); - } - inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) - { - *this = Th.VertexOnBThVertex[i]; - v = ThNew.vertices + Th.Number(v); - - } - inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) - { - *this = Th.subdomains[i]; - if ( head-Th.triangles<0 || head-Th.triangles>=Th.nbt){ - throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt")); - } - head = ThNew.triangles + Th.Number(head) ; - if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{ - throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe")); - } - edge = ThNew.edges+ Th.Number(edge); - } - inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew) - { - *this = rec; - edge = Gh.Number(edge) + GhNew.edges; - } - inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) - { - *this = Th.VertexOnBThEdge[i]; - v = ThNew.vertices + Th.Number(v); - } - - inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew) - { - *this = rec; - mv = ThNew.vertices + Th.Number(mv); - if (gv) - if (abscisse < 0 ) - gv = ThNew.Gh.vertices + Th.Gh.Number(gv); - else - ge = ThNew.Gh.edges + Th.Gh.Number(ge); - - } - inline Real8 Edge::MetricLength() const - { - return LengthInterpole(v[0]->m,v[1]->m,v[1]->r - v[0]->r) ; - } - - inline void Triangles::ReMakeTriangleContainingTheVertex() - { - register Int4 i; - for ( i=0;iat[a]=ta.t; - t->aa[a]=ta.a|l;} - if(ta.t) { - ta.t->at[ta.a] = t ; - ta.t->aa[ta.a] = a| l ; - } - } - - - inline int TriangleAdjacent::Locked() const - { return t->aa[a] &4;} - inline int TriangleAdjacent::Cracked() const - { return t->aa[a] &32;} - inline int TriangleAdjacent::GetAllFlag_UnSwap() const - { return t->aa[a] & 1012;} // take all flag except MarkUnSwap - - inline int TriangleAdjacent::MarkUnSwap() const - { return t->aa[a] &8;} - - inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);} - - inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);} - - inline TriangleAdjacent TriangleAdjacent::Adj() const - { return t->Adj(a);} - - inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const - {return t->ns[VerticesOfTriangularEdge[a][i]]; } - inline Vertex * TriangleAdjacent::OppositeVertex() const - {return t->ns[bamg::OppositeVertex[a]]; } - inline Icoor2 & TriangleAdjacent::det() const - { return t->det;} - inline TriangleAdjacent Adj(const TriangleAdjacent & a) - { return a.Adj();} - - inline TriangleAdjacent Next(const TriangleAdjacent & ta) - { return TriangleAdjacent(ta.t,NextEdge[ta.a]);} - - inline TriangleAdjacent Previous(const TriangleAdjacent & ta) - { return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);} - - inline void Adj(GeometricalEdge * & on,int &i) - {int j=i;i=on->SensAdj[i];on=on->Adj[j];} - - inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc) - { - Real4 ret; - I2 ia=va,ib=vb,ic=vc; - I2 ab=ib-ia,bc=ic-ib,ac=ic-ia; - Icoor2 deta=Det(ab,ac); - if (deta <=0) ret = -1; - else { - Real8 a = sqrt((Real8) (ac,ac)), - b = sqrt((Real8) (bc,bc)), - c = sqrt((Real8) (ab,ab)), - p = a+b+c; - Real8 h= Max(Max(a,b),c),ro=deta/p; - ret = ro/h;} - return ret; - } - - - inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) { - Vertex *v=Th->vertices; - Int4 nbv = Th->nbv; - if (i<0 || j<0 || k<0){ - throw ErrorException(__FUNCT__,exprintf("i<0 || j<0 || k<0")); - } - if (i>=nbv || j>=nbv || k>=nbv){ - throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv")); - } - ns[0]=v+i; - ns[1]=v+j; - ns[2]=v+k; - at[0]=at[1]=at[2]=0; - aa[0]=aa[1]=aa[2]=0; - det=0; - } - - inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){ - ns[0]=v0; - ns[1]=v1; - ns[2]=v2; - at[0]=at[1]=at[2]=0; - aa[0]=aa[1]=aa[2]=0; - if (v0) det=0; - else { - det=-1; - link=NULL;}; - } - - inline Real4 Triangle::qualite() - { - return det < 0 ? -1 : bamg::qualite(*ns[0],*ns[1],*ns[2]); - } - - Int4 inline Vertex::Optim(int i,int koption) - { - Int4 ret=0; - if ( t && (vint >= 0 ) && (vint <3) ) - { - ret = t->Optim(vint,koption); - if(!i) - { - t =0; // for no future optime - vint= 0; } - } - return ret; - } - - Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c) - { - register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y; - register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y; - return bax*cay - bay*cax;} - - - void swap(Triangle *t1,Int1 a1, - Triangle *t2,Int1 a2, - Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2); - - - - int inline TriangleAdjacent::swap() - { return t->swap(a);} - - - - int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , - TriangleAdjacent & tt1,Icoor2 & dets1, - Icoor2 & detsa,Icoor2 & detsb, int & nbswap); - - int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ; - - // inline bofbof FH - inline TriangleAdjacent FindTriangleAdjacent(Edge &E) - { - Vertex * a = E.v[0]; - Vertex * b = E.v[1]; - - Triangle * t = a->t; - int i = a->vint; - TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge - if (!t || i<0 || i>=3){ - throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); - } - if ( a!=(*t)(i)){ - throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); - } - int k=0; - do { // turn around vertex in direct sens (trigo) - k++; - if (k>=20000){ - throw ErrorException(__FUNCT__,exprintf("k>=20000")); - } - // in no crack => ta.EdgeVertex(1) == a otherwise ??? - if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find - ta = ta.Adj(); - if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find - --ta; - } while (t != (Triangle *)ta); - throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found")); - return TriangleAdjacent(0,0);//for compiler - } - - inline Vertex * TheVertex(Vertex * a) // give a unique vertex with smallest number - { // in case on crack in mesh - Vertex * r(a), *rr; - Triangle * t = a->t; - int i = a->vint; - TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge - if (!t || i<0 || i>=3){ - throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); - } - if ( a!=(*t)(i)){ - throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); - } - int k=0; - do { // turn around vertex in direct sens (trigo) - k++; - if (k>=20000){ - throw ErrorException(__FUNCT__,exprintf("k>=20000")); - } - // in no crack => ta.EdgeVertex(1) == a - if ((rr=ta.EdgeVertex(0)) < r) r = rr; - ta = ta.Adj(); - if ((rr=ta.EdgeVertex(1)) < r) r =rr; - --ta; - } while (t != (Triangle*) ta); - return r; - } - - inline double CPUtime(){ -#ifdef SYSTIMES - struct tms buf; - if (times(&buf)!=-1) - return ((double)buf.tms_utime+(double)buf.tms_stime)/(long) sysconf(_SC_CLK_TCK); - else -#endif - return ((double) clock())/CLOCKS_PER_SEC; - } - -} #endif /* _BAMGX_H */ Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2865) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2865) @@ -0,0 +1,1323 @@ +#ifndef _MESH2_H_ +#define _MESH2_H_ + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" + +//From MeshIo +#include +#include +#include +#include +#include +#include +#include +#include + +#if (defined(unix) || defined(__unix)) && !defined(__AIX) +#define SYSTIMES +#include +#include +#endif + +#include "meshtype.h" +#include "R2.h" + +using namespace std; +namespace bamg { + +const double Pi = 3.14159265358979323846264338328; +const float fPi = 3.14159265358979323846264338328; + +extern int hinterpole; + +typedef P2 I2; + +inline int BinaryRand(){ +#ifdef RAND_MAX + const long HalfRandMax = RAND_MAX/2; + return rand() < HalfRandMax; +#else + return rand() & 16384; // 2^14 (for sun because RAND_MAX is not def in stdlib.h) +#endif + +} +typedef P2 R2; +typedef P2xP2 I2xI2; +typedef P2 R2xR2; +} + +#include "Metric.h" + +namespace bamg { +inline float OppositeAngle(float a) + {return a<0 ? fPi + a :a - fPi ;} +inline double OppositeAngle(double a) + {return a<0 ? Pi + a :a - Pi ;} + +Icoor2 inline det(const I2 &a,const I2 & b,const I2 &c) +{ + register Icoor2 bax = b.x - a.x ,bay = b.y - a.y; + register Icoor2 cax = c.x - a.x ,cay = c.y - a.y; + return bax*cay - bay*cax;} + + +// def de numerotation dans un triangles +static const Int2 VerticesOfTriangularEdge[3][2] = {{1,2},{2,0},{0,1}}; +static const Int2 EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}}; +static const Int2 OppositeVertex[3] = {0,1,2}; +static const Int2 OppositeEdge[3] = {0,1,2}; +static const Int2 NextEdge[3] = {1,2,0}; +static const Int2 PreviousEdge[3] = {2,0,1}; +static const Int2 NextVertex[3] = {1,2,0}; +static const Int2 PreviousVertex[3] = {2,0,1}; + +Int4 AGoodNumberPrimeWith(Int4 n); + +// remark all the angle are in radian beetwen [-Pi,Pi] + + +class Geometry; +//static Geometry *NULLGeometry=0; +class Triangles; +class Triangle; +class QuadTree; +class GeometricalEdge; +class VertexOnGeom; +class VertexOnEdge; +///////////////////////////////////////////////////////////////////////////////////// +const int IsVertexOnGeom = 8; +const int IsVertexOnVertex = 16; +const int IsVertexOnEdge = 32; +///////////////////////////////////////////////////////////////////////////////////// + +//class from MeshGeom +class DoubleAndInt4 { +public: double q; Int4 i3j; + int operator<(DoubleAndInt4 a){return q > a.q;} +};// to sort by decroissant +template inline void HeapSort(T *c,long n){ + long l,j,r,i; + T crit; + c--; // on decale de 1 pour que le tableau commence a 1 + if( n <= 1) return; + l = n/2 + 1; + r = n; + while (1) { // label 2 + if(l <= 1 ) { // label 20 + crit = c[r]; + c[r--] = c[1]; + if ( r == 1 ) { c[1]=crit; return;} + } else crit = c[--l]; + j=l; + while (1) {// label 4 + i=j; + j=2*j; + if (j>r) {c[i]=crit;break;} // L8 -> G2 + if ((j G2 + } + } +} +//End from MeshGeom + + +class Direction { // + private: + Icoor1 dir; + public: + Direction(): dir(MaxICoor){}; // no direction set + 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; + 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;} +}; +///////////////////////////////////////////////////////////////////////////////////// +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 + Metric m; + Int4 ReferenceNumber; + Direction DirOfSearch; + union { + Triangle * t; // one triangle which contained the vertex + Int4 color; + Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed + VertexOnGeom * on; // if vint 8; // set with Triangles::SetVertexFieldOn() + Vertex * onbv; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() + VertexOnEdge * onbe; // 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); + int ref() const { return ReferenceNumber;} + + inline void Set(const Vertex & rec,const Triangles &,Triangles &); +}; + +double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &); + +// extern Vertex *Meshend , *Meshbegin; + +///////////////////////////////////////////////////////////////////////////////////// +class TriangleAdjacent { + +public: + Triangle * t; // le triangle + int a; // le numero de l arete + + TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {}; + TriangleAdjacent() {}; + + 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;} + inline TriangleAdjacent Adj() const ; + int swap(); + inline void SetAdj2(const TriangleAdjacent& , int =0); + inline Vertex * EdgeVertex(const int &) const ; + inline Vertex * OppositeVertex() const ; + inline Icoor2 & det() const; + inline int Locked() const ; + inline int GetAllFlag_UnSwap() const ; + inline void SetLock(); + inline int MarkUnSwap() const; + inline void SetMarkUnSwap(); + inline void SetCracked(); + inline int Cracked() const ; +};// end of Vertex class + + +///////////////////////////////////////////////////////////////////////////////////// +class Edge { public: + Vertex * v[2]; + Int4 ref; + GeometricalEdge * on; + 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] no continuite + GeometricalEdge * Adj [2]; + int SensAdj[2]; +// private: + int flag ; + public: + GeometricalEdge * link; // if Cracked() or Equi() + +// end of data + + 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 &); + + 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;} + void SetCracked() { flag |= 1;} + void SetDup() { flag |= 32;} // not a real edge + void SetEqui() { flag |= 2;} + void SetTgA() { flag|=4;} + void SetTgB() { flag|=8;} + void SetMark() { flag|=16;} + void SetUnMark() { flag &= 1007 /* 1023-16*/;} + void SetRequired() { flag |= 64;} + void SetReverseEqui() {flag |= 128;} + + inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); + +}; + +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 Triangle { + friend class TriangleAdjacent; + + + private: // les arete sont opposes a un sommet + Vertex * ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer + Triangle * at [3]; // nu triangle adjacent + Int1 aa[3]; // les nu des arete dans le triangles (mod 4) + public: + Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) + union { + Triangle * link ; + Int4 color; + }; + void SetDet() { + if(ns[0] && ns[1] && ns[2]) det = bamg::det(*ns[0],*ns[1],*ns[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 ns[0]==v || ns[1]==v || ns[2]==v ;} + TriangleAdjacent FindBoundaryEdge(int ) const; + + void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu) + { + if (link >=tb && link =tb && at[0] =tb && at[1] =tb && at[2] =vb && ns[0] =vb && ns[1] =vb && ns[2] at[aatt]=this; + tt->aa[aatt]=a + (aa[a] & 60 ) ;}// Copy all the mark + } + + void SetAdj2(Int1 a,Triangle *t,Int1 aat) + { at[a]=t;aa[a]=aat; + if(t) {t->at[aat]=this;t->aa[aat]=a;} + } + + void SetTriangleContainingTheVertex() + { + if (ns[0]) (ns[0]->t=this,ns[0]->vint=0); + if (ns[1]) (ns[1]->t=this,ns[1]->vint=1); + if (ns[2]) (ns[2]->t=this,ns[2]->vint=2); + } + + int swap(Int2 a1,int=0); + Int4 Optim(Int2 a,int =0); + + int Locked(int a)const { return aa[a]&4;} + int Hidden(int a)const { return aa[a]&16;} + int Cracked(int a) const { return aa[a] & 32;} + // for optimisation + int GetAllflag(int a){return aa[a] & 1020;} + void SetAllFlag(int a,int f){aa[a] = (aa[a] &3) + (1020 & f);} + + void SetHidden(int a){ + register Triangle * t = at[a]; + if(t) t->aa[aa[a] & 3] |=16; + aa[a] |= 16;} + void SetCracked(int a){ + register Triangle * t = at[a]; + if(t) t->aa[aa[a] & 3] |=32; + aa[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){ + register Triangle * t = at[a]; + t->aa[aa[a] & 3] |=4; + aa[a] |= 4;} + + void SetMarkUnSwap(int a){ + register Triangle * t = at[a]; + t->aa[aa[a] & 3] |=8; + aa[a] |=8 ;} + + + void SetUnMarkUnSwap(int a){ + register Triangle * t = at[a]; + t->aa[aa[a] & 3] &=55; // 23 + 32 + aa[a] &=55 ;} + +}; // end of Triangle class + + + + +class ListofIntersectionTriangles { +///////////////////////////////////////////////////////////////////////////////////// +class IntersectionTriangles { +public: + Triangle *t; + Real8 bary[3]; // use if t != 0 + R2 x; + Metric m; + Real8 s;// abscisse curviline + Real8 sp; // len of the previous seg in m + Real8 sn;// len of the next seg in m +}; +///////////////////////////////////////////////////////////////////////////////////// +class SegInterpolation { + public: + GeometricalEdge * e; + Real8 sBegin,sEnd; // abscisse of the seg on edge parameter + Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length + int last;// last index in ListofIntersectionTriangles for this Sub seg of edge + R2 F(Real8 s){ + Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01; + if (lBegin>s || s>lEnd){ + throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd")); + } + 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) + { + long int verbosity=0; + + if (NbSeg>=MaxNbSeg) { + int mneo= MaxNbSeg; + MaxNbSeg *= 2; + if (verbosity>3){ + printf(" reshape lSegsI from %i to %i\n",mneo,MaxNbSeg); + } + SegInterpolation * lEn = new SegInterpolation[MaxNbSeg]; + if (!lSegsI || NbSeg>=MaxNbSeg){ + throw ErrorException(__FUNCT__,exprintf("!lSegsI || NbSeg>=MaxNbSeg")); + } + for (int i=0;i< NbSeg;i++) + lEn[i] = lSegsI[MaxNbSeg]; // copy old to new + delete [] lSegsI; // remove old + lSegsI = lEn; + } + if (NbSeg) + lSegsI[NbSeg-1].last=Size; + lSegsI[NbSeg].e=e; + lSegsI[NbSeg].sBegin=s0; + lSegsI[NbSeg].sEnd=s1; + NbSeg++; + } + +// void CopyMetric(int i,int j){ lIntTria[j].m=lIntTria[i].m;} +// void CopyMetric(const Metric & mm,int j){ lIntTria[j].m=mm;} + + void ReShape() { + register int newsize = MaxSize*2; + IntersectionTriangles* nw = new IntersectionTriangles[newsize]; + if (!nw){ + throw ErrorException(__FUNCT__,exprintf("!nw")); + } + for (int i=0;i3) printf(" ListofIntersectionTriangles ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg); + MaxSize = newsize; + delete [] lIntTria;// remove old + lIntTria = nw; // copy pointer + } + + void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0); + Real8 Length(); + Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx); +}; + + +///////////////////////////////////////////////////////////////////////////////////// +class GeometricalSubDomain { +public: + GeometricalEdge *edge; + int sens; // -1 or 1 + Int4 ref; + inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &); + +}; +///////////////////////////////////////////////////////////////////////////////////// +class SubDomain { +public: + Triangle * head; + Int4 ref; + int sens; // -1 or 1 + Edge * edge; // to geometrical + inline void Set(const Triangles &,Int4,Triangles &); + +}; +///////////////////////////////////////////////////////////////////////////////////// +class VertexOnGeom { public: + + Vertex * mv; + Real8 abscisse; + union{ + GeometricalVertex * gv; // if abscisse <0; + GeometricalEdge * ge; // if abscisse in [0..1] + }; + inline void Set(const VertexOnGeom&,const Triangles &,Triangles &); + int OnGeomVertex()const {return this? abscisse <0 :0;} + int OnGeomEdge() const {return this? abscisse >=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->on=this;mv->vint=IsVertexOnGeom;} + inline void Set(const Triangles &,Int4,Triangles &); + +}; +///////////////////////////////////////////////////////////////////////////////////// +class VertexOnVertex {public: + Vertex * v, *bv; + VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){} + VertexOnVertex() {}; + inline void Set(const Triangles &,Int4,Triangles &); + void SetOnBTh(){v->onbv=bv;v->vint=IsVertexOnVertex;} +}; +///////////////////////////////////////////////////////////////////////////////////// +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->onbe=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;} +}; + + inline TriangleAdjacent FindTriangleAdjacent(Edge &E); + inline Vertex * TheVertex(Vertex * a); // for remove crak in mesh +///////////////////////////////////////////////////////////////////////////////////// + +class CrackedEdge { // a small class to store on crack an uncrack information + friend class Triangles; + class CrackedTriangle { + friend class Triangles; + friend class CrackedEdge; + Triangle * t; // edge of triangle t + int i; // edge number of in triangle + Edge *edge; // the 2 edge + Vertex *New[2]; // new vertex number + 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;} + void Crack(){ + 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) = 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 + + } + + }; // end of class CrackedTriangle + 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();} +}; + +///////////////////////////////////////////////////////////////////////////////////// +class Triangles { +public: + + enum TypeFileMesh { + AutoMesh=0,BDMesh=1,NOPOMesh=2,amMesh=3,am_fmtMesh=4,amdbaMesh=5, + ftqMesh=6,mshMesh=7}; + + int static counter; // to kown the number of mesh in memory + int OnDisk; // true if on disk + 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 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 NbOfTriangleSearchFind; + Int4 NbOfSwapTriangle; + char * name, *identity; + Vertex * vertices; // data of vertices des sommets + + 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; + + QuadTree *quadtree; + 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) ; + Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts); + + Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) + :Gh(BT.Gh),BTh(BT) { + try {GeomToTriangles1(nbvx,keepBackVertices);} + catch(...) { this->~Triangles(); throw; } } + + Triangles(Int4 nbvx,Geometry & G) + :Gh(G),BTh(*this){ + try { GeomToTriangles0(nbvx);} + 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); // truncature + + void SetIntCoor(const char * from =0); + + // void RandomInit(); + // void CubeInit(int ,int); + + Real8 MinimalHmin() {return 2.0/coefIcoor;} + Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} + 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];}; + I2 toI2(const R2 & P) const { + return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) + ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} + R2 toR2(const I2 & P) const { + return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);} + void Add( Vertex & s,Triangle * t,Icoor2 * =0) ; + void Insert(); + // void InsertOld(); + void ForceBoundary(); + void Heap(); + void FindSubDomain(int ); + Int4 ConsRefTriangle(Int4 *) const; + void ShowHistogram() const; + void ShowRegulaty() const; // Add FH avril 2007 +// void ConsLinkTriangle(); + + void ReMakeTriangleContainingTheVertex(); + void UnMarkUnSwapTriangle(); + void SmoothMetric(Real8 raisonmax) ; + void BoundAnisotropy(Real8 anisomax,double hminaniso= 1e-100) ; + void MaxSubDivision(Real8 maxsubdiv); + Edge** MakeGeometricalEdgeToEdge(); + void SetVertexFieldOn(); + void SetVertexFieldOnBTh(); + Int4 SplitInternalEdgeWithBorderVertices(); + void MakeQuadrangles(double costheta); + int SplitElement(int choice); + void MakeQuadTree(); + void NewPoints( Triangles &,int KeepBackVertex =1 ); + Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; + void NewPoints(int KeepBackVertex=1){ NewPoints(*this,KeepBackVertex);} + void ReNumberingTheTriangleBySubDomain(bool justcompress=false); + void ReNumberingVertex(Int4 * renu); + 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); + + Int4 Number(const Triangle & t) const { return &t - triangles;} + Int4 Number(const Triangle * t) const { return t - triangles;} + Int4 Number(const Vertex & t) const { return &t - vertices;} + Int4 Number(const Vertex * t) const { return t - vertices;} + Int4 Number(const Edge & t) const { return &t - edges;} + 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 * FindTriangleContening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; + + 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); + void IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols, + const Real8 hmin,const Real8 hmax, const Real8 coef, + const Real8 anisomax,const Real8 CutOff=1.e-4,const int NbJacobi=1, + const int DoNormalisation=1, + const double power=1.0, + const int choise=0); + void IntersectGeomMetric(const Real8 err,const int iso); + + + int isCracked() const {return NbCrackedVertices != 0;} + int Crack(); + int UnCrack(); + + void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo + void FillHoleInMesh() ; + int CrackMesh(); + private: + void GeomToTriangles1(Int4 nbvx,int KeepBackVertices=1);// the real constructor mesh adaption + void GeomToTriangles0(Int4 nbvx);// the real constructor mesh generator + void PreInit(Int4,char * =0 ); + +}; // End Class Triangles + + +///////////////////////////////////////////////////////////////////////////////////// +class Geometry { +public: + int OnDisk; + Int4 NbRef; // counter of ref on the this class if 0 we can delete + + char *name; + Int4 nbvx,nbtx; // nombre max de sommets , de Geometry + Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux, + Int4 NbSubDomains; // + Int4 NbEquiEdges; + Int4 NbCrackedEdges; +// Int4 nbtf;// de triangle frontiere + Int4 NbOfCurves; + int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); } + GeometricalVertex * vertices; // data of vertices des sommets + 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 + + R2 pmin,pmax; // extrema + Real8 coefIcoor; // coef to integer Icoor1; + Real8 MaximalAngleOfCorner; + +// end of data + + + 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);} + void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); + void EmptyGeometry(); + Geometry() {EmptyGeometry();}// empty Geometry + void AfterRead(); + Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();OnDisk=1;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;} + Int4 Number(const GeometricalEdge & t) const { return &t - edges;} + Int4 Number(const GeometricalEdge * t) const { return t - edges;} + Int4 Number(const Curve * c) const { return c - curves;} + + void UnMarkEdges() { + for (Int4 i=0;i=0) { + t = at[a]; + // if (t-this<0) return 0; + v2 = ns[VerticesOfTriangularEdge[a][0]]; + v0 = ns[VerticesOfTriangularEdge[a][1]]; + v1 = ns[OppositeEdge[a]]; + v3 = t->ns[OppositeEdge[aa[a]&3]]; + } + } + return t; +} + +inline double Triangle::QualityQuad(int a,int option) const +{ // first do the logique part + double q; + if (!link || aa[a] &4) + q= -1; + else { + Triangle * t = at[a]; + if (t-this<0) q= -1;// because we do 2 times + else if (!t->link ) q= -1; + else if (aa[0] & 16 || aa[1] & 16 || aa[2] & 16 || t->aa[0] & 16 || t->aa[1] & 16 || t->aa[2] & 16 ) + q= -1; + else if(option) + { + const Vertex & v2 = *ns[VerticesOfTriangularEdge[a][0]]; + const Vertex & v0 = *ns[VerticesOfTriangularEdge[a][1]]; + const Vertex & v1 = *ns[OppositeEdge[a]]; + const Vertex & v3 = * t->ns[OppositeEdge[aa[a]&3]]; + q = QuadQuality(v0,v1,v2,v3); // do the float part + } + else q= 1; + } + return q; +} + + +inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ) + { + *this = rec; + } +inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ) + { + *this = rec; + } +inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) + { + *this = Th.edges[i]; + v[0] = ThNew.vertices + Th.Number(v[0]); + v[1] = ThNew.vertices + Th.Number(v[1]); + if (on) + on = ThNew.Gh.edges+Th.Gh.Number(on); + if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]); + if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]); + + } +inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew) + { + *this = rec; + v[0] = GhNew.vertices + Gh.Number(v[0]); + v[1] = GhNew.vertices + Gh.Number(v[1]); + if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]); + if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]); + } + +inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew) +{ + *this = rec; + be = GhNew.edges + Gh.Number(be); + ee = GhNew.edges + Gh.Number(ee); + if(next) next= GhNew.curves + Gh.Number(next); +} + +inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew) + { + *this = rec; + if ( ns[0] ) ns[0] = ThNew.vertices + Th.Number(ns[0]); + if ( ns[1] ) ns[1] = ThNew.vertices + Th.Number(ns[1]); + if ( ns[2] ) ns[2] = ThNew.vertices + Th.Number(ns[2]); + if(at[0]) at[0] = ThNew.triangles + Th.Number(at[0]); + if(at[1]) at[1] = ThNew.triangles + Th.Number(at[1]); + if(at[2]) at[2] = ThNew.triangles + Th.Number(at[2]); + if (link >= Th.triangles && link < Th.triangles + Th.nbt) + link = ThNew.triangles + Th.Number(link); + } +inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) +{ + *this = Th.VertexOnBThVertex[i]; + v = ThNew.vertices + Th.Number(v); + +} +inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) +{ + *this = Th.subdomains[i]; + if ( head-Th.triangles<0 || head-Th.triangles>=Th.nbt){ + throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt")); + } + head = ThNew.triangles + Th.Number(head) ; + if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{ + throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe")); + } + edge = ThNew.edges+ Th.Number(edge); +} + inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew) +{ + *this = rec; + edge = Gh.Number(edge) + GhNew.edges; +} +inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) +{ + *this = Th.VertexOnBThEdge[i]; + v = ThNew.vertices + Th.Number(v); +} + +inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew) +{ + *this = rec; + mv = ThNew.vertices + Th.Number(mv); + if (gv) + if (abscisse < 0 ) + gv = ThNew.Gh.vertices + Th.Gh.Number(gv); + else + ge = ThNew.Gh.edges + Th.Gh.Number(ge); + +} +inline Real8 Edge::MetricLength() const + { + return LengthInterpole(v[0]->m,v[1]->m,v[1]->r - v[0]->r) ; + } + +inline void Triangles::ReMakeTriangleContainingTheVertex() + { + register Int4 i; + for ( i=0;iat[a]=ta.t; + t->aa[a]=ta.a|l;} + if(ta.t) { + ta.t->at[ta.a] = t ; + ta.t->aa[ta.a] = a| l ; + } +} + + +inline int TriangleAdjacent::Locked() const +{ return t->aa[a] &4;} +inline int TriangleAdjacent::Cracked() const +{ return t->aa[a] &32;} +inline int TriangleAdjacent::GetAllFlag_UnSwap() const +{ return t->aa[a] & 1012;} // take all flag except MarkUnSwap + +inline int TriangleAdjacent::MarkUnSwap() const +{ return t->aa[a] &8;} + +inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);} + +inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);} + +inline TriangleAdjacent TriangleAdjacent::Adj() const +{ return t->Adj(a);} + +inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const + {return t->ns[VerticesOfTriangularEdge[a][i]]; } +inline Vertex * TriangleAdjacent::OppositeVertex() const +{return t->ns[bamg::OppositeVertex[a]]; } +inline Icoor2 & TriangleAdjacent::det() const +{ return t->det;} +inline TriangleAdjacent Adj(const TriangleAdjacent & a) +{ return a.Adj();} + +inline TriangleAdjacent Next(const TriangleAdjacent & ta) +{ return TriangleAdjacent(ta.t,NextEdge[ta.a]);} + +inline TriangleAdjacent Previous(const TriangleAdjacent & ta) +{ return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);} + +inline void Adj(GeometricalEdge * & on,int &i) + {int j=i;i=on->SensAdj[i];on=on->Adj[j];} + +inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc) +{ + Real4 ret; + I2 ia=va,ib=vb,ic=vc; + I2 ab=ib-ia,bc=ic-ib,ac=ic-ia; + Icoor2 deta=Det(ab,ac); + if (deta <=0) ret = -1; + else { + Real8 a = sqrt((Real8) (ac,ac)), + b = sqrt((Real8) (bc,bc)), + c = sqrt((Real8) (ab,ab)), + p = a+b+c; + Real8 h= Max(Max(a,b),c),ro=deta/p; + ret = ro/h;} + return ret; +} + + +inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) { + Vertex *v=Th->vertices; + Int4 nbv = Th->nbv; + if (i<0 || j<0 || k<0){ + throw ErrorException(__FUNCT__,exprintf("i<0 || j<0 || k<0")); + } + if (i>=nbv || j>=nbv || k>=nbv){ + throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv")); + } + ns[0]=v+i; + ns[1]=v+j; + ns[2]=v+k; + at[0]=at[1]=at[2]=0; + aa[0]=aa[1]=aa[2]=0; + det=0; +} + +inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){ + ns[0]=v0; + ns[1]=v1; + ns[2]=v2; + at[0]=at[1]=at[2]=0; + aa[0]=aa[1]=aa[2]=0; + if (v0) det=0; + else { + det=-1; + link=NULL;}; +} + +inline Real4 Triangle::qualite() +{ + return det < 0 ? -1 : bamg::qualite(*ns[0],*ns[1],*ns[2]); +} + +Int4 inline Vertex::Optim(int i,int koption) +{ + Int4 ret=0; + if ( t && (vint >= 0 ) && (vint <3) ) + { + ret = t->Optim(vint,koption); + if(!i) + { + t =0; // for no future optime + vint= 0; } + } + return ret; +} + +Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c) +{ + register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y; + register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y; + return bax*cay - bay*cax;} + + +void swap(Triangle *t1,Int1 a1, + Triangle *t2,Int1 a2, + Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2); + + + +int inline TriangleAdjacent::swap() +{ return t->swap(a);} + + + +int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , + TriangleAdjacent & tt1,Icoor2 & dets1, + Icoor2 & detsa,Icoor2 & detsb, int & nbswap); + +int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ; + +// inline bofbof FH +inline TriangleAdjacent FindTriangleAdjacent(Edge &E) + { + Vertex * a = E.v[0]; + Vertex * b = E.v[1]; + + Triangle * t = a->t; + int i = a->vint; + TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge + if (!t || i<0 || i>=3){ + throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); + } + if ( a!=(*t)(i)){ + throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); + } + int k=0; + do { // turn around vertex in direct sens (trigo) + k++; + if (k>=20000){ + throw ErrorException(__FUNCT__,exprintf("k>=20000")); + } + // in no crack => ta.EdgeVertex(1) == a otherwise ??? + if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find + ta = ta.Adj(); + if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find + --ta; + } while (t != (Triangle *)ta); + throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found")); + return TriangleAdjacent(0,0);//for compiler + } + +inline Vertex * TheVertex(Vertex * a) // give a unique vertex with smallest number +{ // in case on crack in mesh + Vertex * r(a), *rr; + Triangle * t = a->t; + int i = a->vint; + TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge + if (!t || i<0 || i>=3){ + throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); + } + if ( a!=(*t)(i)){ + throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); + } + int k=0; + do { // turn around vertex in direct sens (trigo) + k++; + if (k>=20000){ + throw ErrorException(__FUNCT__,exprintf("k>=20000")); + } + // in no crack => ta.EdgeVertex(1) == a + if ((rr=ta.EdgeVertex(0)) < r) r = rr; + ta = ta.Adj(); + if ((rr=ta.EdgeVertex(1)) < r) r =rr; + --ta; + } while (t != (Triangle*) ta); + return r; +} + +inline double CPUtime(){ +#ifdef SYSTIMES + struct tms buf; + if (times(&buf)!=-1) + return ((double)buf.tms_utime+(double)buf.tms_stime)/(long) sysconf(_SC_CLK_TCK); + else +#endif + return ((double) clock())/CLOCKS_PER_SEC; +} + +} +#endif Index: /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2865) @@ -3,5 +3,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 2864) +++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 2865) @@ -4,5 +4,5 @@ #include -#include "../Bamgx.h" +#include "../Mesh2.h" #include "../QuadTree.h" #include "../SetOfE4.h" Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 2864) +++ /issm/trunk/src/c/Makefile.am (revision 2865) @@ -320,4 +320,5 @@ ./Bamgx/Bamgx.h\ ./Bamgx/Bamgx.cpp\ + ./Bamgx/Mesh2.h \ ./Bamgx/meshtype.h \ ./Bamgx/objects/MatVVP2x2.cpp \ @@ -660,4 +661,5 @@ ./Bamgx/Bamgx.h\ ./Bamgx/Bamgx.cpp\ + ./Bamgx/Mesh2.h \ ./Bamgx/meshtype.h \ ./Bamgx/objects/MatVVP2x2.cpp \