amroc/amr/GFMBoundary.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@cacr.caltech.edu

#ifndef AMROC_GFMBOUNDARYBASE_H
#define AMROC_GFMBOUNDARYBASE_H

#define REFLECTION     2
#define EXTRAPOLATION  1

#include "GFMGeom.h"
#include "Interpolation.h"

template <class VectorType, int dim>
class GFMBoundaryBase : public AMRBase<VectorType,dim>, 
                        public Geom<typename VectorType::InternalDataType,dim>,
                        public GFMGeom<typename VectorType::InternalDataType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  
  typedef AMRBase<VectorType,dim> base;
  typedef Geom<DataType,dim> geom_base;
  typedef GFMGeom<DataType,dim> gfm_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 GridData<DataType,dim> grid_data_type;   
  typedef GridFunction<DataType,dim> grid_fct_type;  
  typedef GridData<bool,dim> bool_grid_data_type;
  typedef GridFunction<bool,dim> bool_grid_fct_type;  
  typedef Vector<DataType,dim> point_type;

  typedef Interpolation<VectorType,dim> interpolation_type;

  GFMBoundaryBase() : base(), geom_base(), _BoundaryWidth(1), _bndtype(REFLECTION), _naux(0), _Verbose(1),
                      _Cutoff(0.0), _ErrorValue(-1.e37), _normal_factor(2.0) {
    _interpolation = new interpolation_type(); 
    _far_away = 1.e36;    
  }

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

  //******************************************************************************
  // Abstract class interface
  //******************************************************************************
  virtual void SetGrid(vec_grid_data_type& gdu, grid_data_type& gdphi, 
                       const BBox& bb, const int& Level, double t, 
                       const int& nc, const int* idx, const point_type* xc, 
                       DataType* distance, point_type* normal) = 0;  
  virtual void GetBndryCells(vec_grid_fct_type& u, grid_data_type& gdphi, bool_grid_data_type& gdbf, 
                             const BBox& bb, const int& Time, const int& Level, const int& c, 
                             int& nc, int*& idx, int& ni, int*& idx_wrg) = 0;
  virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi, const BBox& bb, 
                           const VectorType* u, DataType* aux, const int& Level, 
                           double t, const int& nc, const int* idx, 
                           const point_type* xc, const DataType* distance, 
                           const point_type* normal) {}
  //******************************************************************************

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"InternalBoundary");
    RegisterAt(base::LocCtrl,"BoundaryWidth",_BoundaryWidth);
    RegisterAt(base::LocCtrl,"Cutoff",_Cutoff);
    RegisterAt(base::LocCtrl,"ErrorValue",_ErrorValue);
    RegisterAt(base::LocCtrl,"BndType",_bndtype);
    RegisterAt(base::LocCtrl,"NAux",_naux);
    RegisterAt(base::LocCtrl,"Verbose",_Verbose);    
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }
  virtual void init() {}
  virtual void update() {}
  virtual void finish() {} 
  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
    base::SetupData(gh, ghosts);
    if (_bndtype == EXTRAPOLATION) _normal_factor = 1.0;
    else if (_bndtype == REFLECTION) _normal_factor = 2.0;
    else _normal_factor = 0.0;
    _Cutoff -= 1.e-12;
  }

  void SetPhysbd(vec_grid_fct_type& u, grid_fct_type& phi, 
                 bool_grid_fct_type& bf, const int Time, 
                 const int Level, double t, bool extrapolate=false) { 
    START_WATCH
    forall (u,Time,Level,c)
      int nc, ni;
      int* idx = (int *) 0;
      int* idx_wrg = (int *) 0;
      BBox bb = u.interiorbbox(Time,Level,c);
      bb.grow(std::min(NBoundaryWidth(),base::NGhosts()-1));

      START_WATCH
        GetBndryCells(u,phi(Time,Level,c),bf(Time,Level,c),bb,Time,Level,c,
                      nc,idx,ni,idx_wrg);
      END_WATCH(GFM_FINDING_CELLS)

      if (_Verbose && (idx_wrg!=(int *) 0)) {
        for (register int cnt=0; cnt<ni; cnt++) {
          std::cout << "Unrecoverable ambigious ghost cell " 
                    << Coords(base::Dim(),&(idx_wrg[base::Dim()*cnt])) 
                    << " on Level " << Level << std::endl;
#ifdef DEBUG_PRINT      
          ( comm_service::log() << "GFMBoundaryBase: " << bb 
            << " Unrecoverable ambigious ghost cell " 
            << Coords(base::Dim(),&(idx_wrg[base::Dim()*cnt])) 
            << " on Level " << Level << "\n" ).flush();
#endif
        }
      }

      point_type* xc = new point_type[nc];
      point_type* normal = new point_type[nc];
      DataType* distance = new DataType[nc];

      START_WATCH
        geom_base::CellCoords(base::GH(),bb.stepsize(),nc,idx,xc);
        gfm_geom_base::CellNormals(base::GH(),phi(Time,Level,c),nc,idx,normal);
      END_WATCH(GFM_GEOMETRY)

      phi(Time,Level,c).plus(Cutoff());

      START_WATCH
        gfm_geom_base::CellDistances(base::GH(),phi(Time,Level,c),nc,idx,distance);
      END_WATCH(GFM_GEOMETRY)

      if (extrapolate) {
        DataType factor = 1.0;
        VectorType* uv = (VectorType*) 0;
        Extrapolation(u(Time,Level,c),phi(Time,Level,c),uv,Level,nc,idx,xc,
                      distance,normal,factor);
        for (register int cnt=0; cnt<nc; cnt++) 
          u(Time,Level,c)(Coords(dim,idx+dim*cnt)) = uv[cnt];
        if (uv) delete [] uv;
      }
      else
        SetGrid(u(Time,Level,c),phi(Time,Level,c),bb,Level,t,
                nc,idx,xc,distance,normal);  

      phi(Time,Level,c).minus(Cutoff());

      if (idx) delete [] idx;
      if (idx_wrg) delete [] idx_wrg;
      if (xc) delete [] xc;
      if (normal) delete [] normal;
      if (distance) delete [] distance;

    end_forall
    END_WATCH(GFM_SETBNDRY_WHOLE)

    START_WATCH
      Sync(u,Time,Level);
    END_WATCH(BOUNDARIES_SYNC)
    START_WATCH
      ExternalBoundaryUpdate(u,Time,Level);
    END_WATCH(BOUNDARIES_EXTERNAL)
  }

  virtual void Extrapolation(vec_grid_data_type& gdu, grid_data_type& gdphi, 
                             VectorType*& u, const int& Level, 
                             const int& nc, const int* idx, 
                             const point_type* xc, DataType* distance, 
                             point_type* normal, const DataType& normal_factor) {
    
    START_WATCH
    if (! ValidNormalFactor(normal_factor)) {
      std::cout << base::NGhosts() << " ghost cells too few for interpolation. "
                << "At least " << int(NBoundaryWidth()*(normal_factor+1)) 
                << " cells recommended.\n";
#ifdef DEBUG_PRINT
      ( comm_service::log() << "GFMBoundaryBase: " << base::NGhosts() 
        << " ghost cells too few for interpolation. " << "At least " 
        << int(NBoundaryWidth()*(normal_factor+1)) << " cells recommended.\n" ).flush();
#endif
    }

    if (u) delete [] u;
    u = new VectorType[nc];
    point_type* pos = new point_type[nc];

    register int n;
    for (n=0; n<nc; n++) {
      pos[n] = xc[n]+normal_factor*distance[n]*normal[n];
      u[n] = _ErrorValue;
    }

    if (_interpolation)
      Interpolation_().Interpolate(base::GH(),gdu,gdphi,nc,pos,u,_ErrorValue);

    for (n = 0; n<nc; n++) 
      if (u[n](0) == _ErrorValue) {
        Coords co(dim,idx+dim*n);
        u[n] = gdu(co);
#ifdef DEBUG_PRINT
        if (_Verbose) {
          ( comm_service::log() << "GFMBoundaryBase: " << gdu.bbox() 
            << " Interpolation error in " << co 
            << "=" << xc[n] << " from " << pos[n] << " on Level " 
            << Level << "\n" ).flush();
          ( comm_service::log() << "                 " << distance[n] 
            << " " << normal[n] << "\n" ).flush();
        }
#endif
        normal[n] = static_cast<DataType>(0.0);
        distance[n] = static_cast<DataType>(0.0);
      }
    delete [] pos;
    END_WATCH(GFM_EXTRAPOLATION)
  }

  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 const int& NAux() const { return _naux; }
  inline void SetNAux(const int naux) { _naux = naux; }

  inline bool ValidNormalFactor(const DataType& nf) const 
  { return (NBoundaryWidth()*(nf+1) <= base::NGhosts()); }
  inline const DataType& NormalFactor() const { return _normal_factor; }
  inline const int& BndType() const { return _bndtype; }
  inline int isgn(const DataType val)
  { return (val>0 ? 1 : (val<0 ? -1 : 0)); }
  inline const int& NBoundaryWidth() const { return _BoundaryWidth; }
  inline const DataType& Cutoff() const { return _Cutoff; }
  inline const DataType& FarAway() const { return _far_away; }
  inline void SetVerbose(const int v) { _Verbose = v; }  
  inline const int& GetVerbose() const { return _Verbose; }

