/usr/include/fem/femMisc.hpp

// Emacs will be in -*- Mode: c++ -*-
//
// ********** DO NOT REMOVE THIS BANNER **********
//
// SUMMARY: Language for a Finite Element Method
//
//
// AUTHORS:  C. Prud'homme
// ORG    :          
// E-MAIL :  prudhomm@users.sourceforge.net
//
// ORIG-DATE:     June-94
// LAST-MOD: 12-Jul-01 at 09:50:55 by 
//
// DESCRIPTION:  
/*
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

*/
// DESCRIP-END.
//

// To change from comlex to real search everywhere for  /* C2R */

#ifndef __OPCLASS_H
#ifdef __GNUG__
#pragma interface
#endif
#define __OPCLASS_H

#include <cmath>
#include <iostream>


#define mmax(a,b)(a>b?a:b)
#define mmin(a,b)(a<b?a:b)
#define ssqr(x) ((x)*(x))


namespace fem
{
  extern int N;
  extern int N2;

  class Complex;
  class cvect;
  class cmat;

#define sizecmat sizeof(cmat)
#define sizecvect sizeof(cvect)
#define deter2(a11,a12,a21,a22) (a11*a22 - a12*a21)

  const float epsilon =(float)1.0e-20;

  //typedef float ccreal; typedef float creal;
  typedef float ccreal;
  typedef Complex creal;  // Complex : change this (see also OPClass.c) /* C2R */

  typedef Complex ccomplex;
#define sizeccomplex sizeof(ccomplex)

#define sizecreal sizeof(creal)
#define sizeccreal sizeof(ccreal)

  void cgauss( cmat& a, cvect& b);

  float norm2(const float& a);// {   return a>0?a:-a;   }
  float imagpart(const float& a);//{   return 0;   }
  float realpart(const float& a);//{   return a;   }
  cmat id(const cvect& a); 
  Complex id(const Complex& x); 
  void myassert(int what);

  class Complex

  {
  protected:
    float re,im;      // Internal variables
  public:
    Complex()  { re = 0.F; im =0.F; }        // default constructor
    Complex(float r, float i =0.F) {re =r; im =i;}  // copy constructor
  
    float& real () { return re;};
    float    real() const { return re;}
    float& imag() { return im;};
    float    imag() const { return im;}
  
    Complex conjug() {return Complex(re,-im);}
   
