Main Page   Alphabetical List   Compound List   File List   Compound 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