StatCPTLevelSet.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_STATCPTLEVELSET_H
#define AMROC_STATCPTLEVELSET_H

#include "cpt/cpt.h"

template <class DataType, int dim>
class StatCPTLevelSet : 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;
  typedef ads::FixedArray<dim,DataType> point_type;
  typedef ads::FixedArray<dim,int> multi_index_type;

  StatCPTLevelSet() : base(), brep_filetype(0), inverse(0), 
                      unsign(0), _FillWidth(0) {
    max_distance = std::numeric_limits<DataType>::max();
    base::_Stationary = 1;
    far_away = 1.e36;
  }

  virtual ~StatCPTLevelSet() {}

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"StatCPT");
    RegisterAt(base::LocCtrl,"Brep",brep_filename);
    RegisterAt(base::LocCtrl,"FileType",brep_filetype);
    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 SetupData(GridHierarchy* gh, const int& ghosts) {
    base::SetupData(gh, ghosts);

    std::ifstream brep_file( brep_filename.c_str() );
    if ( ! brep_file ) {
      std::cerr << "Could not read the b-rep file.  Exiting...\n";
      std::exit( 1 );
    }

    int num_vertices, num_edges;
    point_type* vertices;
    multi_index_type* edges;

    // ISS file
    if (brep_filetype==1) {
      int fdim, mdim;
      brep_file >> fdim >> mdim;
      if ( fdim!=dim || mdim!=dim-1 ) {
        std::cerr << "Mesh in b-rep file has wrong dimension.  Exiting...\n";
        std::exit( 1 );
      }
      brep_file >> num_vertices;
      vertices = new point_type[ num_vertices ];
      const point_type* point_iter_end = vertices + num_vertices;
      for ( point_type* piter = vertices; piter != point_iter_end; ++piter )
        brep_file >> *piter;
      brep_file >> num_edges;
      edges = new multi_index_type[ num_edges ];
    }
    else {
      brep_file >> num_vertices >> num_edges;
      vertices = new point_type[ num_vertices ];
      edges = new multi_index_type[ num_edges ];
      const point_type* point_iter_end = vertices + num_vertices;
      for ( point_type* piter = vertices; piter != point_iter_end; ++piter )
        brep_file >> *piter;
    }
      
    const multi_index_type* edges_iter_end = edges + num_edges;
    for ( multi_index_type* eiter = edges; eiter != edges_iter_end; ++eiter ) 
      brep_file >> *eiter;

    assert( brep_file.good() );
    brep_file.close();
    
    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;
    max_distance = 0.;
    for (register int n=0; n<dim; n++) {
      max_distance += dist(n);
      domain[n] = lbcorner(n); domain[dim+n] = ubcorner(n);
    }
    max_distance = std::sqrt(max_distance);

    state.setParameters( domain, max_distance );
    state.setBRepWithNoClipping( num_vertices, reinterpret_cast<const void*>(vertices), 
                                 num_edges, reinterpret_cast<const void*>(edges) );

    delete[] vertices;
    delete[] 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 SetFilename(const std::string filename) { brep_filename = filename; }
  inline std::string Filename() { return brep_filename; }
  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;
  std::string brep_filename;
  DataType max_distance, far_away;
  int brep_filetype, inverse, unsign, _FillWidth; 
  DataType domain[2*dim];
};


#endif

---