    Complex& operator=(const Complex& a) { re =a.re; im = a.im; return *this; }   
    //Complex& operator=(const float& a) { re =a; im = 0.F; return *this; }   
    Complex& operator*=(const Complex& a) { re =a.re*re - a.im*im;im= a.re*im + a.im*re; return *this; }   
    Complex& operator/=(const Complex& a) { 
      float xl=ssqr(a.re) + ssqr(a.im); re =(a.re*re+a.im*im)/xl; im =(a.re*im-a.im*re)/xl; 
      return *this; 
    }   
    Complex& operator+=(const Complex& a) { re +=a.re; im += a.im; return *this; }   
    Complex& operator-=(const Complex& a) { re -=a.re; im -= a.im; return *this; }   
    Complex& operator*=(const float a) { re *= a;im *= a; return *this; }   
    Complex& operator/=(const float a) { re /= a;im /= a; return *this; }   
    Complex& operator+=(const float a) { re +=a;  return *this; }   
    Complex& operator-=(const float a) { re -=a;  return *this; }   
    Complex& operator-() { re =-re; im =-im; return *this; }   
    float modul2() const {return ssqr(re)+ssqr(im);} 
    float arg() const;
  
    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);}
    friend Complex operator *(const Complex a, const float b) { return Complex(a.re*b, a.im*b); }  
    friend Complex operator *(const float b, const Complex a) { return Complex(a.re*b, a.im*b); }   
    friend Complex operator /(const Complex a, const float b) { return Complex(a.re/b, a.im/b); }  
    friend Complex operator /(const float b, const Complex a) { return Complex(b,0.F)/a; }   
    friend Complex operator /(const Complex a, const Complex b) { 
      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; 
      return c;
    }
    friend Complex operator +(const Complex a, const Complex b)  {return Complex(a.re+b.re, a.im+b.im);}
    friend Complex operator +(const Complex a, const float b) { return Complex(a.re+b, a.im); }  
    friend Complex operator +(const float b, const Complex a) { return Complex(a.re+b, a.im); }   
    friend Complex operator -(const Complex a, const Complex b)  {return Complex(a.re-b.re, a.im-b.im);}
    friend Complex operator -(const Complex a, const float b) { return Complex(a.re-b, a.im); }  
    friend Complex operator -(const float b, const Complex a) { return Complex(b-a.re, -a.im); }  
    friend float norm2( const Complex& a) { return a.modul2();   }
    friend float realpart(const Complex& a){   return a.re;   }
    friend float imagpart(const Complex& a){   return a.im;   }
    friend std::ostream& operator<<(std::ostream& os, const Complex& a);
    friend std::istream& operator>>(std::istream& os,  Complex& a);
    friend Complex id(const Complex& x); 
  
    friend Complex pow(const Complex&, const float&);
    friend float cos(const Complex&);
    friend float sin(const Complex&);
  
  };
  inline float
  Complex::arg() const
  {
    if (modul2() > 1e-8)
      if (im > 0)
	return acos(re/sqrt(modul2()));
      else
	return 8*atan(1.) - acos(re/sqrt(modul2()));
    else
      return 0.;
  
  }
  inline Complex
  pow(const Complex& c, const float& f)
  {
    Complex z(std::cos(f*c.arg()), std::sin(f*c.arg())); 
    return std::pow(std::sqrt(c.modul2()),f)*z;
  }
  inline float
  cos(const Complex& c)
  {
      return std::cos(c.real())*std::cosh(c.imag()) + std::sin(c.real())*std::sinh(c.imag());\
  }
  inline float
  sin(const Complex& c)
  {
    return std::sin(c.real())*std::cosh(c.imag()) - std::cos(c.real())*std::sinh(c.imag());
  }
  //____________________________________________________________________________
  class cvect {public:  ccreal  val[2];       //  Internal cvector to n values
    //----------------------------------
    cvect()  {  val[0]=0.F; val[1]=0.F; }         // constructor
    cvect(const cvect& r) { val[0]=r.val[0];val[1]=r.val[1];} // copy constructor
    ///cvect( cvect& r) {  val[0]=r.val[0];val[1]=r.val[1];} // copy constructor
    cvect(ccreal r) {   val[0]=r;val[1]=r;} // copy constructor from a creal
    ccreal& operator[] ( int i)   { /*myassert((i< 2)&&(i>=0)); */ return val[i];}
    const ccreal& operator[] ( int i) const   { /*myassert((i< 2)&&(i>=0));*/  return val[i];}
    float modul2()   { return norm2(val[0])+norm2(val[1]);} // public fnct
    ccreal operator*=(const cvect& a)   { return val[0]*a.val[0] + val[1]*a.val[1]; }   
    cvect& operator*=(const ccreal a)   { val[0] *= a; val[1] *= a; return *this; }   
    friend ccreal operator *(const cvect& a, const cvect& b) { cvect c(a); return c *= b; }  
    friend cvect operator *(const cvect& a, const ccreal b) { cvect c(a); return c *= b; }  
    friend cvect operator *(const ccreal b, const cvect& a) {  cvect c(a); return c *= b; }   
    cvect& operator/=(const ccreal a)   {val[0] /= a;val[1] /= a; return *this; }   
    friend cvect operator /(const cvect& a, const ccreal b) { cvect c(a); return c /= b; }  
    cvect& operator+=(const cvect& a)   { val[0] += a.val[0];val[1] += a.val[1]; return *this; }   
    cvect& operator+=(const ccreal a)   { val[0] += a;val[1] += a; return *this; }   
    friend cvect operator +(const cvect& a, const cvect& b) { cvect c(a); return c += b; }  
    friend cvect operator +(const cvect& a, const ccreal b) { cvect c(a); return c += b; }  
    friend cvect operator +(const ccreal b, const cvect& a) {  cvect c(a); return c += b; }   
    cvect& operator-=(const cvect& a)  { val[0] -= a.val[0];val[1] -= a.val[1]; return *this; }   
    cvect& operator-=(const ccreal a)  { val[0] -= a;val[1] -= a; return *this; }   
    friend cvect operator -(const cvect& a, const cvect& b) { cvect c(a); return c -= b; }  
    friend cvect operator -(const cvect& a, const ccreal b) { cvect c(a); return c -= b; }  
    friend cvect operator -(const ccreal b, const cvect& a) {  cvect c(b); return c -= a; }   
    cvect& operator=(const cvect& a)   { val[0] = a.val[0];val[1] = a.val[1]; return *this; }     
    cvect& operator=(const ccreal a)   {val[0] = a; val[1] = a; return *this; }     
    cvect& operator-()   { val[0] = -val[0]; val[1] = -val[1]; return *this; }     
    friend float norm2( cvect& a) { return a.modul2();}
    friend float realpart(const cvect& a){   return realpart(a.val[0]);   }

    friend std::ostream& operator<<(std::ostream& os,  cvect& a);
    friend std::istream& operator>>(std::istream& os,  cvect& a);
  };

  //_______________________________
  class cmat { 
    //---------
  public: ccreal  val[4];       //  Internal cvector to n values
    cmat()  {val[0]=0.F;val[1]=0.F;val[2]=0.F;val[3]=0.F; }         // constructor
    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
    ///cmat(cmat& r) { val[0]=r.val[0];  val[1]=r.val[1];  val[2]=r.val[2];  val[3]=r.val[3]; }  // copy constructor
    cmat(ccreal r) {  val[0] = r;val[1] = r;val[2] = r;val[3] = r;} // copy constructor from a creal produces a diag mat
    ccreal& operator() (int i, int j)   { /*myassert((i< N)&&(i>=0)&&(j<N)&&(j>=0));*/ return val[i+i+j];}
    ccreal& operator[] (int i)    { /*myassert((i< N2)&&(i>=0));*/  return val[i];}
    const ccreal& operator[] (int i) const    { /*myassert((i< N2)&&(i>=0));*/  return val[i];}
    float modul2()    {return val[0]+val[1]+val[2]+val[3];} // public fnct
    cmat operator*=(const cmat& a)  
    { cmat b;  
    b.val[0] =val[0]*a.val[0]+val[1]*a.val[2];
    b.val[1] =val[0]*a.val[1]+val[1]*a.val[3];
    b.val[2] =val[2]*a.val[0]+val[3]*a.val[2];
    b.val[3] =val[2]*a.val[1]+val[3]*a.val[3];
    return b; 
    }   
    cmat& operator*=(const ccreal a) { val[0] *= a; val[1] *= a; val[2] *= a; val[3] *= a; return *this; }   
    friend cmat operator *(const cmat& a, const cmat& b) { cmat c(a); return c *= b; }  
    friend cmat operator *(const cmat& a, const ccreal b) { cmat c(a); return c *= b; }  
    friend cmat operator *(const ccreal b, const cmat& a) { cmat c(a); return c *= b; }   
    cmat& operator/=(const ccreal a) { val[0] /= a; val[1] /= a; val[2] /= a; val[3] /= a; return *this; }   
    friend cmat operator /(const cmat& a, const ccreal b) { cmat c(a); return c /= b; }  
    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; }   
    cmat& operator+=(const ccreal a) { val[0] += a; val[1] += a; val[2] += a; val[3] += a; return *this; }   
    friend cmat operator +(const cmat& a, const cmat& b) { cmat c(a); return c += b; }  
    friend cmat operator +(const cmat& a, const ccreal b) { cmat c(a); return c += b; }  
    friend cmat operator +(const ccreal b, const cmat& a) {  cmat c(a); return c += b; }   

    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; }   
    cmat& operator-=(const ccreal a)  { val[0] -= a; val[2] -= a; val[1] -= a; val[3] -= a; return *this; }   
    friend cmat operator -(const cmat& a, const cmat& b) { cmat c(a); return c -= b; }  
    friend cmat operator -(const cmat& a, const ccreal b) { cmat c(a); return c -= b; }  
    friend cmat operator -(const ccreal b, const cmat& a) {  cmat c(b); return c -= a; }   
    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; }     
    cmat& operator=(const ccreal a)   { val[0] = a; val[1] = a; val[2] = a; val[3] = a; return *this; }     
    cmat& operator-()   { val[0] = -val[0];val[1] = -val[1];val[2] = -val[2];val[3] = -val[3]; return *this; }     
    friend float norm2( cmat& a) { return a.modul2();}
    /**/
    friend cvect operator *(const cmat& a, const cvect& b) 
    {   cvect c;  c.val[0] = a.val[0]*b.val[0]+a.val[1]*b.val[1];
    c.val[1] = a.val[2]*b.val[0]+a.val[3]*b.val[1]; 
    return c ; }  
    friend cvect operator *(const cvect& b, const cmat& a)
    { cvect c;    c.val[0] = a.val[0]*b.val[0]+a.val[2]*b.val[1];
    c.val[1] = a.val[1]*b.val[0]+a.val[3]*b.val[1]; 
    return c ; }  //transposed
    friend cvect operator /(const  cvect& b,const cmat& a) 
    {       cvect c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
    if(norm2(s)<epsilon) { s = epsilon;}
    c.val[0] =  deter2(b.val[0],a.val[1], b.val[1],a.val[3]) / s;
    c.val[1] =  deter2(a.val[0],b.val[0], a.val[2],b.val[1]) / s;
    return c;
    }  
    friend cmat operator /(const cmat& b,const cmat& a) 
    {
      cmat c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
      if(norm2(s)<epsilon){ s = epsilon;}
      c.val[0] =  deter2(b.val[0],a.val[1], b.val[2],a.val[3]) / s;
      c.val[2] =  deter2(a.val[0],b.val[0], a.val[2],b.val[2]) / s;
      c.val[1] =  deter2(b.val[1],a.val[1], b.val[3],a.val[3]) / s;
      c.val[3] =  deter2(a.val[0],b.val[1], a.val[2],b.val[3]) / s;
      return c;
    } 
    friend std::ostream& operator<<(std::ostream& os,  cmat& a);
    friend std::istream& operator>>(std::istream& os,  cmat& a);
    friend cmat id(const cvect& a); 
  };

  /*
    class Afloat
    { public:     
    int szz; 
    float* cc; 
  
    Afloat():szz(0),cc(NULL){}
    Afloat(int sz=0)  {cc = 0;   
    if(sz>0){ cc = new float[sz]; 
    if(!cc) erreur("Out of Memory");
    for(int i=0;i<sz;i++)cc[i]=0; } 
    szz = sz;
    }
    Afloat(Afloat& a)       {   cc = 0;  szz = 0; 
    if(a.szz>0) { szz = a.szz; cc = new float[szz]; 
    if(!cc) erreur("Out of Memory");
    else for(int i=0;i<szz;i++)cc[i]=a.cc[i];} 
    }
    ~Afloat()           { delete  [] cc;cc=0;szz = 0;} 
    float& operator[] (int i)   { myassert((i< szz)&&(i>=0)); return cc[i];}
    void init(int newSize)    { myassert(!(szz || cc)); szz =newSize; 
    cc =new float[szz]; if(!cc) erreur("Out of Memory");
    for(int i=0;i<szz;i++)cc[i]=0;;   
    } 
    };
  */
  class Acreal
  { 
  public:   
    long szz; 
    creal* cc; 
  
    Acreal(long=0);
    Acreal(const Acreal&);
    ~Acreal()  { delete  [] cc;cc=0;szz = 0;} 
    void destroy()  {delete  [] cc;cc=0;szz = 0;} 
    creal& operator[] (long i)  { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
    creal* operator&()      { return cc;} 
    void init(long newSize);/*    { myassert(!(szz || cc)); szz =newSize; 
				  cc =new creal[szz];   if(!cc) erreur("Out of Memory");
				  for(int i=0;i<szz;i++)cc[i]=0; 
				  } */
  };
  class Aint
  { 
  public:
    long szz; 
    int* cc; 
  
    Aint(long=0);
    Aint(const Aint&);
    ~Aint()           { delete  [] cc;cc=0;szz = 0;}
    void destroy()  {delete  [] cc;cc=0;szz = 0;} 
    int& operator[] (long i)  { myassert((i< szz)&&(i>=0)); return cc[i];}
    int* operator&()      { return cc;} 
    void init (long);
  };


  class Acvect
  {
  public:
    long szz; 
    cvect* cc; 
  
    Acvect(long=0);
    Acvect(const Acvect&);
    ~Acvect()           { delete  [] cc;cc=0;szz = 0;} 
    void destroy()  { delete  [] cc;cc=0;szz = 0;} 
    cvect& operator[] (long i)  { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
    cvect* operator&()      { return cc;} 
    void init (long);
  };

  class Acmat
  { 
  public: 
    long szz; 
    cmat* cc; 
  
    Acmat(long=0);
    Acmat(const Acmat&);
    ~Acmat()          
    {
      delete  [] cc;
      cc=0;
      szz = 0;
    } 
    void destroy()  {delete  [] cc;cc=0;szz = 0;} 
    cmat& operator[] (long i)   { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
    cmat* operator&()     { return cc;} 
    void init (long);
  };

  class AAcmat
  {
  public:
    long szz; 
    Acmat* cc; 

    AAcmat(long=0);
    AAcmat(const AAcmat&);
    ~AAcmat()           { delete  [] cc;cc=0;szz = 0;} 
    void destroy()  {delete  [] cc;cc=0;szz = 0;} 
    Acmat& operator[] (long i)  { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
    Acmat* operator&()      { return cc;} 
    void init(long );
  };

  class AAcreal
  { 
  public:
    long szz; 
    Acreal* cc; 
  
    AAcreal(long=0);
    AAcreal(const AAcreal& );
    ~AAcreal()          { delete  [] cc;cc=0;szz = 0;} 
    void destroy()  {delete  [] cc;cc=0;szz = 0;} 
    Acreal& operator[] (long i)   { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
    Acreal* operator&()     { return cc;} 
    void init (long);
  };
}

#endif /* __OPCLASS_H */
libfreefem-dev 3.5.8-5.1+b2 / usr / include / fem / femMisc.hpp
