Main Page | Alphabetical List | Class List | File List | Class Members | File Members

femMisc.hpp

00001 // Emacs will be in -*- Mode: c++ -*- 00002 // 00003 // ********** DO NOT REMOVE THIS BANNER ********** 00004 // 00005 // SUMMARY: Language for a Finite Element Method 00006 // 00007 // 00008 // AUTHORS: C. Prud'homme 00009 // ORG : 00010 // E-MAIL : prudhomm@users.sourceforge.net 00011 // 00012 // ORIG-DATE: June-94 00013 // LAST-MOD: 12-Jul-01 at 09:50:55 by 00014 // 00015 // DESCRIPTION: 00016 /* 00017 This program is free software; you can redistribute it and/or modify 00018 it under the terms of the GNU General Public License as published by 00019 the Free Software Foundation; either version 2 of the License, or 00020 (at your option) any later version. 00021 00022 This program is distributed in the hope that it will be useful, 00023 but WITHOUT ANY WARRANTY; without even the implied warranty of 00024 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00025 GNU General Public License for more details. 00026 00027 You should have received a copy of the GNU General Public License 00028 along with this program; if not, write to the Free Software 00029 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 00030 00031 */ 00032 // DESCRIP-END. 00033 // 00034 00035 // To change from comlex to real search everywhere for /* C2R */ 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 //typedef float ccreal; typedef float creal; 00068 typedef float ccreal; 00069 typedef Complex creal; // Complex : change this (see also OPClass.c) /* C2R */ 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);// { return a>0?a:-a; } 00080 float imagpart(const float& a);//{ return 0; } 00081 float realpart(const float& a);//{ return a; } 00082 void myassert(int what); 00083 00084 class Complex 00085 00086 { 00087 protected: 00088 float re,im; // Internal variables 00089 public: 00090 Complex() { re = 0.F; im =0.F; } // default constructor 00091 Complex(float r, float i =0.F) {re =r; im =i;} // copy constructor 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 //Complex& operator=(const float& a) { re =a; im = 0.F; return *this; } 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]; // Internal cvector to n values 00174 //---------------------------------- 00175 cvect() { val[0]=0.F; val[1]=0.F; } // constructor 00176 cvect(const cvect& r) { val[0]=r.val[0];val[1]=r.val[1];} // copy constructor 00178 cvect(ccreal r) { val[0]=r;val[1]=r;} // copy constructor from a creal 00179 ccreal& operator[] ( int i) { /*myassert((i< 2)&&(i>=0)); */ return val[i];} 00180 const ccreal& operator[] ( int i) const { /*myassert((i< 2)&&(i>=0));*/ return val[i];} 00181 float modul2() { return norm2(val[0])+norm2(val[1]);} // public fnct 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]; // Internal cvector to n values 00213 cmat() {val[0]=0.F;val[1]=0.F;val[2]=0.F;val[3]=0.F; } // constructor 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];} // copy constructor 00216 cmat(ccreal r) { val[0] = r;val[1] = r;val[2] = r;val[3] = r;} // copy constructor from a creal produces a diag mat 00217 ccreal& operator() (int i, int j) { /*myassert((i< N)&&(i>=0)&&(j<N)&&(j>=0));*/ return val[i+i+j];} 00218 ccreal& operator[] (int i) { /*myassert((i< N2)&&(i>=0));*/ return val[i];} 00219 const ccreal& operator[] (int i) const { /*myassert((i< N2)&&(i>=0));*/ return val[i];} 00220 float modul2() {return val[0]+val[1]+val[2]+val[3];} // public fnct 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 ; } //transposed 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 class Afloat 00283 { public: 00284 int szz; 00285 float* cc; 00286 00287 Afloat():szz(0),cc(NULL){} 00288 Afloat(int sz=0) {cc = 0; 00289 if(sz>0){ cc = new float[sz]; 00290 if(!cc) erreur("Out of Memory"); 00291 for(int i=0;i<sz;i++)cc[i]=0; } 00292 szz = sz; 00293 } 00294 Afloat(Afloat& a) { cc = 0; szz = 0; 00295 if(a.szz>0) { szz = a.szz; cc = new float[szz]; 00296 if(!cc) erreur("Out of Memory"); 00297 else for(int i=0;i<szz;i++)cc[i]=a.cc[i];} 00298 } 00299 ~Afloat() { delete [] cc;cc=0;szz = 0;} 00300 float& operator[] (int i) { myassert((i< szz)&&(i>=0)); return cc[i];} 00301 void init(int newSize) { myassert(!(szz || cc)); szz =newSize; 00302 cc =new float[szz]; if(!cc) erreur("Out of Memory"); 00303 for(int i=0;i<szz;i++)cc[i]=0;; 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) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];} 00318 creal* operator&() { return cc;} 00319 void init(long newSize);/* { myassert(!(szz || cc)); szz =newSize; 00320 cc =new creal[szz]; if(!cc) erreur("Out of Memory"); 00321 for(int i=0;i<szz;i++)cc[i]=0; 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) { /*myassert((i< szz)&&(i>=0));*/ 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) { /*myassert((i< szz)&&(i>=0));*/ 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) { /*myassert((i< szz)&&(i>=0));*/ 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) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];} 00401 Acreal* operator&() { return cc;} 00402 void init (long); 00403 }; 00404 } 00405 00406 #endif /* __OPCLASS_H */

This is the FreeFEM reference manual
Provided by The KFEM project