![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2984 for issm/trunk/src/c/Bamgx/Metric.h
- Timestamp:
- 02/08/10 13:33:59 (15 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/Bamgx/Metric.h (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Metric.h
r2981 r2984 6 6 namespace bamg { 7 7 8 typedef P2<double,double> D2;8 typedef P2<double,double> D2; 9 9 typedef P2xP2<double,double> D2xD2; 10 10 … … 13 13 14 14 15 class Metric{ public: 16 friend class MatVVP2x2; 17 Real8 a11,a21,a22; 18 Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){} 19 Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){} 20 Metric() {}; // 21 Metric(const Real8 a[3],const Metric m0, 22 const Metric m1,const Metric m2 ); 23 R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);} 24 Real8 det() const {return a11*a22-a21*a21;} 25 R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 26 R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 27 // D2 Orthogonal(const D2 x){return D2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 28 Metric( Real8 a,const Metric ma, 29 Real8 b,const Metric mb); 30 int IntersectWith(const Metric M2); 31 Metric operator*(Real8 c) const {Real8 c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);} 32 Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} 33 operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} 15 class Metric{ 34 16 35 Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; 36 // Real8 operator()(D2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; 37 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;}; 38 inline void Box(Real4 &hx,Real4 &hy) const ; 39 Metric(const MatVVP2x2); 40 void Echo(); 17 public: 18 //fields 19 Real8 a11,a21,a22; 20 //friends 21 friend class MatVVP2x2; 22 //functions 23 Metric(){}; 24 Metric(const MatVVP2x2); 25 Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){} 26 Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){} 27 Metric( Real8 a,const Metric ma, Real8 b,const Metric mb); 28 Metric(const Real8 a[3],const Metric m0,const Metric m1,const Metric m2 ); 29 void Echo(); 30 R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);} 31 Real8 det() const {return a11*a22-a21*a21;} 32 R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 33 R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 34 int IntersectWith(const Metric M2); 35 inline void Box(Real4 &hx,Real4 &hy) const ; 36 //operators 37 Metric operator*(Real8 c) const {Real8 c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);} 38 Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} 39 operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} 40 Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; 41 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;}; 42 41 43 }; 42 44 43 45 class MatVVP2x2{ 44 friend class Metric;45 46 public: 46 double lambda1,lambda2; 47 D2 v; 48 49 void Echo(); 50 MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){} 51 52 void Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);} 53 void pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);} 54 void Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;} 55 void Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;} 56 57 void Minh(double h) {Min(1.0/(h*h));} 58 void Maxh(double h) {Max(1.0/(h*h));} 59 void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);} 60 MatVVP2x2(const Metric ); 61 Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));} 62 Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));} 63 Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);} 64 Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);} 65 Real8 Aniso2() const { return lmax()/lmin();} 66 inline void BoundAniso2(const Real8 coef); 67 Real8 Aniso() const { return sqrt( Aniso2());} 68 void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));} 69 void operator *=(double coef){ lambda1*=coef;lambda2*=coef;} 47 //fields 48 double lambda1,lambda2; 49 D2 v; 50 //friends 51 friend class Metric; 52 //functions 53 MatVVP2x2(const Metric ); 54 MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){} 55 void Echo(); 56 void Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);} 57 void pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);} 58 void Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;} 59 void Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;} 60 void Minh(double h) {Min(1.0/(h*h));} 61 void Maxh(double h) {Max(1.0/(h*h));} 62 void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);} 63 Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));} 64 Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));} 65 Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);} 66 Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);} 67 Real8 Aniso2() const { return lmax()/lmin();} 68 Real8 Aniso() const { return sqrt( Aniso2());} 69 void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));} 70 inline void BoundAniso2(const Real8 coef); 71 //operators 72 void operator *=(double coef){ lambda1*=coef;lambda2*=coef;} 70 73 }; 71 74 75 class SaveMetricInterpole { 76 friend Real8 LengthInterpole(const Metric ,const Metric , R2 ); 77 friend Real8 abscisseInterpole(const Metric ,const Metric , R2 ,Real8 ,int ); 78 int opt; 79 Real8 lab; 80 Real8 L[1024],S[1024]; 81 }; 82 83 extern SaveMetricInterpole LastMetricInterpole; // for optimization 84 //Functions 85 void SimultaneousMatrixReduction( Metric M1, Metric M2,D2xD2 &V); 86 Real8 LengthInterpole(const Metric Ma,const Metric Mb, R2 AB); 87 Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8 s,int optim=0); 88 89 //inlines 72 90 inline void MatVVP2x2::BoundAniso2(const Real8 coef){ 73 91 if (coef<=1.00000000001){ … … 84 102 } 85 103 } 86 87 void SimultaneousMatrixReduction( Metric M1, Metric M2,D2xD2 &V);88 89 104 inline Metric::Metric(const MatVVP2x2 M) { 90 // recompose M in: M = V^t lambda V91 // V = ( v,v^\perp)92 // where v^\perp = (-v_1,v_0)93 105 double v00=M.v.x*M.v.x; 94 106 double v11=M.v.y*M.v.y; … … 97 109 a21=v01*(M.lambda1-M.lambda2); 98 110 a22=v00*M.lambda2+v11*M.lambda1; 99 } 100 111 } 101 112 inline void Metric::Box(Real4 &hx,Real4 &hy) const { 102 113 Real8 d= a11*a22-a21*a21; … … 104 115 hy = sqrt(a11/d); 105 116 } 106 107 class SaveMetricInterpole { 108 friend Real8 LengthInterpole(const Metric ,const Metric , R2 ); 109 friend Real8 abscisseInterpole(const Metric ,const Metric , R2 ,Real8 ,int ); 110 int opt; 111 Real8 lab; 112 Real8 L[1024],S[1024]; 113 }; 114 115 extern SaveMetricInterpole LastMetricInterpole; // for optimization 116 117 Real8 LengthInterpole(const Metric Ma,const Metric Mb, R2 AB); 118 Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8 s,int optim=0); 119 120 inline Real8 LengthInterpole(Real8 la,Real8 lb) 121 { return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ? (la+lb)/2 : la*lb*log(la/lb)/(la-lb);} 122 123 inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s) 124 { 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);} 117 inline Real8 LengthInterpole(Real8 la,Real8 lb) { 118 return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ? (la+lb)/2 : la*lb*log(la/lb)/(la-lb); 119 } 120 inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s){ 121 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); 122 } 125 123 126 124 }
Note:
See TracChangeset
for help on using the changeset viewer.
