Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Function.hpp

//
// Beschreibung von Funktionen
//

#ifndef _FUNCTION_
# define _FUNCTION_

#include <math.h>
#include "Math.hpp"
#include <iostream>
using namespace std;

struct Function {
  virtual double Eval         (double *x)=0 ; 
  virtual void   SetParameters(double *p)=0 ;


  double GaussQuadrature(double *XA, double *XE, int D, int order=4) ;

private:
  double GQ(double *x, int i, double *XA, double *XE, int D , int order) ;
  // definite integral of *this over (XA[i]..XE[i])x..x(XA[D-1]..XE[D-1]) ;
  // where the argument of F consists of (x[0],...,x[i-1]  ,  y[i],...,y[D-1]) ;
  // and y[k] \in XA[k]..XE[k], k>=i 
} ;


struct ConstFunction : public Function {
  double  par ;
  void   SetParameters(double *p) {par=*p     ;}
  double Eval(double *)           {return par ;} 
} ;


template<int N>
struct SparseGridFunction : public Function {
  int    MaxLevel[N] ;
         SparseGridFunction() ;
  void   SetParameters(double *p) ;
  double Eval(double *x) ;
} ;

template<int N>
struct FullGridFunction : public SparseGridFunction<N> {
  double Eval(double *x) ;
} ;


template<int N>
struct TensorNPolynomFunction : public Function {
  int    p[N] ;
  double c[N][100] ;
  double q         ;
  void   SetParameters(double *p) ;
  double Eval(double *x) ;
} ;

struct SubFunction : public Function {
  Function *F ;
  int dir,dim ;
  double xdir ;
  void   SetParameters(double *p) ;
  double Eval         (double *x) ;
} ;

template <int DIM>
struct GeneralRadialFunction : public Function {
  double a,e,y[DIM] ;
  void   SetParameters(double *p) ; // p[0]=a, p[1]=, p[2]..p[DIM+1]=y[0]..y[DIM-1] 
  double Eval         (double *x) ;
} ;

template <int DIM>
struct GeneralRadialFunctionDX : public GeneralRadialFunction<DIM> {
  double Eval         (double *x) ;
} ;
template <int DIM>
struct GeneralRadialFunctionDY : public GeneralRadialFunction<DIM> {
  double Eval         (double *x) ;
} ;

template <int DIM>
struct LaplaceGeneralRadialFunction : public GeneralRadialFunction<DIM> {
  double Eval         (double *x) ;
} ;


//
//  implementations
//

// Sparse/FullGridFunction 
template<int N>
SparseGridFunction<N>::SparseGridFunction() {
 for (int i=0; i<N; i++) MaxLevel[i]=MAXINT ;
}
template<int N>
void SparseGridFunction<N>::SetParameters(double *p) {
 for (int i=0; i<N; i++) MaxLevel[i]=(int)(p[i]+1e-9) ;
}
template<int N>
double SparseGridFunction<N>::Eval(double *x) {
 double s=0 ;
 for (int i=0; i<N; i++) {
   s += x[i]; 
   if (x[i]>MaxLevel[i]) return MAXINT ;
 }
 return s ;
}


template<int N>
double FullGridFunction<N>::Eval(double *x) {
 double s=0 ;
 for (int i=0; i<N; i++) {
   s = s>fabs(x[i]) ? s : fabs(x[i]) ; 
   if (x[i]>MaxLevel[i]) return MAXINT ;
 }
 return s ;
}


template<int N>
void TensorNPolynomFunction<N>::SetParameters(double *par) {
 int j=0,n,i ;
 q=par[j++] ; 
 for (n=0; n<N; n++) p[n]=(int)par[j++] ;
 for (n=0; n<N; n++) 
   for (i=0; i<=p[n]; i++) c[n][i]=par[j++] ;
}


template<int N>
double TensorNPolynomFunction<N>::Eval(double *x) {
 int n,i ;
 double w=1,ww  ;
 for (n=0; n<N; n++) {
   ww=0 ;
   for (i=0; i<=p[n]; i++) ww +=c[n][i]*pow(x[n],i) ;
   w *= ww ;
  }
 return q*w ;
}

void SubFunction::SetParameters(double *p) {
  if (sizeof(double)<sizeof(Function *)) {
    std::cout<<"SubFunction::SetParameters: works only, if sizeof(double) >= sizeof(Function *)\n" ;
    exit(-1) ;
  }
  F   =*((Function **)p) ;
  dim = (int)(p[1]+1e-10);
  dir = (int)(p[2]+1e-10);
  xdir=      p[3];
}

double SubFunction::Eval(double *x) {
  double y[100] ;
  int q=0,i ;
  for (i=0; i<dim; i++) {
    if (i==dir) continue ;
    y[i]=x[q++] ;
  }
  y[dir]=xdir ;
  return F->Eval(y) ;
}


template <int DIM>
void GeneralRadialFunction<DIM>::SetParameters(double *p) {
  a=p[0] ; e=p[1] ;
  for (int i=0; i<DIM; i++) y[i]=p[2+i] ; 
}

template <int DIM>
double GeneralRadialFunction<DIM>::Eval(double *x) {
  double r=e ;
  for (int i=0; i<DIM; i++) r += sqr(x[i]-y[i]) ;
  return pow(r,a/2) ;
}

template <int DIM>
double GeneralRadialFunctionDX<DIM>::Eval(double *x) {
  double r=e ;
  for (int i=0; i<DIM; i++) r += sqr(x[i]-y[i]) ;
  return a*pow(r,a/2-1)*(x[0]-y[0]) ;
}

template <int DIM>
double GeneralRadialFunctionDY<DIM>::Eval(double *x) {
  double r=e ;
  for (int i=0; i<DIM; i++) r += sqr(x[i]-y[i]) ;
  return a*pow(r,a/2-1)*(x[1]-y[1]) ;
}

template <int DIM>
double LaplaceGeneralRadialFunction<DIM>::Eval(double *x) {
  double r=e ;
  for (int i=0; i<DIM; i++) r += sqr(x[i]-y[i]) ;
  return 2*a*pow(r,a/2-1)*( a/2 - (a/2-1)*e/r ) ;
}

#endif

---