Interpolation.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROC_INTERPOLATION_H
#define AMROC_INTERPOLATION_H

#include "numerical/grid_interp_extrap.h"

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

  Interpolation() : _Equations(VectorType::Length()) {}

  virtual ~Interpolation() {}

  virtual void Interpolate(GridHierarchy& GH, vec_grid_data_type& gdu, const int& nc, 
                           const point_type* xc, VectorType* uv, const DataType& ErrorValue) {

    grid_data_type phihelp(grow(gdu.bbox(),2));
    phihelp.equals(static_cast<DataType>(-1.0));
    phihelp.equals(static_cast<DataType>(1.0),gdu.bbox());

    Coords ex = phihelp.extents();
    DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
    DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
      
    DataType domain[2*dim];
    for (register int n=0; n<dim; n++) {
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }

    const DataType default_values[1] = { ErrorValue };
    DataType* values = new DataType[nc];
    register int cnt;
    for (int neq=0; neq<NEquations(); neq++) {
      grid_data_type field(grow(gdu.bbox(),2));
      equals_from(field,gdu.bbox(),gdu,gdu.bbox(),neq);
      const DataType* fields[1] = { field.data() };
      for (cnt=0; cnt<nc; cnt++) values[cnt] = ErrorValue;

      numerical::grid_interp_extrap<dim,1,DataType>( nc, values, xc[0].data(), 
                                                     default_values, ex(), domain, 
                                                     phihelp.data(), fields );  
      for (cnt = 0; cnt<nc; cnt++) 
        if (values[cnt] != ErrorValue)
          uv[cnt](neq) = values[cnt];
    }
    delete [] values;
  }

  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) {

    grid_data_type phihelp(grow(gdphi.bbox(),2));
    phihelp.equals(static_cast<DataType>(-1.0));
    phihelp.equals(gdphi,gdphi.bbox());

    Coords ex = phihelp.extents();
    DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
    DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
      
    DataType domain[2*dim];
    for (register int n=0; n<dim; n++) {
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }

    const DataType default_values[1] = { ErrorValue };
    DataType* values = new DataType[nc];
    register int cnt;
    for (int neq=0; neq<NEquations(); neq++) {
      grid_data_type field(grow(gdu.bbox(),2));
      equals_from(field,gdu.bbox(),gdu,gdu.bbox(),neq);
      const DataType* fields[1] = { field.data() };
      for (cnt=0; cnt<nc; cnt++) values[cnt] = ErrorValue;

      numerical::grid_interp_extrap<dim,1,DataType>( nc, values, xc[0].data(), 
                                                     default_values, ex(), domain, 
                                                     phihelp.data(), fields );  
      for (cnt = 0; cnt<nc; cnt++) 
        if (values[cnt] != ErrorValue)
          uv[cnt](neq) = values[cnt];
    }
    delete [] values;
  }

  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) {

    grid_data_type phihelp(grow(gdflg.bbox(),2));
    SetPhiFromBool(phihelp,gdflg);

    Coords ex = phihelp.extents();
    DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
    DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
      
    DataType domain[2*dim];
    for (register int n=0; n<dim; n++) {
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }

    const DataType default_values[1] = { ErrorValue };
    DataType* values = new DataType[nc];
    register int cnt;
    for (int neq=0; neq<NEquations(); neq++) {
      grid_data_type field(grow(gdu.bbox(),2));
      equals_from(field,gdu.bbox(),gdu,gdu.bbox(),neq);
      const DataType* fields[1] = { field.data() };
      for (cnt=0; cnt<nc; cnt++) values[cnt] = ErrorValue;

      numerical::grid_interp_extrap<dim,1,DataType>( nc, values, xc[0].data(), 
                                                     default_values, ex(), domain, 
                                                     phihelp.data(), fields );  
      for (cnt = 0; cnt<nc; cnt++) 
        if (values[cnt] != ErrorValue)
          uv[cnt](neq) = values[cnt];
    }
    delete [] values;
  }

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

    grid_data_type phihelp(grow(gd.bbox(),2));
    phihelp.equals(static_cast<DataType>(-1.0));
    phihelp.equals(static_cast<DataType>(1.0),gd.bbox());

    Coords ex = phihelp.extents();
    DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
    DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
      
    DataType domain[2*dim];
    for (register int n=0; n<dim; n++) {
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }

    const DataType default_values[1] = { ErrorValue };
    grid_data_type field(grow(gd.bbox(),2));
    field.equals(gd,gd.bbox());
    const DataType* fields[1] = { field.data() };

    numerical::grid_interp_extrap<dim,1,DataType>( nc, u, xc[0].data(), 
                                                   default_values, ex(), domain, 
                                                   phihelp.data(), fields );
  }

  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) {

    grid_data_type phihelp(grow(gdphi.bbox(),2));
    phihelp.equals(static_cast<DataType>(-1.0));
    phihelp.equals(gdphi,gdphi.bbox());

    Coords ex = phihelp.extents();
    DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
    DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
      
    DataType domain[2*dim];
    for (register int n=0; n<dim; n++) {
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }

    const DataType default_values[1] = { ErrorValue };
    grid_data_type field(grow(gd.bbox(),2));
    field.equals(gd,gd.bbox());
    const DataType* fields[1] = { field.data() };

    numerical::grid_interp_extrap<dim,1,DataType>( nc, u, xc[0].data(), 
                                                   default_values, ex(), domain, 
                                                   phihelp.data(), fields );
  }

  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) {

    grid_data_type phihelp(grow(gdflg.bbox(),2));
    SetPhiFromBool(phihelp,gdflg);

    Coords ex = phihelp.extents();
    DCoords dxh2(GH.worldStep(phihelp.stepsize())); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = GH.worldCoords(phihelp.lower(), phihelp.stepsize()) + dxh2;
    DCoords ubcorner = GH.worldCoords(phihelp.upper(), phihelp.stepsize()) + dxh2;
      
    DataType domain[2*dim];
    for (register int n=0; n<dim; n++) {
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }

    const DataType default_values[1] = { ErrorValue };
    grid_data_type field(grow(gd.bbox(),2));
    field.equals(gd,gd.bbox());
    const DataType* fields[1] = { field.data() };

    numerical::grid_interp_extrap<dim,1,DataType>( nc, u, xc[0].data(), 
                                                   default_values, ex(), domain, 
                                                   phihelp.data(), fields );
  }

  const int& NEquations() const { return _Equations; }

