femMisc.hpp
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 void myassert(int what);
00083
00084 class Complex
00085
00086 {
00087 protected:
00088 float re,im;
00089 public:
00090 Complex() { re = 0.F; im =0.F; }
00091 Complex(float r, float i =0.F) {re =r; im =i;}
00092
00093 float& real () { return re;};
00094 float real() const { return re;}
00095 float& imag() { return im;};
00096 float imag() const { return im;}
00097
00098 Complex conjug() {return Complex(re,-im);}
00099
00100 Complex& operator=(const Complex& a) { re =a.re; im = a.im; return *this; }
00101
00102 Complex& operator*=(const Complex& a) { re =a.re*re - a.im*im;im= a.re*im + a.im*re; return *this; }
00103 Complex& operator/=(const Complex& a) {
00104 float xl=ssqr(a.re) + ssqr(a.im); re =(a.re*re+a.im*im)/xl; im =(a.re*im-a.im*re)/xl;
00105 return *this;
00106 }
00107 Complex& operator+=(const Complex& a) { re +=a.re; im += a.im; return *this; }
00108 Complex& operator-=(const Complex& a) { re -=a.re; im -= a.im; return *this; }
00109 Complex& operator*=(const float a) { re *= a;im *= a; return *this; }
00110 Complex& operator/=(const float a) { re /= a;im /= a; return *this; }
00111 Complex& operator+=(const float a) { re +=a; return *this; }
00112 Complex& operator-=(const float a) { re -=a; return *this; }
00113 Complex& operator-() { re =-re; im =-im; return *this; }
00114 float modul2() const {return ssqr(re)+ssqr(im);}
00115 float arg() const;
00116
00117 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);}
00118 friend Complex operator *(const Complex a, const float b) { return Complex(a.re*b, a.im*b); }
00119 friend Complex operator *(const float b, const Complex a) { return Complex(a.re*b, a.im*b); }
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(b,0.F)/a; }
00122 friend Complex operator /(const Complex a, const Complex b) {
00123 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;
00124 return c;
00125 }
00126 friend Complex operator +(const Complex a, const Complex b) {return Complex(a.re+b.re, a.im+b.im);}
00127 friend Complex operator +(const Complex a, const float b) { return Complex(a.re+b, a.im); }
00128 friend Complex operator +(const float b, const Complex a) { return Complex(a.re+b, a.im); }
00129 friend Complex operator -(const Complex a, const Complex b) {return Complex(a.re-b.re, a.im-b.im);}
00130 friend Complex operator -(const Complex a, const float b) { return Complex(a.re-b, a.im); }
00131 friend Complex operator -(const float b, const Complex a) { return Complex(b-a.re, -a.im); }
00132 friend float norm2( const Complex& a) { return a.modul2(); }
00133 friend float realpart(const Complex& a){ return a.re; }
00134 friend float imagpart(const Complex& a){ return a.im; }
00135 friend std::ostream& operator<<(std::ostream& os, const Complex& a);
00136 friend std::istream& operator>>(std::istream& os, Complex& a);
00137 friend Complex id(const Complex& x);
00138
00139 friend Complex pow(const Complex&, const float&);
00140 friend float cos(const Complex&);
00141 friend float sin(const Complex&);
00142
00143 };
00144 inline float
00145 Complex::arg() const
00146 {
00147 if (modul2() > 1e-8)
00148 if (im > 0)
00149 return acos(re/sqrt(modul2()));
00150 else
00151 return 8*atan(1.) - acos(re/sqrt(modul2()));
00152 else
00153 return 0.;
00154
00155 }
00156 inline Complex
00157 pow(const Complex& c, const float& f)
00158 {
00159 Complex z(std::cos(f*c.arg()), std::sin(f*c.arg()));
00160 return std::pow(std::sqrt(c.modul2()),f)*z;
00161 }
00162 inline float
00163 cos(const Complex& c)
00164 {
00165 return std::cos(c.real())*std::cosh(c.imag()) + std::sin(c.real())*std::sinh(c.imag());\
00166 }
00167 inline float
00168 sin(const Complex& c)
00169 {
00170 return std::sin(c.real())*std::cosh(c.imag()) - std::cos(c.real())*std::sinh(c.imag());
00171 }
00172
00173 class cvect {public: ccreal val[2];
00174
00175 cvect() { val[0]=0.F; val[1]=0.F; }
00176 cvect(const cvect& r) { val[0]=r.val[0];val[1]=r.val[1];}
00178 cvect(ccreal r) { val[0]=r;val[1]=r;}
00179 ccreal& operator[] ( int i) { return val[i];}
00180 const ccreal& operator[] ( int i) const { return val[i];}
00181 float modul2() { return norm2(val[0])+norm2(val[1]);}
00182 ccreal operator*=(const cvect& a) { return val[0]*a.val[0] + val[1]*a.val[1]; }
00183 cvect& operator*=(const ccreal a) { val[0] *= a; val[1] *= a; return *this; }
00184 friend ccreal operator *(const cvect& a, const cvect& b) { cvect c(a); return c *= b; }
00185 friend cvect operator *(const cvect& a, const ccreal b) { cvect c(a); return c *= b; }
00186 friend cvect operator *(const ccreal b, const cvect& a) { cvect c(a); return c *= b; }
00187 cvect& operator/=(const ccreal a) {val[0] /= a;val[1] /= a; return *this; }
00188 friend cvect operator /(const cvect& a, const ccreal b) { cvect c(a); return c /= b; }
00189 cvect& operator+=(const cvect& a) { val[0] += a.val[0];val[1] += a.val[1]; return *this; }
00190 cvect& operator+=(const ccreal a) { val[0] += a;val[1] += a; return *this; }
00191 friend cvect operator +(const cvect& a, const cvect& b) { cvect c(a); return c += b; }
00192 friend cvect operator +(const cvect& a, const ccreal b) { cvect c(a); return c += b; }
00193 friend cvect operator +(const ccreal b, const cvect& a) { cvect c(a); return c += b; }
00194 cvect& operator-=(const cvect& a) { val[0] -= a.val[0];val[1] -= a.val[1]; return *this; }
00195 cvect& operator-=(const ccreal a) { val[0] -= a;val[1] -= a; return *this; }
00196 friend cvect operator -(const cvect& a, const cvect& b) { cvect c(a); return c -= b; }
00197 friend cvect operator -(const cvect& a, const ccreal b) { cvect c(a); return c -= b; }
00198 friend cvect operator -(const ccreal b, const cvect& a) { cvect c(b); return c -= a; }
00199 cvect& operator=(const cvect& a) { val[0] = a.val[0];val[1] = a.val[1]; return *this; }
00200 cvect& operator=(const ccreal a) {val[0] = a; val[1] = a; return *this; }
00201 cvect& operator-() { val[0] = -val[0]; val[1] = -val[1]; return *this; }
00202 friend float norm2( cvect& a) { return a.modul2();}
00203 friend float realpart(const cvect& a){ return realpart(a.val[0]); }
00204
00205 friend std::ostream& operator<<(std::ostream& os, cvect& a);
00206 friend std::istream& operator>>(std::istream& os, cvect& a);
00207 };
00208
00209
00210 class cmat {
00211
00212 public: ccreal val[4];
00213 cmat() {val[0]=0.F;val[1]=0.F;val[2]=0.F;val[3]=0.F; }
00214 cmat(const cmat& r) {val[0]=r.val[0];val[1]=r.val[1];val[2]=r.val[2];val[3]=r.val[3];}
00216 cmat(ccreal r) { val[0] = r;val[1] = r;val[2] = r;val[3] = r;}
00217 ccreal& operator() (int i, int j) { return val[i+i+j];}
00218 ccreal& operator[] (int i) { return val[i];}
00219 const ccreal& operator[] (int i) const { return val[i];}
00220 float modul2() {return val[0]+val[1]+val[2]+val[3];}
00221 cmat operator*=(const cmat& a)
00222 { cmat b;
00223 b.val[0] =val[0]*a.val[0]+val[1]*a.val[2];
00224 b.val[1] =val[0]*a.val[1]+val[1]*a.val[3];
00225 b.val[2] =val[2]*a.val[0]+val[3]*a.val[2];
00226 b.val[3] =val[2]*a.val[1]+val[3]*a.val[3];
00227 return b;
00228 }
00229 cmat& operator*=(const ccreal a) { val[0] *= a; val[1] *= a; val[2] *= a; val[3] *= a; return *this; }
00230 friend cmat operator *(const cmat& a, const cmat& b) { cmat c(a); return c *= b; }
00231 friend cmat operator *(const cmat& a, const ccreal b) { cmat c(a); return c *= b; }
00232 friend cmat operator *(const ccreal b, const cmat& a) { cmat c(a); return c *= b; }
00233 cmat& operator/=(const ccreal a) { val[0] /= a; val[1] /= a; val[2] /= a; val[3] /= a; return *this; }
00234 friend cmat operator /(const cmat& a, const ccreal b) { cmat c(a); return c /= b; }
00235 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; }
00236 cmat& operator+=(const ccreal a) { val[0] += a; val[1] += a; val[2] += a; val[3] += a; return *this; }
00237 friend cmat operator +(const cmat& a, const cmat& b) { cmat c(a); return c += b; }
00238 friend cmat operator +(const cmat& a, const ccreal b) { cmat c(a); return c += b; }
00239 friend cmat operator +(const ccreal b, const cmat& a) { cmat c(a); return c += b; }
00240
00241 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; }
00242 cmat& operator-=(const ccreal a) { val[0] -= a; val[2] -= a; val[1] -= a; val[3] -= a; return *this; }
00243 friend cmat operator -(const cmat& a, const cmat& b) { cmat c(a); return c -= b; }
00244 friend cmat operator -(const cmat& a, const ccreal b) { cmat c(a); return c -= b; }
00245 friend cmat operator -(const ccreal b, const cmat& a) { cmat c(b); return c -= a; }
00246 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; }
00247 cmat& operator=(const ccreal a) { val[0] = a; val[1] = a; val[2] = a; val[3] = a; return *this; }
00248 cmat& operator-() { val[0] = -val[0];val[1] = -val[1];val[2] = -val[2];val[3] = -val[3]; return *this; }
00249 friend float norm2( cmat& a) { return a.modul2();}
00250
00251 friend cvect operator *(const cmat& a, const cvect& b)
00252 { cvect c; c.val[0] = a.val[0]*b.val[0]+a.val[1]*b.val[1];
00253 c.val[1] = a.val[2]*b.val[0]+a.val[3]*b.val[1];
00254 return c ; }
00255 friend cvect operator *(const cvect& b, const cmat& a)
00256 { cvect c; c.val[0] = a.val[0]*b.val[0]+a.val[2]*b.val[1];
00257 c.val[1] = a.val[1]*b.val[0]+a.val[3]*b.val[1];
00258 return c ; }
00259 friend cvect operator /(const cvect& b,const cmat& a)
00260 { cvect c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
00261 if(norm2(s)<epsilon) { s = epsilon;}
00262 c.val[0] = deter2(b.val[0],a.val[1], b.val[1],a.val[3]) / s;
00263 c.val[1] = deter2(a.val[0],b.val[0], a.val[2],b.val[1]) / s;
00264 return c;
00265 }
00266 friend cmat operator /(const cmat& b,const cmat& a)
00267 {
00268 cmat c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
00269 if(norm2(s)<epsilon){ s = epsilon;}
00270 c.val[0] = deter2(b.val[0],a.val[1], b.val[2],a.val[3]) / s;
00271 c.val[2] = deter2(a.val[0],b.val[0], a.val[2],b.val[2]) / s;
00272 c.val[1] = deter2(b.val[1],a.val[1], b.val[3],a.val[3]) / s;
00273 c.val[3] = deter2(a.val[0],b.val[1], a.val[2],b.val[3]) / s;
00274 return c;
00275 }
00276 friend std::ostream& operator<<(std::ostream& os, cmat& a);
00277 friend std::istream& operator>>(std::istream& os, cmat& a);
00278 friend cmat id(const cvect& a);
00279 };
00280
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 class Acreal
00308 {
00309 public:
00310 long szz;
00311 creal* cc;
00312
00313 Acreal(long=0);
00314 Acreal(const Acreal&);
00315 ~Acreal() { delete [] cc;cc=0;szz = 0;}
00316 void destroy() {delete [] cc;cc=0;szz = 0;}
00317 creal& operator[] (long i) { return cc[i];}
00318 creal* operator&() { return cc;}
00319 void init(long newSize);
00320
00321
00322
00323 };
00324 class Aint
00325 {
00326 public:
00327 long szz;
00328 int* cc;
00329
00330 Aint(long=0);
00331 Aint(const Aint&);
00332 ~Aint() { delete [] cc;cc=0;szz = 0;}
00333 void destroy() {delete [] cc;cc=0;szz = 0;}
00334 int& operator[] (long i) { myassert((i< szz)&&(i>=0)); return cc[i];}
00335 int* operator&() { return cc;}
00336 void init (long);
00337 };
00338
00339
00340 class Acvect
00341 {
00342 public:
00343 long szz;
00344 cvect* cc;
00345
00346 Acvect(long=0);
00347 Acvect(const Acvect&);
00348 ~Acvect() { delete [] cc;cc=0;szz = 0;}
00349 void destroy() { delete [] cc;cc=0;szz = 0;}
00350 cvect& operator[] (long i) { return cc[i];}
00351 cvect* operator&() { return cc;}
00352 void init (long);
00353 };
00354
00355 class Acmat
00356 {
00357 public:
00358 long szz;
00359 cmat* cc;
00360
00361 Acmat(long=0);
00362 Acmat(const Acmat&);
00363 ~Acmat()
00364 {
00365 delete [] cc;
00366 cc=0;
00367 szz = 0;
00368 }
00369 void destroy() {delete [] cc;cc=0;szz = 0;}
00370 cmat& operator[] (long i) { return cc[i];}
00371 cmat* operator&() { return cc;}
00372 void init (long);
00373 };
00374
00375 class AAcmat
00376 {
00377 public:
00378 long szz;
00379 Acmat* cc;
00380
00381 AAcmat(long=0);
00382 AAcmat(const AAcmat&);
00383 ~AAcmat() { delete [] cc;cc=0;szz = 0;}
00384 void destroy() {delete [] cc;cc=0;szz = 0;}
00385 Acmat& operator[] (long i) { return cc[i];}
00386 Acmat* operator&() { return cc;}
00387 void init(long );
00388 };
00389
00390 class AAcreal
00391 {
00392 public:
00393 long szz;
00394 Acreal* cc;
00395
00396 AAcreal(long=0);
00397 AAcreal(const AAcreal& );
00398 ~AAcreal() { delete [] cc;cc=0;szz = 0;}
00399 void destroy() {delete [] cc;cc=0;szz = 0;}
00400 Acreal& operator[] (long i) { return cc[i];}
00401 Acreal* operator&() { return cc;}
00402 void init (long);
00403 };
00404 }
00405
00406 #endif
This is the FreeFEM reference manual
Provided by The KFEM project