amroc/amr/CPTLevelSet.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_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), 
                  used_num_vertices(0), used_num_connections(0), used_vertices(0), used_connections(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_connections, const int* connections) {
    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_connections, reinterpret_cast<const void*>(connections) );
    
    used_num_vertices = num_vertices;
    used_vertices = vertices;
    used_num_connections = num_connections;
    used_connections = connections;
  }

  bool GetBrep(int& num_vertices, const DataType*& vertices, int& num_connections, const int*& connections) {
    num_vertices = used_num_vertices;
    vertices = used_vertices;
    num_connections = used_num_connections;
    connections = used_connections;
    return ((num_vertices>0) && (num_connections>0) && (vertices!=0) && (connections!=0));
  }

  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)

    // Correct potential mistakes from heuristic algorithm inside floodFillDetermineSign()
    if (!unsign) {
      if (base::Dim()==1) {
        forall (phi,Time,Level,c)
          grid_data_type& gdphi = phi(Time,Level,c);
          BBox bb = gdphi.bbox();
          Coords ss = bb.stepsize();
          BeginFastIndex3(ph, bb, gdphi.data(), DataType);
          bb.grow(-base::NGhosts());
          for_1 (n, bb, ss)
            if (FastIndex1(ph,n)==-far_away) {
              if (FastIndex1(ph,n-ss(0))>=0 || FastIndex1(ph,n+ss(0))>=0)
                FastIndex1(ph,n)=far_away;
            }
          end_for      
          EndFastIndex1(ph);
        end_forall
      }
      else if (base::Dim()==2) {
        forall (phi,Time,Level,c)
          grid_data_type& gdphi = phi(Time,Level,c);
          BBox bb = gdphi.bbox();
          Coords ss = bb.stepsize();
          BeginFastIndex2(ph, bb, gdphi.data(), DataType);
          bb.grow(-base::NGhosts());
          for_2 (n, m, bb, ss)
            if (FastIndex2(ph,n,m)==-far_away) {
              if (FastIndex2(ph,n-ss(0),m)>=0 || FastIndex2(ph,n,m-ss(1))>=0 || 
                  FastIndex2(ph,n+ss(0),m)>=0 || FastIndex2(ph,n,m+ss(1))>=0)
                FastIndex2(ph,n,m)=far_away;
            }
          end_for      
          EndFastIndex2(ph);
        end_forall
      }
      else if (base::Dim()==3) {
        forall (phi,Time,Level,c)
          grid_data_type& gdphi = phi(Time,Level,c);
          BBox bb = gdphi.bbox();
          Coords ss = bb.stepsize();
          BeginFastIndex3(ph, bb, gdphi.data(), DataType);
          bb.grow(-base::NGhosts());
          for_3 (n, m, l, bb, ss)
            if (FastIndex3(ph,n,m,l)==-far_away) {
              if (FastIndex3(ph,n-ss(0),m,l)>=0 || FastIndex3(ph,n,m-ss(1),l)>=0 || FastIndex3(ph,n,m,l-ss(2))>=0 ||
                  FastIndex3(ph,n+ss(0),m,l)>=0 || FastIndex3(ph,n,m+ss(1),l)>=0 || FastIndex3(ph,n,m,l+ss(2))>=0)
                FastIndex3(ph,n,m,l)=far_away;
            }
          end_for      
          EndFastIndex3(ph);
        end_forall
      }
    }
    
    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; 
  int used_num_vertices, used_num_connections;
  const DataType* used_vertices;
  const int* used_connections;
  DataType far_away;
  DataType domain[2*dim];
};


#endif

---