AMRInterpolation.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_AMRINTERPOLATION_H
#define AMROC_AMRINTERPOLATION_H

#include "Geom.h"
#include "Interpolation.h"

template <class VectorType, int dim>
class AMRInterpolation : public Geom<typename VectorType::InternalDataType,dim>, 
                         public AMRBase<VectorType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  typedef AMRBase<VectorType,dim> base;
  typedef Geom<DataType,dim> geom_base;

public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::vec_grid_data_type vec_grid_data_type;  
  typedef GridFunction<DataType,dim> grid_fct_type;   
  typedef GridData<DataType,dim> grid_data_type;   
  typedef Vector<DataType,dim> point_type;
  typedef Interpolation<VectorType,dim> interpolation_type;

  AMRInterpolation() : base(), _ErrorValue(-1.e37) {
    _interpolation = new interpolation_type(); 
  }

  virtual ~AMRInterpolation() {
    if (_interpolation) delete _interpolation;
  }

  void PointsValues(vec_grid_fct_type& u, const double& t, const int& Npoints, 
                    const point_type* xc, VectorType* uv, const point_type* xct=0) {

    if (Npoints<=0 || !uv) return;
    for (register int n=0; n<Npoints; n++) 
      uv[n] = _ErrorValue;

    for (register int l=0; l<=FineLevel(base::GH()); l++) {
      int Time = CurrentTime(base::GH(),l);
      PointsValues(u,Time,l,t,Npoints,xc,uv,xct);
    }
  }

  void PointsValues(vec_grid_fct_type& u, const int& Time, const int& Level, const double& t,
                    const int& Npoints, const point_type* xc, VectorType* uv, const point_type* xct=0) {

    if (Npoints<=0 || !uv) return;
    bool* lused = new bool[Npoints];

    forall (u,Time,Level,c)
      BBox bb = u.interiorbbox(Time,Level,c);
      int ncl;
      if (xct) 
        ncl = geom_base::LocalList(base::GH(),bb,Npoints,xct,lused);
      else 
        ncl = geom_base::LocalList(base::GH(),bb,Npoints,xc,lused);
      if (ncl>0) {
        register int n, nl;
        point_type* xcl = new point_type[ncl];
        VectorType* uvl = new VectorType[ncl];
        for (n=0, nl=0; n<Npoints; n++) 
          if (lused[n]) xcl[nl++] = xc[n];

        GetGrid(u(Time,Level,c),Time,Level,c,t,ncl,xcl,uvl);
          
        for (n=0, nl=0; n<Npoints; n++) 
          if (lused[n]) uv[n] = uvl[nl++];

        delete [] xcl;
        delete [] uvl;
      }
    end_forall

    delete [] lused;
  }    

  int LocalPoints(vec_grid_fct_type& u, const int& Time, const int& Level, 
                  const int& Npoints, const point_type* xc, bool* lused) {
    int ncl = 0;
    if (Npoints<=0 || !xc) return ncl;

    for (register int n=0; n<Npoints; n++) 
      lused[n] = false;

    bool* lused_help = new bool[Npoints];
    forall (u,Time,Level,c)
      BBox bb = u.interiorbbox(Time,Level,c);
      int ncl_help = geom_base::LocalList(base::GH(),bb,Npoints,xc,lused_help);
      if (ncl_help>0)
        for (register int n=0; n<Npoints; n++) 
          if (lused_help[n]) lused[n] = true;
    end_forall
    delete [] lused_help;

    for (register int n=0; n<Npoints; n++) 
      if (lused[n]) ncl++;
    return ncl;
  }    

  virtual void GetGrid(vec_grid_data_type& gdu, const int& Time, const int& Level, 
                       const int& c, const double& t, const int& Npoints, 
                       const point_type* xcl, VectorType* uvl) {
    if (_interpolation)
      Interpolation_().Interpolate(base::GH(),gdu,Npoints,xcl,uvl,_ErrorValue);
  }

  void PointsValuesPar(vec_grid_fct_type& u, const double& t, const int& Npoints, 
                       const point_type* xc, VectorType* uv) {
    if (!uv) return;
    for (register int n=0; n<Npoints; n++) 
      uv[n] = _ErrorValue;

    int Np,Npr;
    int num = comm_service::proc_num();
    point_type *xcp = 0;
    PointsAllGather(Npoints,xc,Np,xcp);
    VectorType* uvp = new VectorType[Np*num]; 

    PointsValues(u,t,Np*num,xcp,uvp);
    delete [] xcp;

    VectorType* uvpr = 0; 
    DataAllScatter(Np,uvp,Npr,uvpr);
    for (int p=0; p<num; p++) 
      for (register int n=0; n<Npoints; n++) 
        if (uvpr[p*Npr+n] != _ErrorValue) 
          uv[n] = uvpr[p*Npr+n];

    delete [] uvpr;
    delete [] uvp;
  }

  void PointsValues(grid_fct_type& u, const double& t, const int& Npoints, 
                    const point_type* xc, DataType* uv, const point_type* xct=0) {

    if (Npoints<=0 || !uv) return;
    for (register int n=0; n<Npoints; n++) 
      uv[n] = _ErrorValue;

    for (register int l=0; l<=FineLevel(base::GH()); l++) {
      int Time = CurrentTime(base::GH(),l);
      PointsValues(u,Time,l,t,Npoints,xc,uv,xct);
    }
  }

  void PointsValues(grid_fct_type& u, const int& Time, const int& Level, const double& t,
                    const int& Npoints, const point_type* xc, DataType* uv, const point_type* xct=0) {

    if (Npoints<=0 || !uv) return;
    bool* lused = new bool[Npoints];

    forall (u,Time,Level,c)
      BBox bb = u.interiorbbox(Time,Level,c);
      int ncl;
      if (xct) 
        ncl = geom_base::LocalList(base::GH(),bb,Npoints,xct,lused);
      else 
        ncl = geom_base::LocalList(base::GH(),bb,Npoints,xc,lused);
      if (ncl>0) {
        register int n, nl;
        point_type* xcl = new point_type[ncl];
        DataType* uvl = new DataType[ncl];
        for (n=0, nl=0; n<Npoints; n++) 
          if (lused[n]) xcl[nl++] = xc[n];

        GetGrid(u(Time,Level,c),Time,Level,c,t,ncl,xcl,uvl);
          
        for (n=0, nl=0; n<Npoints; n++) 
          if (lused[n]) uv[n] = uvl[nl++];

        delete [] xcl;
        delete [] uvl;
      }
    end_forall

    delete [] lused;
  }    

  int LocalPoints(grid_fct_type& u, const int& Time, const int& Level, 
                  const int& Npoints, const point_type* xc, bool* lused) {
    int ncl = 0;
    if (Npoints<=0 || !xc) return ncl;

    for (register int n=0; n<Npoints; n++) 
      lused[n] = false;

    bool* lused_help = new bool[Npoints];
    forall (u,Time,Level,c)
      BBox bb = u.interiorbbox(Time,Level,c);
      int ncl_help = geom_base::LocalList(base::GH(),bb,Npoints,xc,lused_help);
      if (ncl_help>0)
        for (register int n=0; n<Npoints; n++) 
          if (lused_help[n]) lused[n] = true;
    end_forall
    delete [] lused_help;

    for (register int n=0; n<Npoints; n++) 
      if (lused[n]) ncl++;
    return ncl;
  }    

  virtual void GetGrid(grid_data_type& gdu, const int& Time, const int& Level, 
                       const int& c, const double& t, const int& Npoints, 
                       const point_type* xcl, DataType* uvl) {
    if (_interpolation)
      Interpolation_().Interpolate(base::GH(),gdu,Npoints,xcl,uvl,_ErrorValue);
  }

  void PointsValuesPar(grid_fct_type& u, const double& t, const int& Npoints, 
                       const point_type* xc, DataType* uv) {
    if (!uv) return;
    for (register int n=0; n<Npoints; n++) 
      uv[n] = _ErrorValue;

    int Np,Npr;
    int num = comm_service::proc_num();
    point_type* xcp = 0;
    PointsAllGather(Npoints,xc,Np,xcp);
    DataType* uvp = new DataType[Np*num]; 

    PointsValues(u,t,Np*num,xcp,uvp);
    delete [] xcp;

    DataType* uvpr = 0; 
    DataAllScatter(Np,uvp,Npr,uvpr);
    for (int p=0; p<num; p++) 
      for (register int n=0; n<Npoints; n++) 
        if (uvpr[p*Npr+n] != _ErrorValue) 
          uv[n] = uvpr[p*Npr+n];

    delete [] uvpr;
    delete [] uvp;
  }

  void PointsAllGather(int Nsnd, const point_type* snd, int& Nrcv, point_type*& rcv) {
    if (rcv) delete [] rcv;

    if (!comm_service::dce() || comm_service::proc_world() == 1) {
      Nrcv = Nsnd;
      rcv = new point_type[Nrcv]; 
      for (register int i=0; i<Nsnd; i++) rcv[i] = snd[i];
      return;
    }

#ifdef DAGH_NO_MPI
#else
    int num = comm_service::proc_num();
    int R = MPI_Allreduce(&Nsnd, &Nrcv, 1, MPI_INT, MPI_MAX, comm_service::comm());
    if ( MPI_SUCCESS != R ) 
      comm_service::error_die("AMRInterpolation::ArrayAllGather","MPI_Allreduce",R);
    rcv = new point_type[Nrcv*num]; 
    point_type* tmpsnd = new point_type[Nrcv*num];
    for (register int i=0; i<Nsnd; i++) tmpsnd[i] = snd[i];
    for (register int i=Nsnd; i<Nrcv; i++) tmpsnd[i] = _ErrorValue;
    int sndsize = Nrcv*sizeof(point_type);
    int rcvsize = Nrcv*sizeof(point_type);
    R = MPI_Allgather((void *)tmpsnd, sndsize, MPI_BYTE, (void *)rcv, rcvsize, 
                      MPI_BYTE, comm_service::comm()); 
    if ( MPI_SUCCESS != R ) 
        comm_service::error_die( "AMRInterpolation::PointsAllGather", "MPI_Allgather", R );
    delete [] tmpsnd;
#endif
  }

  void DataAllScatter(int Nsnd, const VectorType* snd, int& Nrcv, VectorType*& rcv) {
    Nrcv = Nsnd;
    if (rcv) delete [] rcv;

    if (!comm_service::dce() || comm_service::proc_world() == 1) {
      rcv = new VectorType[Nrcv]; 
      for (register int i=0; i<Nsnd; i++) rcv[i] = snd[i];
      return;
    }

#ifdef DAGH_NO_MPI
#else
    int num = comm_service::proc_num();
    rcv = new VectorType[Nrcv*num]; 
    int sndsize = Nsnd*sizeof(VectorType);
    int rcvsize = Nrcv*sizeof(VectorType);
    int R = MPI_Alltoall((void*) snd, sndsize, MPI_BYTE, (void *)rcv, rcvsize, 
                         MPI_BYTE, comm_service::comm()); 
    if ( MPI_SUCCESS != R ) 
      comm_service::error_die("AMRInterpolation::DataAllScatter","MPI_Alltoall",R);
#endif
  }

  void DataAllScatter(int Nsnd, const DataType* snd, int& Nrcv, DataType*& rcv) {
    Nrcv = Nsnd;
    if (rcv) delete [] rcv;

    if (!comm_service::dce() || comm_service::proc_world() == 1) {
      rcv = new DataType[Nrcv]; 
      for (register int i=0; i<Nsnd; i++) rcv[i] = snd[i];
      return;
    }

#ifdef DAGH_NO_MPI
#else
    int num = comm_service::proc_num();
    rcv = new DataType[Nrcv*num]; 
    int sndsize = Nsnd*sizeof(DataType);
    int rcvsize = Nrcv*sizeof(DataType);
    int R = MPI_Alltoall((void*) snd, sndsize, MPI_BYTE, (void *)rcv, rcvsize, 
                         MPI_BYTE, comm_service::comm()); 
    if ( MPI_SUCCESS != R ) 
      comm_service::error_die("AMRInterpolation::DataAllScatter","MPI_Alltoall",R);
#endif
  }

  void ArrayCombine(const int& TargetNode, const int& Nvalues, DataType* dat) {
#ifdef DAGH_NO_MPI
#else
    if (comm_service::dce() && comm_service::proc_num() > 1) {  
      int num = comm_service::proc_num();
      MPI_Status status;
      int R;
      int size = sizeof(DataType)*Nvalues;
      int tag = 10000;
      if (MY_PROC!=TargetNode) {
        R = MPI_Send((void *)dat, size, MPI_BYTE, TargetNode, tag, comm_service::comm());
        if ( MPI_SUCCESS != R ) 
          comm_service::error_die( "AMRInterpolation::ArrayCombine", "MPI_Send", R );
      }
      else 
        for (int proc=0; proc<num; proc++) {
          if (proc==TargetNode) continue;
          DataType* datp = new DataType[Nvalues];
          R = MPI_Recv((void *)datp, size, MPI_BYTE, proc, tag, comm_service::comm(), &status);
          if ( MPI_SUCCESS != R ) 
            comm_service::error_die( "AMRInterpolation::ArrayCombine", "MPI_Recv", R );
          for (register int n=0; n<Nvalues; n++) 
            if (datp[n] != _ErrorValue) 
              dat[n] = datp[n];
          delete [] datp;
        }
    }
#endif
  }

  inline void SetInterpolation(interpolation_type* inter) { 
    if (_interpolation) delete _interpolation;
    _interpolation = inter; 
  }
  inline interpolation_type& Interpolation_() { return *_interpolation; }
  inline const interpolation_type& Interpolation_() const { return *_interpolation; }

  inline void SetErrorValue(const DataType val) { _ErrorValue = val; }
  inline DataType ErrorValue() const { return _ErrorValue; }

protected:
  DataType _ErrorValue;
  interpolation_type* _interpolation;
};

#endif

---