femmisc.hpp

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

// 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 :  [e-mail]
//
// 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]);   }