F77Interpolation.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_F77_INTERPOLATION_H
#define AMROC_F77_INTERPOLATION_H

#include "Interpolation.h"

template <class VectorType, int dim>
class F77Interpolation : public Interpolation<VectorType,dim> {
  typedef Interpolation<VectorType,dim> base;
  typedef typename VectorType::InternalDataType DataType;
public:
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::grid_data_type grid_data_type;  
  typedef typename base::point_type point_type;
  typedef typename base::bool_grid_data_type bool_grid_data_type;

  typedef generic_fortran_func generic_func_type;

  typedef void (*int_1_func_type) ( const DOUBLE q[], const INTEGER &mx, 
                                    const INTEGER &meqn, const INTEGER &nc, 
                                    DOUBLE qint[], const DOUBLE xc[],
                                    const DOUBLE corn[], const DOUBLE dx[], 
                                    const DOUBLE& ev);

  typedef void (*int_2_func_type) ( const DOUBLE q[], const INTEGER &mx, 
                                    const INTEGER &my, 
                                    const INTEGER &meqn, const INTEGER &nc, 
                                    DOUBLE qint[], const DOUBLE xc[],
                                    const DOUBLE corn[], const DOUBLE dx[], 
                                    const DOUBLE& ev);

  typedef void (*int_3_func_type) ( const DOUBLE q[], const INTEGER &mx, 
                                    const INTEGER &my, const INTEGER &mz, 
                                    const INTEGER &meqn, const INTEGER &nc, 
                                    DOUBLE qint[], const DOUBLE xc[],
                                    const DOUBLE corn[], const DOUBLE dx[], 
                                    const DOUBLE& ev);

  F77Interpolation() : base(), f_intpol(0) {}
  F77Interpolation(generic_func_type intpol) : base(), f_intpol(intpol){}
  
  virtual ~F77Interpolation() {}
  
  virtual void Interpolate(GridHierarchy& GH, vec_grid_data_type& gdu, const int& nc, 
                           const point_type* xc, VectorType* uv, const DataType& ErrorValue) {
    if (!f_intpol) return;
    DCoords dx = GH.worldStep(gdu.stepsize());
    DCoords corn = GH.worldCoords(gdu.lower(),gdu.stepsize());
    int mx[3];
    for (register int d=0; d<dim; d++) 
      mx[d] = gdu.extents(d);

    if (dim == 1)
      ((int_1_func_type) f_intpol)((DataType*)(gdu.data()),AA(1,mx),base::NEquations(),nc,
                                   uv[0].data(),xc[0].data(),corn(),dx(),ErrorValue); 
    else if (dim == 2) 
      ((int_2_func_type) f_intpol)((DataType*)(gdu.data()),AA(2,mx),base::NEquations(),nc,
                                   uv[0].data(),xc[0].data(),corn(),dx(),ErrorValue); 
    else if (dim == 3) 
      ((int_3_func_type) f_intpol)((DataType*)(gdu.data()),AA(3,mx),base::NEquations(),nc,
                                   uv[0].data(),xc[0].data(),corn(),dx(),ErrorValue); 
  }    

  virtual void Interpolate(GridHierarchy& GH, vec_grid_data_type& gdu, 
                           grid_data_type& gdphi, const int& nc, const point_type* xc, 
                           VectorType* uv, const DataType& ErrorValue) {
    base::Interpolate(GH,gdu,gdphi,nc,xc,uv,ErrorValue);
  }

  virtual void Interpolate(GridHierarchy& GH, vec_grid_data_type& gdu, 
                           bool_grid_data_type& gdflg, const int& nc, const point_type* xc, 
                           VectorType* uv, const DataType& ErrorValue) {
    base::Interpolate(GH,gdu,gdflg,nc,xc,uv,ErrorValue);
  }

  virtual void Interpolate(GridHierarchy& GH, grid_data_type& gd, const int& nc, 
                           const point_type* xc, DataType* u, const DataType& ErrorValue) {
    if (!f_intpol) return;
    DCoords dx = GH.worldStep(gd.stepsize());
    DCoords corn = GH.worldCoords(gd.lower(),gd.stepsize());
    int mx[3];
    for (register int d=0; d<dim; d++) 
      mx[d] = gd.extents(d);
    int meqn = 1;

    if (dim == 1) 
      ((int_1_func_type) f_intpol)(gd.data(),AA(1,mx),meqn,nc,
                                   u,xc[0].data(),corn(),dx(),ErrorValue); 
    else if (dim == 2) 
      ((int_2_func_type) f_intpol)(gd.data(),AA(2,mx),meqn,nc,
                                   u,xc[0].data(),corn(),dx(),ErrorValue); 
    else if (dim == 3) 
      ((int_3_func_type) f_intpol)(gd.data(),AA(3,mx),meqn,nc,
                                   u,xc[0].data(),corn(),dx(),ErrorValue); 
  }

  virtual void Interpolate(GridHierarchy& GH, grid_data_type& gd, 
                           grid_data_type& gdphi, const int& nc, const point_type* xc, 
                           DataType* u, const DataType& ErrorValue) {
    base::Interpolate(GH,gd,gdphi,nc,xc,u,ErrorValue);
  }

  virtual void Interpolate(GridHierarchy& GH, grid_data_type& gd, 
                           bool_grid_data_type& gdflg, const int& nc, const point_type* xc, 
                           DataType* u, const DataType& ErrorValue) {
    base::Interpolate(GH,gd,gdflg,nc,xc,u,ErrorValue);
  }

  inline void SetIntPolFunc(generic_func_type intpol) { f_intpol = intpol; }
  generic_func_type GetIntPolFunc() const { return f_intpol; }

protected:
  generic_func_type f_intpol;
};


#endif

---