protected:
  void SetPhiFromBool(grid_data_type& gdphi, bool_grid_data_type& gdflg) {
    gdphi.equals(static_cast<DataType>(-1.0));
    if (dim == 1) {
      BeginFastIndex1(phi, gdphi.bbox(), gdphi.data(), DataType);
      BeginFastIndex1(flg, gdflg.bbox(), gdflg.data(), bool);
      for_1 (n, gdflg.bbox(), gdflg.bbox().stepsize())
        if (!FastIndex1(flg,n)) 
          FastIndex1(phi,n) = static_cast<DataType>(1.0);
      end_for
      EndFastIndex1(phi);
      EndFastIndex1(flg);
    }      
    else if (dim == 2) {
      BeginFastIndex2(phi, gdphi.bbox(), gdphi.data(), DataType);
      BeginFastIndex2(flg, gdflg.bbox(), gdflg.data(), bool);
      for_2 (n, m, gdflg.bbox(), gdflg.bbox().stepsize())
        if (!FastIndex2(flg,n,m)) 
          FastIndex2(phi,n,m) = static_cast<DataType>(1.0);
      end_for
      EndFastIndex2(phi);
      EndFastIndex2(flg);
    }
    else if (dim == 3) {
      BeginFastIndex3(phi, gdphi.bbox(), gdphi.data(), DataType);
      BeginFastIndex3(flg, gdflg.bbox(), gdflg.data(), bool);
      for_3 (n, m, l, gdflg.bbox(), gdflg.bbox().stepsize())
        if (!FastIndex3(flg,n,m,l)) 
          FastIndex3(phi,n,m,l) = static_cast<DataType>(1.0);
      end_for
      EndFastIndex3(phi);
      EndFastIndex3(flg);
    }
  }
    
protected:
  int _Equations;
};

#endif

---