FreeFem 3.5.x
|
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 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; // Internal variables 00091 public: 00092 Complex() { re = 0.F; im =0.F; } // default constructor 00093 Complex(float r, float i =0.F) {re =r; im =i;} // copy constructor 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 //Complex& operator=(const float& a) { re =a; im = 0.F; return *this; } 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]; // Internal cvector to n values 00176 //---------------------------------- 00177 cvect() { val[0]=0.F; val[1]=0.F; } // constructor 00178 cvect(const cvect& r) { val[0]=r.val[0];val[1]=r.val[1];} // copy constructor 00180 cvect(ccreal r) { val[0]=r;val[1]=r;} // copy constructor from a creal 00181 ccreal& operator[] ( int i) { /*myassert((i< 2)&&(i>=0)); */ return val[i];} 00182 const ccreal& operator[] ( int i) const { /*myassert((i< 2)&&(i>=0));*/ return val[i];} 00183 float modul2() { return norm2(val[0])+norm2(val[1]);} // public fnct 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]; // Internal cvector to n values 00215 cmat() {val[0]=0.F;val[1]=0.F;val[2]=0.F;val[3]=0.F; } // constructor 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];} // copy constructor 00218 cmat(ccreal r) { val[0] = r;val[1] = r;val[2] = r;val[3] = r;} // copy constructor from a creal produces a diag mat 00219 ccreal& operator() (int i, int j) { /*myassert((i< N)&&(i>=0)&&(j<N)&&(j>=0));*/ return val[i+i+j];} 00220 ccreal& operator[] (int i) { /*myassert((i< N2)&&(i>=0));*/ return val[i];} 00221 const ccreal& operator[] (int i) const { /*myassert((i< N2)&&(i>=0));*/ return val[i];} 00222 float modul2() {return val[0]+val[1]+val[2]+val[3];} // public fnct 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 ; } //transposed 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 class Afloat 00285 { public: 00286 int szz; 00287 float* cc; 00288 00289 Afloat():szz(0),cc(NULL){} 00290 Afloat(int sz=0) {cc = 0; 00291 if(sz>0){ cc = new float[sz]; 00292 if(!cc) erreur("Out of Memory"); 00293 for(int i=0;i<sz;i++)cc[i]=0; } 00294 szz = sz; 00295 } 00296 Afloat(Afloat& a) { cc = 0; szz = 0; 00297 if(a.szz>0) { szz = a.szz; cc = new float[szz]; 00298 if(!cc) erreur("Out of Memory"); 00299 else for(int i=0;i<szz;i++)cc[i]=a.cc[i];} 00300 } 00301 ~Afloat() { delete [] cc;cc=0;szz = 0;} 00302 float& operator[] (int i) { myassert((i< szz)&&(i>=0)); return cc[i];} 00303 void init(int newSize) { myassert(!(szz || cc)); szz =newSize; 00304 cc =new float[szz]; if(!cc) erreur("Out of Memory"); 00305 for(int i=0;i<szz;i++)cc[i]=0;; 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) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];} 00320 creal* operator&() { return cc;} 00321 void init(long newSize);/* { myassert(!(szz || cc)); szz =newSize; 00322 cc =new creal[szz]; if(!cc) erreur("Out of Memory"); 00323 for(int i=0;i<szz;i++)cc[i]=0; 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) { /*myassert((i< szz)&&(i>=0));*/ 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) { /*myassert((i< szz)&&(i>=0));*/ 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) { /*myassert((i< szz)&&(i>=0));*/ 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) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];} 00403 Acreal* operator&() { return cc;} 00404 void init (long); 00405 }; 00406 } 00407 00408 #endif /* __OPCLASS_H */