FreeFem 3.5.x
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   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 */

This is the FreeFEM reference manual
Provided by The KFEM project