CPTLevelSet.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROC_CPTLEVELSET_H
#define AMROC_CPTLEVELSET_H

#include "cpt/cpt.h"

template <class DataType, int dim>
class CPTLevelSet : public GFMLevelSet<DataType,dim> {
  typedef GFMLevelSet<DataType,dim> base;

public:
  typedef typename base::grid_fct_type grid_fct_type;   
  typedef typename base::grid_data_type grid_data_type;   
  typedef typename cpt::State<dim,DataType> cpt_type;

  CPTLevelSet() : base(), inverse(0), unsign(0), _FillWidth(0) {
    far_away = 1.e36;
  }

  virtual ~CPTLevelSet() {}

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"CPT");
    RegisterAt(base::LocCtrl,"PlotPhi",base::_PlotPhi);
    RegisterAt(base::LocCtrl,"Inverse",inverse);
    RegisterAt(base::LocCtrl,"Unsigned",unsign);
    RegisterAt(base::LocCtrl,"FillWidth",_FillWidth);
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual void SetBrep(const int num_vertices, const DataType* vertices,
                       const int num_edges, const int* edges) {
    Coords ex(dim,0);
    BBox wholebb = grow(coarsen(base::GH().wholebbox(),DAGHShadowFactor),base::NGhosts());
    DCoords lbcorner = base::GH().worldCoords(wholebb.lower(), wholebb.stepsize());
    DCoords ubcorner = base::GH().worldCoords(wholebb.upper()+wholebb.stepsize(), 
                                              wholebb.stepsize());
    DCoords dist = ubcorner-lbcorner; dist *= dist;
    DataType distance = 0.;
    for (register int n=0; n<dim; n++) {
      distance += dist(n);
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }
    distance = std::sqrt(distance);
    state.setParameters( domain, distance );
    state.setBRepWithNoClipping( num_vertices, reinterpret_cast<const void*>(vertices), 
                                 num_edges, reinterpret_cast<const void*>(edges) );
  }

  virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) { 
    if (state.getNumberOfGrids()>0) state.clearGrids();

    if (phi.len(Level)<=0) return;

    Coords ss = phi(Time,Level,0).stepsize();
    DCoords dx = base::GH().worldStep(ss);
    BBox bb = base::GH().wholebbox();
    if (ss(0)>bb.stepsize(0)) bb.coarsen(ss/bb.stepsize());
    else bb.refine(bb.stepsize()/ss);
    bb.grow(base::NGhosts());

    DCoords dxh2(dx); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = base::GH().worldCoords(bb.lower(), bb.stepsize()) + dxh2;
    DCoords ubcorner = base::GH().worldCoords(bb.upper(), bb.stepsize()) + dxh2;
    DataType ldomain[2*dim];
    for (register int n=0; n<dim; n++) {
      ldomain[n] = lbcorner(n); ldomain[dim+n] = ubcorner(n);
    }

    state.setLattice(bb.extents(),ldomain);

    forall (phi,Time,Level,c)
      grid_data_type& gdphi = phi(Time,Level,c);
      Coords lb = gdphi.lower(); 
      Coords ub = gdphi.upper()+gdphi.stepsize(); 
      lb /= gdphi.stepsize(); ub /= gdphi.stepsize();
      lb += base::NGhosts();  ub += base::NGhosts();      
      
      state.insertGrid( reinterpret_cast<const int*>(lb()), 
                        reinterpret_cast<const int*>(ub()), 
                        reinterpret_cast<DataType*> (gdphi.data()),     
                        reinterpret_cast<DataType*> (0), 
                        reinterpret_cast<DataType*> (0), 
                        reinterpret_cast<int*> (0) );
    end_forall

    if (state.getNumberOfGrids()>0) {
      if (_FillWidth>0) {
        double distance = 0.;
        for (register int n=0; n<dim; n++)
          distance += dx[n]*dx[n];
        distance = std::sqrt(distance)*_FillWidth;
        state.setParameters( domain, distance );
      }

      if (unsign) {
        if (state.areGridsValidUnsigned()) {
          START_WATCH
            state.computeClosestPointTransformUnsigned();
          END_WATCH(LS_CPT_TRANSFORM)
          START_WATCH
            state.floodFillUnsigned( far_away );
          END_WATCH(LS_CPT_FLOODFILL)
        }
        else 
          std::cout << "Invalid grid for closest point transform on lev=" << Level 
                    << " time=" << t << std::endl;
      }
      else {
        if (state.areGridsValid()) {
          START_WATCH
            state.computeClosestPointTransform();
          END_WATCH(LS_CPT_TRANSFORM)
          START_WATCH
            state.floodFillDetermineSign( far_away );
          END_WATCH(LS_CPT_FLOODFILL)
          if (inverse) {
            forall (phi,Time,Level,c)
              grid_data_type& gdphi = phi(Time,Level,c);
              gdphi.multiply(-1.0);
            end_forall
          }
        }
        else 
          std::cout << "Invalid grid for closest point transform on lev=" << Level 
                    << " time=" << t << std::endl;
      }
    }

    START_WATCH
      Sync(phi,Time,Level);
    END_WATCH(LS_SYNC)
    START_WATCH
      ExternalBoundaryUpdate(phi,Time,Level);
    END_WATCH(BOUNDARIES_EXTERNAL)
  }

  virtual void SetGrid(grid_data_type& gdphi, const int& Level, const double& t) {
    if (state.getNumberOfGrids()>0) state.clearGrids();

    BBox bb = gdphi.bbox();
    Coords ss = bb.stepsize();
    DCoords dx = base::GH().worldStep(ss);
    DCoords dxh2(dx); dxh2 *= static_cast<DataType>(0.5);
    DCoords lbcorner = base::GH().worldCoords(bb.lower(), ss) + dxh2;
    DCoords ubcorner = base::GH().worldCoords(bb.upper(), ss) + dxh2;
    DataType ldomain[2*dim];
    for (register int n=0; n<dim; n++) {
      ldomain[n] = lbcorner(n); ldomain[dim+n] = ubcorner(n);
    }

    state.setLattice(bb.extents(),ldomain);

    Coords lb = gdphi.lower(); 
    Coords ub = gdphi.upper()+ss; 
    lb /= ss; ub /= ss;
    lb += base::NGhosts();  ub += base::NGhosts();      
      
    state.insertGrid( reinterpret_cast<const int*>(lb()), 
                      reinterpret_cast<const int*>(ub()), 
                      reinterpret_cast<DataType*> (gdphi.data()),       
                      reinterpret_cast<DataType*> (0), 
                      reinterpret_cast<DataType*> (0), 
                      reinterpret_cast<int*> (0) );

    if (state.getNumberOfGrids()>0) {
      if (_FillWidth>0) {
        double distance = 0.;
        for (register int n=0; n<dim; n++)
          distance += dx[n]*dx[n];
        distance = std::sqrt(distance)*_FillWidth;
        state.setParameters( domain, distance );
      }

      if (unsign) {
        if (state.areGridsValidUnsigned()) {
          START_WATCH
            state.computeClosestPointTransformUnsigned();
          END_WATCH(LS_CPT_TRANSFORM)
          START_WATCH
            state.floodFillUnsigned( far_away );
          END_WATCH(LS_CPT_FLOODFILL)
        }
        else 
          std::cout << "Invalid grid for closest point transform on lev=" << Level 
                    << " time=" << t << std::endl;
      }
      else {
        if (state.areGridsValid()) {
          START_WATCH
            state.computeClosestPointTransform();
          END_WATCH(LS_CPT_TRANSFORM)
          START_WATCH
            state.floodFillDetermineSign( far_away );
          END_WATCH(LS_CPT_FLOODFILL)
          if (inverse)
            gdphi.multiply(-1.0);
        }
        else 
          std::cout << "Invalid grid for closest point transform on lev=" << Level 
                    << " time=" << t << std::endl;
      }
    }
  }

  inline cpt_type& State() { return state; }
  inline const cpt_type& State() const { return state; }
  inline void SetUnsigned(const int us) { unsign = us; }
  inline int Unsigned() const { return unsign; }
  inline void SetInverse(const int iv) { inverse = iv; }
  inline int Inverse() const { return inverse; }

protected:
  cpt_type state;
  int inverse, unsign, _FillWidth; 
  DataType far_away;
  DataType domain[2*dim];
};


#endif

---