00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #ifndef __OPCLASS_H
00038 #ifdef __GNUG__
00039 #pragma interface
00040 #endif
00041 #define __OPCLASS_H
00042
00043 #include <cmath>
00044 #include <iostream>
00045
00046
00047 #define mmax(a,b)(a>b?a:b)
00048 #define mmin(a,b)(a<b?a:b)
00049 #define ssqr(x) ((x)*(x))
00050
00051
00052 namespace fem
00053 {
00054 extern int N;
00055 extern int N2;
00056
00057 class Complex;
00058 class cvect;
00059 class cmat;
00060
00061 #define sizecmat sizeof(cmat)
00062 #define sizecvect sizeof(cvect)
00063 #define deter2(a11,a12,a21,a22) (a11*a22 - a12*a21)
00064
00065 const float epsilon =(float)1.0e-20;
00066
00067
00068 typedef float ccreal;
00069 typedef Complex creal;
00070
00071 typedef Complex ccomplex;
00072 #define sizeccomplex sizeof(ccomplex)
00073
00074 #define sizecreal sizeof(creal)
00075 #define sizeccreal sizeof(ccreal)
00076
00077 void cgauss( cmat& a, cvect& b);
00078
00079 float norm2(const float& a);
00080 float imagpart(const float& a);
00081 float realpart(const float& a);
00082 cmat id(const cvect& a);
00083 Complex id(const Complex& x);
00084 void myassert(int what);
00085
00086 class Complex
00087
00088 {
00089 protected:
00090 float re,im;
00091 public:
00092 Complex() { re = 0.F; im =0.F; }
00093 Complex(float r, float i =0.F) {re =r; im =i;}
00094
00095 float& real () { return re;};
00096 float real() const { return re;}
00097 float& imag() { return im;};
00098 float imag() const { return im;}
00099
00100 Complex conjug() {return Complex(re,-im);}
00101
00102 Complex& operator=(const Complex& a) { re =a.re; im = a.im; return *this; }
00103
00104 Complex& operator*=(const Complex& a) { re =a.re*re - a.im*im;im= a.re*im + a.im*re; return *this; }
00105 Complex& operator/=(const Complex& a) {
00106 float xl=ssqr(a.re) + ssqr(a.im); re =(a.re*re+a.im*im)/xl; im =(a.re*im-a.im*re)/xl;
00107 return *this;
00108 }
00109 Complex& operator+=(const Complex& a) { re +=a.re; im += a.im; return *this; }
00110 Complex& operator-=(const Complex& a) { re -=a.re; im -= a.im; return *this; }
00111 Complex& operator*=(const float a) { re *= a;im *= a; return *this; }
00112 Complex& operator/=(const float a) { re /= a;im /= a; return *this; }
00113 Complex& operator+=(const float a) { re +=a; return *this; }
00114 Complex& operator-=(const float a) { re -=a; return *this; }
00115 Complex& operator-() { re =-re; im =-im; return *this; }
00116 float modul2() const {return ssqr(re)+ssqr(im);}
00117 float arg() const;
00118
00119 friend Complex operator *(const Complex a, const Complex b) {return Complex(a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re);}
00120 friend Complex operator *(const Complex a, const float b) { return Complex(a.re*b, a.im*b); }
00121 friend Complex operator *(const float b, const Complex a) { return Complex(a.re*b, a.im*b); }
00122 friend Complex operator /(const Complex a, const float b) { return Complex(a.re/b, a.im/b); }
00123 friend Complex operator /(const float b, const Complex a) { return Complex(b,0.F)/a; }
00124 friend Complex operator /(const Complex a, const Complex b) {
00125 Complex c; float xl=ssqr(b.re) + ssqr(b.im); c.re =(a.re*b.re+a.im*b.im)/xl; c.im =(b.re*a.im-b.im*a.re)/xl;
00126 return c;
00127 }
00128 friend Complex operator +(const Complex a, const Complex b) {return Complex(a.re+b.re, a.im+b.im);}
00129 friend Complex operator +(const Complex a, const float b) { return Complex(a.re+b, a.im); }
00130 friend Complex operator +(const float b, const Complex a) { return Complex(a.re+b, a.im); }
00131 friend Complex operator -(const Complex a, const Complex b) {return Complex(a.re-b.re, a.im-b.im);}
00132 friend Complex operator -(const Complex a, const float b) { return Complex(a.re-b, a.im); }
00133 friend Complex operator -(const float b, const Complex a) { return Complex(b-a.re, -a.im); }
00134 friend float norm2( const Complex& a) { return a.modul2(); }
00135 friend float realpart(const Complex& a){ return a.re; }
00136 friend float imagpart(const Complex& a){ return a.im; }
00137 friend std::ostream& operator<<(std::ostream& os, const Complex& a);
00138 friend std::istream& operator>>(std::istream& os, Complex& a);
00139 friend Complex id(const Complex& x);
00140
00141 friend Complex pow(const Complex&, const float&);
00142 friend float cos(const Complex&);
00143 friend float sin(const Complex&);
00144
00145 };
00146 inline float
00147 Complex::arg() const
00148 {
00149 if (modul2() > 1e-8)
00150 if (im > 0)
00151 return acos(re/sqrt(modul2()));
00152 else
00153 return 8*atan(1.) - acos(re/sqrt(modul2()));
00154 else
00155 return 0.;
00156
00157 }
00158 inline Complex
00159 pow(const Complex& c, const float& f)
00160 {
00161 Complex z(std::cos(f*c.arg()), std::sin(f*c.arg()));
00162 return std::pow(std::sqrt(c.modul2()),f)*z;
00163 }
00164 inline float
00165 cos(const Complex& c)
00166 {
00167 return std::cos(c.real())*std::cosh(c.imag()) + std::sin(c.real())*std::sinh(c.imag());\
00168 }
00169 inline float
00170 sin(const Complex& c)
00171 {
00172 return std::sin(c.real())*std::cosh(c.imag()) - std::cos(c.real())*std::sinh(c.imag());
00173 }
00174
00175 class cvect {public: ccreal val[2];
00176
00177 cvect() { val[0]=0.F; val[1]=0.F; }
00178 cvect(const cvect& r) { val[0]=r.val[0];val[1]=r.val[1];}
00180 cvect(ccreal r) { val[0]=r;val[1]=r;}
00181 ccreal& operator[] ( int i) { return val[i];}
00182 const ccreal& operator[] ( int i) const { return val[i];}
00183 float modul2() { return norm2(val[0])+norm2(val[1]);}
00184 ccreal operator*=(const cvect& a) { return val[0]*a.val[0] + val[1]*a.val[1]; }
00185 cvect& operator*=(const ccreal a) { val[0] *= a; val[1] *= a; return *this; }
00186 friend ccreal operator *(const cvect& a, const cvect& b) { cvect c(a); return c *= b; }
00187 friend cvect operator *(const cvect& a, const ccreal b) { cvect c(a); return c *= b; }
00188 friend cvect operator *(const ccreal b, const cvect& a) { cvect c(a); return c *= b; }
00189 cvect& operator/=(const ccreal a) {val[0] /= a;val[1] /= a; return *this; }
00190 friend cvect operator /(const cvect& a, const ccreal b) { cvect c(a); return c /= b; }
00191 cvect& operator+=(const cvect& a) { val[0] += a.val[0];val[1] += a.val[1]; return *this; }
00192 cvect& operator+=(const ccreal a) { val[0] += a;val[1] += a; return *this; }
00193 friend cvect operator +(const cvect& a, const cvect& b) { cvect c(a); return c += b; }
00194 friend cvect operator +(const cvect& a, const ccreal b) { cvect c(a); return c += b; }
00195 friend cvect operator +(const ccreal b, const cvect& a) { cvect c(a); return c += b; }
00196 cvect& operator-=(const cvect& a) { val[0] -= a.val[0];val[1] -= a.val[1]; return *this; }
00197 cvect& operator-=(const ccreal a) { val[0] -= a;val[1] -= a; return *this; }
00198 friend cvect operator -(const cvect& a, const cvect& b) { cvect c(a); return c -= b; }
00199 friend cvect operator -(const cvect& a, const ccreal b) { cvect c(a); return c -= b; }
00200 friend cvect operator -(const ccreal b, const cvect& a) { cvect c(b); return c -= a; }
00201 cvect& operator=(const cvect& a) { val[0] = a.val[0];val[1] = a.val[1]; return *this; }
00202 cvect& operator=(const ccreal a) {val[0] = a; val[1] = a; return *this; }
00203 cvect& operator-() { val[0] = -val[0]; val[1] = -val[1]; return *this; }
00204 friend float norm2( cvect& a) { return a.modul2();}
00205 friend float realpart(const cvect& a){ return realpart(a.val[0]); }
00206
00207 friend std::ostream& operator<<(std::ostream& os, cvect& a);
00208 friend std::istream& operator>>(std::istream& os, cvect& a);
00209 };
00210
00211
00212 class cmat {
00213
00214 public: ccreal val[4];
00215 cmat() {val[0]=0.F;val[1]=0.F;val[2]=0.F;val[3]=0.F; }
00216 cmat(const cmat& r) {val[0]=r.val[0];val[1]=r.val[1];val[2]=r.val[2];val[3]=r.val[3];}
00218 cmat(ccreal r) { val[0] = r;val[1] = r;val[2] = r;val[3] = r;}
00219 ccreal& operator() (int i, int j) { return val[i+i+j];}
00220 ccreal& operator[] (int i) { return val[i];}
00221 const ccreal& operator[] (int i) const { return val[i];}
00222 float modul2() {return val[0]+val[1]+val[2]+val[3];}
00223 cmat operator*=(const cmat& a)
00224 { cmat b;
00225 b.val[0] =val[0]*a.val[0]+val[1]*a.val[2];
00226 b.val[1] =val[0]*a.val[1]+val[1]*a.val[3];
00227 b.val[2] =val[2]*a.val[0]+val[3]*a.val[2];
00228 b.val[3] =val[2]*a.val[1]+val[3]*a.val[3];
00229 return b;
00230 }
00231 cmat& operator*=(const ccreal a) { val[0] *= a; val[1] *= a; val[2] *= a; val[3] *= a; return *this; }
00232 friend cmat operator *(const cmat& a, const cmat& b) { cmat c(a); return c *= b; }
00233 friend cmat operator *(const cmat& a, const ccreal b) { cmat c(a); return c *= b; }
00234 friend cmat operator *(const ccreal b, const cmat& a) { cmat c(a); return c *= b; }
00235 cmat& operator/=(const ccreal a) { val[0] /= a; val[1] /= a; val[2] /= a; val[3] /= a; return *this; }
00236 friend cmat operator /(const cmat& a, const ccreal b) { cmat c(a); return c /= b; }
00237 cmat& operator+=(const cmat& a){val[0]+=a.val[0];val[1]+=a.val[1];val[2]+=a.val[2];val[3]+=a.val[3]; return *this; }
00238 cmat& operator+=(const ccreal a) { val[0] += a; val[1] += a; val[2] += a; val[3] += a; return *this; }
00239 friend cmat operator +(const cmat& a, const cmat& b) { cmat c(a); return c += b; }
00240 friend cmat operator +(const cmat& a, const ccreal b) { cmat c(a); return c += b; }
00241 friend cmat operator +(const ccreal b, const cmat& a) { cmat c(a); return c += b; }
00242
00243 cmat& operator-=(const cmat& a){val[0]-=a.val[0];val[1]-=a.val[1];val[2]-=a.val[2];val[3]-=a.val[3]; return *this; }
00244 cmat& operator-=(const ccreal a) { val[0] -= a; val[2] -= a; val[1] -= a; val[3] -= a; return *this; }
00245 friend cmat operator -(const cmat& a, const cmat& b) { cmat c(a); return c -= b; }
00246 friend cmat operator -(const cmat& a, const ccreal b) { cmat c(a); return c -= b; }
00247 friend cmat operator -(const ccreal b, const cmat& a) { cmat c(b); return c -= a; }
00248 cmat& operator=(const cmat& a){val[0]=a.val[0];val[1]=a.val[1];val[2]=a.val[2]; val[3] = a.val[3]; return *this; }
00249 cmat& operator=(const ccreal a) { val[0] = a; val[1] = a; val[2] = a; val[3] = a; return *this; }
00250 cmat& operator-() { val[0] = -val[0];val[1] = -val[1];val[2] = -val[2];val[3] = -val[3]; return *this; }
00251 friend float norm2( cmat& a) { return a.modul2();}
00252
00253 friend cvect operator *(const cmat& a, const cvect& b)
00254 { cvect c; c.val[0] = a.val[0]*b.val[0]+a.val[1]*b.val[1];
00255 c.val[1] = a.val[2]*b.val[0]+a.val[3]*b.val[1];
00256 return c ; }
00257 friend cvect operator *(const cvect& b, const cmat& a)
00258 { cvect c; c.val[0] = a.val[0]*b.val[0]+a.val[2]*b.val[1];
00259 c.val[1] = a.val[1]*b.val[0]+a.val[3]*b.val[1];
00260 return c ; }
00261 friend cvect operator /(const cvect& b,const cmat& a)
00262 { cvect c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
00263 if(norm2(s)<epsilon) { s = epsilon;}
00264 c.val[0] = deter2(b.val[0],a.val[1], b.val[1],a.val[3]) / s;
00265 c.val[1] = deter2(a.val[0],b.val[0], a.val[2],b.val[1]) / s;
00266 return c;
00267 }
00268 friend cmat operator /(const cmat& b,const cmat& a)
00269 {
00270 cmat c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
00271 if(norm2(s)<epsilon){ s = epsilon;}
00272 c.val[0] = deter2(b.val[0],a.val[1], b.val[2],a.val[3]) / s;
00273 c.val[2] = deter2(a.val[0],b.val[0], a.val[2],b.val[2]) / s;
00274 c.val[1] = deter2(b.val[1],a.val[1], b.val[3],a.val[3]) / s;
00275 c.val[3] = deter2(a.val[0],b.val[1], a.val[2],b.val[3]) / s;
00276 return c;
00277 }
00278 friend std::ostream& operator<<(std::ostream& os, cmat& a);
00279 friend std::istream& operator>>(std::istream& os, cmat& a);
00280 friend cmat id(const cvect& a);
00281 };
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 class Acreal
00310 {
00311 public:
00312 long szz;
00313 creal* cc;
00314
00315 Acreal(long=0);
00316 Acreal(const Acreal&);
00317 ~Acreal() { delete [] cc;cc=0;szz = 0;}
00318 void destroy() {delete [] cc;cc=0;szz = 0;}
00319 creal& operator[] (long i) { return cc[i];}
00320 creal* operator&() { return cc;}
00321 void init(long newSize);
00322
00323
00324
00325 };
00326 class Aint
00327 {
00328 public:
00329 long szz;
00330 int* cc;
00331
00332 Aint(long=0);
00333 Aint(const Aint&);
00334 ~Aint() { delete [] cc;cc=0;szz = 0;}
00335 void destroy() {delete [] cc;cc=0;szz = 0;}
00336 int& operator[] (long i) { myassert((i< szz)&&(i>=0)); return cc[i];}
00337 int* operator&() { return cc;}
00338 void init (long);
00339 };
00340
00341
00342 class Acvect
00343 {
00344 public:
00345 long szz;
00346 cvect* cc;
00347
00348 Acvect(long=0);
00349 Acvect(const Acvect&);
00350 ~Acvect() { delete [] cc;cc=0;szz = 0;}
00351 void destroy() { delete [] cc;cc=0;szz = 0;}
00352 cvect& operator[] (long i) { return cc[i];}
00353 cvect* operator&() { return cc;}
00354 void init (long);
00355 };
00356
00357 class Acmat
00358 {
00359 public:
00360 long szz;
00361 cmat* cc;
00362
00363 Acmat(long=0);
00364 Acmat(const Acmat&);
00365 ~Acmat()
00366 {
00367 delete [] cc;
00368 cc=0;
00369 szz = 0;
00370 }
00371 void destroy() {delete [] cc;cc=0;szz = 0;}
00372 cmat& operator[] (long i) { return cc[i];}
00373 cmat* operator&() { return cc;}
00374 void init (long);
00375 };
00376
00377 class AAcmat
00378 {
00379 public:
00380 long szz;
00381 Acmat* cc;
00382
00383 AAcmat(long=0);
00384 AAcmat(const AAcmat&);
00385 ~AAcmat() { delete [] cc;cc=0;szz = 0;}
00386 void destroy() {delete [] cc;cc=0;szz = 0;}
00387 Acmat& operator[] (long i) { return cc[i];}
00388 Acmat* operator&() { return cc;}
00389 void init(long );
00390 };
00391
00392 class AAcreal
00393 {
00394 public:
00395 long szz;
00396 Acreal* cc;
00397
00398 AAcreal(long=0);
00399 AAcreal(const AAcreal& );
00400 ~AAcreal() { delete [] cc;cc=0;szz = 0;}
00401 void destroy() {delete [] cc;cc=0;szz = 0;}
00402 Acreal& operator[] (long i) { return cc[i];}
00403 Acreal* operator&() { return cc;}
00404 void init (long);
00405 };
00406 }
00407
00408 #endif