Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2983) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2984) @@ -32,10 +32,10 @@ //some functions inline int BinaryRand(){ - #ifdef RAND_MAX +#ifdef RAND_MAX const long HalfRandMax = RAND_MAX/2; return rand() < HalfRandMax; - #else +#else return rand() & 16384; //2^14 (for sun because RAND_MAX is not def in stdlib.h) - #endif +#endif } @@ -45,6 +45,6 @@ typedef P2 R2; typedef P2 R2xR2; + } - #include "Metric.h" Index: /issm/trunk/src/c/Bamgx/Metric.h =================================================================== --- /issm/trunk/src/c/Bamgx/Metric.h (revision 2983) +++ /issm/trunk/src/c/Bamgx/Metric.h (revision 2984) @@ -6,5 +6,5 @@ namespace bamg { - typedef P2 D2; + typedef P2 D2; typedef P2xP2 D2xD2; @@ -13,61 +13,79 @@ - class Metric{ public: - friend class MatVVP2x2; - Real8 a11,a21,a22; - Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){} - Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){} - Metric() {}; // - Metric(const Real8 a[3],const Metric m0, - const Metric m1,const Metric m2 ); - R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);} - Real8 det() const {return a11*a22-a21*a21;} - R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} - R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} - // D2 Orthogonal(const D2 x){return D2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} - Metric( Real8 a,const Metric ma, - Real8 b,const Metric mb); - int IntersectWith(const Metric M2); - Metric operator*(Real8 c) const {Real8 c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);} - Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} - operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} + class Metric{ - Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; - // Real8 operator()(D2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; - Real8 operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;}; - inline void Box(Real4 &hx,Real4 &hy) const ; - Metric(const MatVVP2x2); - void Echo(); + public: + //fields + Real8 a11,a21,a22; + //friends + friend class MatVVP2x2; + //functions + Metric(){}; + Metric(const MatVVP2x2); + Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){} + Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){} + Metric( Real8 a,const Metric ma, Real8 b,const Metric mb); + Metric(const Real8 a[3],const Metric m0,const Metric m1,const Metric m2 ); + void Echo(); + R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);} + Real8 det() const {return a11*a22-a21*a21;} + R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} + R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} + int IntersectWith(const Metric M2); + inline void Box(Real4 &hx,Real4 &hy) const ; + //operators + Metric operator*(Real8 c) const {Real8 c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);} + Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} + operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} + Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; + Real8 operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;}; + }; class MatVVP2x2{ - friend class Metric; public: - double lambda1,lambda2; - D2 v; - - void Echo(); - MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){} - - void Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);} - void pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);} - void Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;} - void Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;} - - void Minh(double h) {Min(1.0/(h*h));} - void Maxh(double h) {Max(1.0/(h*h));} - void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);} - MatVVP2x2(const Metric ); - Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));} - Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));} - Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);} - Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);} - Real8 Aniso2() const { return lmax()/lmin();} - inline void BoundAniso2(const Real8 coef); - Real8 Aniso() const { return sqrt( Aniso2());} - void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));} - void operator *=(double coef){ lambda1*=coef;lambda2*=coef;} + //fields + double lambda1,lambda2; + D2 v; + //friends + friend class Metric; + //functions + MatVVP2x2(const Metric ); + MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){} + void Echo(); + void Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);} + void pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);} + void Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;} + void Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;} + void Minh(double h) {Min(1.0/(h*h));} + void Maxh(double h) {Max(1.0/(h*h));} + void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);} + Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));} + Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));} + Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);} + Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);} + Real8 Aniso2() const { return lmax()/lmin();} + Real8 Aniso() const { return sqrt( Aniso2());} + void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));} + inline void BoundAniso2(const Real8 coef); + //operators + void operator *=(double coef){ lambda1*=coef;lambda2*=coef;} }; + class SaveMetricInterpole { + friend Real8 LengthInterpole(const Metric ,const Metric , R2 ); + friend Real8 abscisseInterpole(const Metric ,const Metric , R2 ,Real8 ,int ); + int opt; + Real8 lab; + Real8 L[1024],S[1024]; + }; + + extern SaveMetricInterpole LastMetricInterpole; // for optimization + //Functions + void SimultaneousMatrixReduction( Metric M1, Metric M2,D2xD2 &V); + Real8 LengthInterpole(const Metric Ma,const Metric Mb, R2 AB); + Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8 s,int optim=0); + + //inlines inline void MatVVP2x2::BoundAniso2(const Real8 coef){ if (coef<=1.00000000001){ @@ -84,11 +102,5 @@ } } - - void SimultaneousMatrixReduction( Metric M1, Metric M2,D2xD2 &V); - inline Metric::Metric(const MatVVP2x2 M) { - // recompose M in: M = V^t lambda V - // V = ( v,v^\perp) - // where v^\perp = (-v_1,v_0) double v00=M.v.x*M.v.x; double v11=M.v.y*M.v.y; @@ -97,6 +109,5 @@ a21=v01*(M.lambda1-M.lambda2); a22=v00*M.lambda2+v11*M.lambda1; - } - + } inline void Metric::Box(Real4 &hx,Real4 &hy) const { Real8 d= a11*a22-a21*a21; @@ -104,23 +115,10 @@ hy = sqrt(a11/d); } - - class SaveMetricInterpole { - friend Real8 LengthInterpole(const Metric ,const Metric , R2 ); - friend Real8 abscisseInterpole(const Metric ,const Metric , R2 ,Real8 ,int ); - int opt; - Real8 lab; - Real8 L[1024],S[1024]; - }; - - extern SaveMetricInterpole LastMetricInterpole; // for optimization - - Real8 LengthInterpole(const Metric Ma,const Metric Mb, R2 AB); - Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8 s,int optim=0); - - inline Real8 LengthInterpole(Real8 la,Real8 lb) - { return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ? (la+lb)/2 : la*lb*log(la/lb)/(la-lb);} - - inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s) - { return ( Abs(la - lb) <1.0e-6*Max3(la,lb,1.0e-20)) ? s : (exp(s*lab*(la-lb)/(la*lb))-1)*lb/(la-lb);} + inline Real8 LengthInterpole(Real8 la,Real8 lb) { + return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ? (la+lb)/2 : la*lb*log(la/lb)/(la-lb); + } + inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s){ + return ( Abs(la - lb) <1.0e-6*Max3(la,lb,1.0e-20)) ? s : (exp(s*lab*(la-lb)/(la*lb))-1)*lb/(la-lb); + } } Index: /issm/trunk/src/c/Bamgx/QuadTree.h =================================================================== --- /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2983) +++ /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2984) @@ -5,11 +5,13 @@ namespace bamg { + typedef long IntQuad; + const int MaxDeep = 30; - typedef long IntQuad; - //long int 1, bitwise operation: 8L = 00001000 << 2L -> 00100000 shifted left by 2 const IntQuad MaxISize = ( 1L << MaxDeep); + class Triangles; class Vertex; - class QuadTree { + + class QuadTree{ private: class QuadTreeBox { @@ -39,9 +41,7 @@ public: - //fields QuadTreeBox* root; Triangles* th; - //functions ~QuadTree(); Index: /issm/trunk/src/c/Bamgx/R2.h =================================================================== --- /issm/trunk/src/c/Bamgx/R2.h (revision 2983) +++ /issm/trunk/src/c/Bamgx/R2.h (revision 2984) @@ -4,7 +4,11 @@ template class P2 { - public: + //fields R x,y; + //functions + P2 () :x(0),y(0) {}; + P2 (R a,R b) :x(a),y(b) {} + P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {} void Echo(){ printf("Member of P2:\n"); @@ -12,7 +16,5 @@ printf(" y: %g\n",y); } - P2 () :x(0),y(0) {}; - P2 (R a,R b) :x(a),y(b) {} - P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {} + //operators RR operator,(const P2 & cc) const {return (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} //scalar product P2 operator+(const P2 & cc) const {return P2(x+cc.x,y+cc.y);} @@ -29,37 +31,37 @@ template class P2xP2 { - friend P2 operator*(P2 c,P2xP2 cc){ - return P2(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y); - } + private: + friend P2 operator*(P2 c,P2xP2 cc){ + return P2(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y); + } public: - P2 x,y; - - //functions - void Echo(){ - printf("Member of P2xP2:\n"); - printf(" x.x: %g x.y: %g\n",x.x,x.y); - printf(" y.x: %g y.x: %g\n",y.x,y.y); - } - P2xP2 (): x(),y() {} - P2xP2 (P2 a,P2 b): x(a),y(b) {} - P2xP2 (P2 a,P2 b,P2 c ): x(b-a),y(c-a) {} - P2xP2 (R xx,R xy,R yx,R yy) :x(xx,xy),y(yx,yy) {} - RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;} - P2xP2 inv() const{ - RR d = (*this).det(); - return P2xP2((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) ); - }; - P2xP2 t() {return P2xP2(x.x,y.x,x.y,y.y);} //transposer - P2 tx() {return P2(x.x,y.x);} - P2 ty() {return P2(x.y,y.y);} - - //Operators - P2 operator*(const P2 c) const {return P2(x.x*c.x + x.y*c.y, y.x*c.x + y.y*c.y);} - P2xP2 operator*(P2xP2 c) const{ - return P2xP2(x.x*c.x.x + x.y*c.y.x, - x.x*c.x.y + x.y*c.y.y, - y.x*c.x.x + y.y*c.y.x, - y.x*c.x.y + y.y*c.y.y); - } + //fields + P2 x,y; + //functions + P2xP2 (): x(),y() {} + P2xP2 (P2 a,P2 b): x(a),y(b) {} + P2xP2 (P2 a,P2 b,P2 c ): x(b-a),y(c-a) {} + P2xP2 (R xx,R xy,R yx,R yy) :x(xx,xy),y(yx,yy) {} + void Echo(){ + printf("Member of P2xP2:\n"); + printf(" x.x: %g x.y: %g\n",x.x,x.y); + printf(" y.x: %g y.x: %g\n",y.x,y.y); + } + RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;} + P2xP2 inv() const{ + RR d = (*this).det(); + return P2xP2((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) ); + }; + P2xP2 t() {return P2xP2(x.x,y.x,x.y,y.y);} //transposer + P2 tx() {return P2(x.x,y.x);} + P2 ty() {return P2(x.y,y.y);} + //Operators + P2 operator*(const P2 c) const {return P2(x.x*c.x + x.y*c.y, y.x*c.x + y.y*c.y);} + P2xP2 operator*(P2xP2 c) const{ + return P2xP2(x.x*c.x.x + x.y*c.y.x, + x.x*c.x.y + x.y*c.y.y, + y.x*c.x.x + y.y*c.y.x, + y.x*c.x.y + y.y*c.y.y); + } }; @@ -69,30 +71,24 @@ return (RR) x.x * (RR) y.y - (RR) x.y * (RR) y.x ; } - template inline RR Area2 (const P2 a,const P2 b,const P2 c) { return Det(b-a,c-a) ; } - template inline R Norme1 (const P2 x) { return (Abs(x.x)+Abs(x.y)) ; } - template inline R NormeInfini (const P2 x) { return Max(Abs(x.x),Abs(x.y)) ; } - template inline RR Norme2_2 (const P2 x) { return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ; } - template inline RR Norme2 (const P2 x) { return sqrt((RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y) ; } - template inline P2 Orthogonal (const P2 x) { Index: /issm/trunk/src/c/Bamgx/SetOfE4.h =================================================================== --- /issm/trunk/src/c/Bamgx/SetOfE4.h (revision 2983) +++ /issm/trunk/src/c/Bamgx/SetOfE4.h (revision 2984) @@ -14,22 +14,24 @@ class SetOfEdges4 { - Int4 nx,nbax,NbOfEdges; - Int4* head; - Int4Edge* Edges; + + private: + Int4 nx,nbax,NbOfEdges; + Int4* head; + Int4Edge* Edges; public: - SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet - ~SetOfEdges4() {delete [] head; delete [] Edges;} - Int4 add (Int4 ii,Int4 jj); - Int4 SortAndAdd (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;} - Int4 nb(){return NbOfEdges;} - Int4 find (Int4 ii,Int4 jj); - Int4 SortAndFind (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;} - Int4 i(Int4 k){return Edges[k].i;} - Int4 j(Int4 k){return Edges[k].j;} - Int4 newarete(Int4 k){return NbOfEdges == k+1;} + SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet + ~SetOfEdges4() {delete [] head; delete [] Edges;} + Int4 add (Int4 ii,Int4 jj); + Int4 SortAndAdd (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;} + Int4 nb(){return NbOfEdges;} + Int4 find (Int4 ii,Int4 jj); + Int4 SortAndFind (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;} + Int4 i(Int4 k){return Edges[k].i;} + Int4 j(Int4 k){return Edges[k].j;} + Int4 newarete(Int4 k){return NbOfEdges == k+1;} - //operators - Int4Edge & operator[](Int4 k){return Edges[k];} + //operators + Int4Edge & operator[](Int4 k){return Edges[k];} }; } Index: /issm/trunk/src/c/Bamgx/meshtype.h =================================================================== --- /issm/trunk/src/c/Bamgx/meshtype.h (revision 2983) +++ /issm/trunk/src/c/Bamgx/meshtype.h (revision 2984) @@ -3,33 +3,36 @@ #include namespace bamg { -template inline T Square (const T &a) { return a*a;} -template inline T Min (const T &a,const T &b){return a < b ? a : b;} -template inline T Max (const T &a,const T & b){return a > b ? a : b;} -template inline T Abs (const T &a){return a <0 ? -a : a;} -template inline double Norme (const T &a){return sqrt(a*a);} -template inline void Exchange (T& a,T& b) {T c=a;a=b;b=c;} -// for pb on microsoft compiler -template inline T Max3 (const T &a,const T & b,const T & c){return Max(Max(a,b),c);} -template inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);} -typedef float Real4; -typedef double Real8; -typedef short Int1; -typedef short Int2; -typedef long Int4; + //template functions + template inline T Square (const T &a) { return a*a;} + template inline T Min (const T &a,const T &b){return a < b ? a : b;} + template inline T Max (const T &a,const T & b){return a > b ? a : b;} + template inline T Abs (const T &a){return a <0 ? -a : a;} + template inline double Norme (const T &a){return sqrt(a*a);} + template inline void Exchange (T& a,T& b) {T c=a;a=b;b=c;} + + // for pb on microsoft compiler + template inline T Max3 (const T &a,const T & b,const T & c){return Max(Max(a,b),c);} + template inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);} + + typedef float Real4; + typedef double Real8; + typedef short Int1; + typedef short Int2; + typedef long Int4; #if LONG_BIT > 63 -// for alpha and silicon - typedef int Icoor1; - typedef long Icoor2; - const Icoor1 MaxICoor = 1073741823; // 2^30-1 - const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; + // for alpha and silicon + typedef int Icoor1; + typedef long Icoor2; + const Icoor1 MaxICoor = 1073741823; // 2^30-1 + const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; #else - typedef int Icoor1; - typedef double Icoor2; - const Icoor1 MaxICoor = 8388608; // 2^23 - const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; + typedef int Icoor1; + typedef double Icoor2; + const Icoor1 MaxICoor = 8388608; // 2^23 + const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; #endif - class Triangles; + class Triangles; } #endif