protected:
  int _BoundaryWidth, _bndtype, _naux, _Verbose;
  DataType _Cutoff, _ErrorValue, _normal_factor, _far_away;
  interpolation_type* _interpolation;
};

template <class VectorType, int dim> class GFMBoundary;


template <class VectorType>
class GFMBoundary<VectorType,2> : public GFMBoundaryBase<VectorType,2> {
  typedef typename VectorType::InternalDataType DataType;
  typedef GFMBoundaryBase<VectorType,2> base; 

public:
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::grid_data_type grid_data_type;
  typedef typename base::bool_grid_data_type bool_grid_data_type;
  typedef typename base::point_type point_type;

  GFMBoundary() : base() {}

  virtual ~GFMBoundary() {}

  virtual void GetBndryCells(vec_grid_fct_type& u, grid_data_type& gdphi, bool_grid_data_type& gdbf, 
                             const BBox& bb, const int& Time, const int& Level, const int& c,
                             int& nc, int*& idx, int& Nillegal, int*& idx_wrg) {

    GridData<char,2> flag(gdphi.bbox());
    flag.equals(0);
    bool_grid_data_type pf(gdphi.bbox());
    BeginFastIndex2(phi, gdphi.bbox(), gdphi.data(), DataType);
    BeginFastIndex2(flag, flag.bbox(), flag.data(), char);
    BeginFastIndex2(bflag, gdbf.bbox(), gdbf.data(), bool);
    BeginFastIndex2(pf, pf.bbox(), pf.data(), bool);
    
    Sgn(gdphi,pf,base::Cutoff());

    Coords gdss(gdphi.bbox().stepsize());
    DCoords dx = base::GH().worldStep(gdss);    
    double distance = 0.;
    for (register int n=0; n<base::Dim(); n++)
      distance += dx[n]*dx[n];
    distance = std::sqrt(distance)-base::Cutoff();
    for_2 (n, m, gdphi.bbox(), gdss)
      if (FastIndex2(phi,n,m)<-base::Cutoff() || 
          FastIndex2(phi,n,m)>distance) continue;
      BBox mark(2,n,m,n,m,gdss(0),gdss(1));
      mark.grow(base::NBoundaryWidth());
      flag.equals(1,mark);
    end_for

    Coords ss(bb.stepsize());
    nc = 0;   
    Nillegal = 0;   
    for_2 (n, m, bb, ss)
      if (FastIndex2(pf,n,m))
        FastIndex2(flag,n,m) = 0;
      else {
        if (FastIndex2(flag,n,m)==1) {
          if (FastIndex2(pf,n+ss(0),m) && FastIndex2(pf,n-ss(0),m) &&
              FastIndex2(pf,n,m+ss(1)) && FastIndex2(pf,n,m-ss(1))) {
            if (u.interiorbbox(Time,Level,c).inside(n,m)) {
              bool illegal_point = true;
              for (register int ci=0; ci<u.children(Time,Level,c); ci++) 
                if (u.interiorbbox(CurrentTime(base::GH(),Level+1),Level+1,
                                   u.childidx(Time,Level,c,ci)).inside(n,m)) {
                  illegal_point = false;
                  break;
                }      
              if (illegal_point) FastIndex2(flag,n,m) = -1;
            }
            else 
              FastIndex2(flag,n,m) = -2;
          }
          if (FastIndex2(flag,n,m)==-1)
            Nillegal++;
          else if (FastIndex2(flag,n,m)==1)
            nc++;
        }
        else 
          if (FastIndex2(phi,n,m)>=-base::FarAway() && 
              FastIndex2(phi,n,m)<-base::Cutoff())
          FastIndex2(flag,n,m) = 2;
      }
    end_for
      
    if (idx) delete [] idx;
    if (idx_wrg) delete [] idx_wrg;
    
    idx = new int[2*nc];
    if (base::_Verbose) 
      idx_wrg = new int[2*Nillegal];
    int cnti = 0, cntill = 0;
   
    for_2 (n, m, gdbf.bbox(), ss)
      if (bb.inside(n,m)) {
        if (FastIndex2(flag,n,m)>0)
          FastIndex2(bflag,n,m) = true; 
        if (FastIndex2(flag,n,m)==1) {
          idx[2*cnti] = n; idx[2*cnti+1] = m; cnti++;
        }
        else if (FastIndex2(flag,n,m)==-1 && base::_Verbose) { 
          idx_wrg[2*cntill] = n; idx_wrg[2*cntill+1] = m; cntill++;
        }
      }
      else
        FastIndex2(bflag,n,m) = true; 
    end_for
 
    EndFastIndex2(pf);
    EndFastIndex2(bflag);
    EndFastIndex2(flag);
    EndFastIndex2(phi);
  }

