amroc/amr/AMRGFMInterpolation.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_AMRGFMINTERPOLATION_H
#define AMROC_AMRGFMINTERPOLATION_H

#include "AMRGFMSolver.h"
#include "AMRInterpolation.h"
#include "GFMGeom.h"

template <class VectorType, class FixupType, class FlagType, int dim>
class AMRGFMInterpolation : public AMRInterpolation<VectorType,dim>,
                            public GFMGeom<typename VectorType::InternalDataType,dim> {
  typedef typename VectorType::InternalDataType DataType;
  typedef AMRInterpolation<VectorType,dim> base;
  typedef GFMGeom<DataType,dim> gfm_geom_base;

public:
  typedef AMRGFMSolver<VectorType,FixupType,FlagType,dim> gfm_solver_type;
  typedef GhostFluidMethod<VectorType,dim> gfm_type;
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::vec_grid_data_type vec_grid_data_type;  
  typedef typename base::grid_fct_type grid_fct_type;   
  typedef typename base::grid_data_type grid_data_type;   
  typedef typename base::point_type point_type;
  typedef typename base::interpolation_type interpolation_type;
  typedef GridFunction<bool,dim> bool_grid_fct_type;  

  AMRGFMInterpolation(gfm_solver_type& solver) : base(), _solver(solver) {}

  virtual ~AMRGFMInterpolation() {}

  void PointsGeom(grid_fct_type& phi, const double& t, const int& Npoints, 
                  const point_type* xc, point_type* normal, DataType* distance) {

    if (Npoints<=0 || (!normal && !distance)) return;
    for (register int n=0; n<Npoints; n++) {
      if (normal) normal[n] = base::ErrorValue();
      if (distance) distance[n] = base::ErrorValue();
    }

    for (register int l=0; l<=FineLevel(base::GH()); l++) {
      int Time = CurrentTime(base::GH(),l);
      PointsGeom(phi,Time,l,t,Npoints,xc,normal,distance);
    }
  }

  void PointsGeom(grid_fct_type& phi, const int& Time, const int& Level, 
                  const double& t, const int& Npoints, const point_type* xc, 
                  point_type* normal, DataType* distance) {

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

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

        base::CellIndices(base::GH(),bb.stepsize(),ncl,xcl,idxl);
        if (normal) {
          point_type* normall = new point_type[ncl];
          gfm_geom_base::CellNormals(base::GH(),phi(Time,Level,c),ncl,idxl,normall);
          for (n=0, nl=0; n<Npoints; n++) 
            if (lused[n]) normal[n] = normall[nl++];
          delete [] normall;
        }
        if (distance) {
          DataType* distancel = new DataType[ncl];
          gfm_geom_base::CellDistances(base::GH(),phi(Time,Level,c),ncl,idxl,distancel);
          for (n=0, nl=0; n<Npoints; n++) 
            if (lused[n]) distance[n] = distancel[nl++];
          delete [] distancel;
        }
        delete [] idxl;
        delete [] xcl;
      }
    end_forall

    delete [] lused;
  }    

  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 (!base::_interpolation) return;
    if (NGFM()>0) 
      base::Interpolation_().Interpolate(base::GH(),gdu,BF()(Time,Level,c),
                                         Npoints,xcl,uvl,base::ErrorValue());
    else
      base::Interpolation_().Interpolate(base::GH(),gdu,Npoints,xcl,uvl,base::ErrorValue());      
  }

  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 (!base::_interpolation) return;
    if (NGFM()>0)
      base::Interpolation_().Interpolate(base::GH(),gdu,BF()(Time,Level,c),
                                         Npoints,xcl,uvl,base::ErrorValue());
    else
      base::Interpolation_().Interpolate(base::GH(),gdu,Npoints,xcl,uvl,base::ErrorValue());      
  }

  void PointsGeomPar(grid_fct_type& phi, const double& t, const int& Npoints, 
                     const point_type* xc, point_type* normal, DataType* distance) {
    if (!normal && !distance) return;
    for (register int n=0; n<Npoints; n++) {
      if (normal) normal[n] = base::ErrorValue();
      if (distance) distance[n] = base::ErrorValue();
    }

    int Np,Npr;
    int num = comm_service::proc_num();
    point_type *xcp = 0;
    base::PointsAllGather(Npoints,xc,Np,xcp);
    DataType* geomp = new DataType[(dim+1)*Np*num];

    PointsGeom(phi,t,Np*num,xcp,(point_type *)geomp,geomp+dim*Np*num);
    delete [] xcp;

    DataType* geompr = 0; 
    base::DataAllScatter((dim+1)*Np,geomp,Npr,geompr);
    for (int p=0; p<num; p++) 
      for (register int n=0; n<Npoints; n++) {
        if (normal) {
          point_type* normalpr = (point_type *)(geompr+p*Npr+dim*n);
          if ((*normalpr) != base::ErrorValue()) 
            normal[n] = (*normalpr);
        }
        
        if (distance) {
          DataType* distancepr = geompr+p*Npr+dim*(Npr/(dim+1))+n;
          if ((*distancepr) != base::ErrorValue()) 
            distance[n] = (*distancepr);
        }
      }
    
    delete [] geompr;
    delete [] geomp;
  }

  inline gfm_type& GFM(const int& n) { return _solver.GFM(n); }
  inline const int& NGFM() { return _solver.NGFM(); }
  inline bool_grid_fct_type& BF() { return _solver.BF(); }
  inline const bool_grid_fct_type& BF() const{ return _solver.BF(); }

protected:
  gfm_solver_type& _solver;
};

#endif

---