  virtual void Sgn(grid_data_type& gdphi, bool_grid_data_type& gdbf, const DataType& c) {
    BBox OpBox = gdphi.bbox();
    Coords& OpBox_stepsize = OpBox.stepsize();
    BeginFastIndex2(phi, OpBox, gdphi.data(), DataType);
    BeginFastIndex2(bf, OpBox, gdbf.data(), bool);
    for_2 (n, m, OpBox, OpBox_stepsize)
      FastIndex2(bf,n,m) = (FastIndex2(phi,n,m)>=-c); 
    end_for
    EndFastIndex2(bf);
    EndFastIndex2(phi);
  }
};

template <class VectorType>
class GFMBoundary<VectorType,3> : public GFMBoundaryBase<VectorType,3> {
  typedef typename VectorType::InternalDataType DataType;
  typedef GFMBoundaryBase<VectorType,3> base; 

public:
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::grid_data_type grid_data_type;
  typedef typename base::bool_grid_data_type bool_grid_data_type;
  typedef typename base::point_type point_type;

  GFMBoundary() : base() {}

  virtual ~GFMBoundary() {}

  virtual void GetBndryCells(vec_grid_fct_type& u, grid_data_type& gdphi, bool_grid_data_type& gdbf, 
                             const BBox& bb, const int& Time, const int& Level, const int& c,
                             int& nc, int*& idx, int& Nillegal, int*& idx_wrg) {

    GridData<char,3> flag(gdphi.bbox());
    flag.equals(0);
    bool_grid_data_type pf(gdphi.bbox());
    BeginFastIndex3(phi, gdphi.bbox(), gdphi.data(), DataType);
    BeginFastIndex3(flag, flag.bbox(), flag.data(), char);
    BeginFastIndex3(bflag, gdbf.bbox(), gdbf.data(), bool);
    BeginFastIndex3(pf, pf.bbox(), pf.data(), bool);

    Sgn(gdphi,pf,base::Cutoff());

    Coords gdss(gdphi.bbox().stepsize());
    DCoords dx = base::GH().worldStep(gdss);    
    double distance = 0.;
    for (register int n=0; n<base::Dim(); n++)
      distance += dx[n]*dx[n];
    distance = std::sqrt(distance)-base::Cutoff();
    for_3 (n, m, l, gdphi.bbox(), gdss)
      if (FastIndex3(phi,n,m,l)<-base::Cutoff() ||
          FastIndex3(phi,n,m,l)>distance) continue;
      BBox mark(3,n,m,l,n,m,l,gdss(0),gdss(1),gdss(2));
      mark.grow(base::NBoundaryWidth());
      flag.equals(1,mark);
    end_for

    Coords ss(bb.stepsize());
    nc = 0;   
    Nillegal = 0;   
    for_3 (n, m, l, bb, ss)
      if (FastIndex3(pf,n,m,l))
        FastIndex3(flag,n,m,l) = 0;
      else {
        if (FastIndex3(flag,n,m,l)==1) {
          if (FastIndex3(pf,n+ss(0),m,l) && FastIndex3(pf,n-ss(0),m,l) &&
              FastIndex3(pf,n,m+ss(1),l) && FastIndex3(pf,n,m-ss(1),l) &&
              FastIndex3(pf,n,m,l+ss(2)) && FastIndex3(pf,n,m,l-ss(2))) {
            if (u.interiorbbox(Time,Level,c).inside(n,m,l)) {
              bool illegal_point = true;
              for (register int ci=0; ci<u.children(Time,Level,c); ci++) 
                if (u.interiorbbox(CurrentTime(base::GH(),Level+1),Level+1,
                                   u.childidx(Time,Level,c,ci)).inside(n,m,l)) {
                  illegal_point = false;
                  break;
                }      
              if (illegal_point) FastIndex3(flag,n,m,l) = -1;
            }
            else 
              FastIndex3(flag,n,m,l) = -2;
          }
          if (FastIndex3(flag,n,m,l)==-1)
            Nillegal++;
          else if (FastIndex3(flag,n,m,l)==1)
            nc++;
        }
        else if (FastIndex3(phi,n,m,l)>=-base::FarAway() && 
                 FastIndex3(phi,n,m,l)<-base::Cutoff())
          FastIndex3(flag,n,m,l) = 2;
      }
    end_for
      
    if (idx) delete [] idx;
    if (idx_wrg) delete [] idx_wrg;
    
    idx = new int[3*nc];
    if (base::_Verbose)
      idx_wrg = new int[3*Nillegal];
    int cnti = 0, cntill = 0;
   
    for_3 (n, m, l, gdbf.bbox(), ss)
      if (bb.inside(n,m,l)) {
        if (FastIndex3(flag,n,m,l)>0)
          FastIndex3(bflag,n,m,l) = true; 
        if (FastIndex3(flag,n,m,l)==1) {
          idx[3*cnti] = n; idx[3*cnti+1] = m; idx[3*cnti+2] = l; cnti++;
        }
        else if (FastIndex3(flag,n,m,l)==-1 && base::_Verbose) { 
          idx_wrg[3*cntill] = n; idx_wrg[3*cntill+1] = m; 
          idx_wrg[3*cntill+2] = l; cntill++;
        }
      }
      else
        FastIndex3(bflag,n,m,l) = true; 
    end_for

    EndFastIndex3(pf);
    EndFastIndex3(bflag);
    EndFastIndex3(flag);
    EndFastIndex3(phi);
  }

  virtual void Sgn(grid_data_type& gdphi, bool_grid_data_type& gdbf, const DataType& c) {
    BBox OpBox = gdphi.bbox();
    Coords& OpBox_stepsize = OpBox.stepsize();
    BeginFastIndex3(phi, OpBox, gdphi.data(), DataType);
    BeginFastIndex3(bf, OpBox, gdbf.data(), bool);
    for_3 (n, m, l, OpBox, OpBox_stepsize)
      FastIndex3(bf,n,m,l) = (FastIndex3(phi,n,m,l)>=-c); 
    end_for
    EndFastIndex3(bf);
    EndFastIndex3(phi);
  }
};

#endif

---