GFMGeom.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROC_GFM_GEOM_H
#define AMROC_GFM_GEOM_H

#include "Geom.h"

template <class DataType, int dim>
class GFMGeomBase {
public:
  typedef GridData<DataType,dim> grid_data_type;   
  typedef Vector<DataType,dim> point_type;

  virtual ~GFMGeomBase() {}

  //******************************************************************************
  // Abstract class interface
  //******************************************************************************
  virtual void CellNormals(GridHierarchy& GH, grid_data_type& gdphi, const int& nc, 
                           const int* idx, point_type* normals) = 0;
  virtual void CellDistances(GridHierarchy& GH, grid_data_type& gdphi, const int& nc, 
                             const int* idx, DataType* distance) = 0;
  //******************************************************************************
};

template <class DataType, int dim> class GFMGeom;

template <class DataType>
class GFMGeom<DataType,2> : public GFMGeomBase<DataType,2> {
  typedef GFMGeomBase<DataType,2> base; 

public:
  typedef typename base::grid_data_type grid_data_type;
  typedef typename base::point_type point_type;

  GFMGeom() : base() {}

  virtual ~GFMGeom() {}

  virtual void CellNormals(GridHierarchy& GH, grid_data_type& gdphi, const int& nc, 
                           const int* idx, point_type* normals) {
    Coords ss(gdphi.stepsize());
    DCoords dx2(GH.worldStep(ss));  
    dx2 *= static_cast<DataType>(2.0);

    BeginFastIndex2(phi, gdphi.bbox(), gdphi.data(), DataType);
    for (register int i = 0; i<nc; i++) {
      const int& n = idx[2*i]; const int& m = idx[2*i+1];
      normals[i](0) = -(FastIndex2(phi,n+ss(0),m)-FastIndex2(phi,n-ss(0),m))/dx2(0);
      normals[i](1) = -(FastIndex2(phi,n,m+ss(1))-FastIndex2(phi,n,m-ss(1)))/dx2(1);
      normals[i] /= normals[i].abs();
    }
    EndFastIndex2(phi);
  }

  virtual void CellDistances(GridHierarchy& GH, grid_data_type& gdphi, const int& nc, 
                             const int* idx, DataType* distance) {
    BeginFastIndex2(phi, gdphi.bbox(), gdphi.data(), DataType);
    for (register int n = 0; n<nc; n++)
      distance[n] = FastIndex2(phi,idx[2*n],idx[2*n+1]);
    EndFastIndex2(phi);
  }
};


template <class DataType>
class GFMGeom<DataType,3> : public GFMGeomBase<DataType,3> {
  typedef GFMGeomBase<DataType,3> base; 

public:
  typedef typename base::grid_data_type grid_data_type;
  typedef typename base::point_type point_type;

  GFMGeom() : base() {}

  virtual ~GFMGeom() {}

  virtual void CellNormals(GridHierarchy& GH, grid_data_type& gdphi, const int& nc, 
                           const int* idx, point_type* normals) {
    Coords ss(gdphi.stepsize());
    DCoords dx2(GH.worldStep(ss));  
    dx2 *= static_cast<DataType>(2.0);

    BeginFastIndex3(phi, gdphi.bbox(), gdphi.data(), DataType);
    for (register int i = 0; i<nc; i++) {
      const int& n = idx[3*i]; const int& m = idx[3*i+1]; const int& l = idx[3*i+2];
      normals[i](0) = -(FastIndex3(phi,n+ss(0),m,l)-FastIndex3(phi,n-ss(0),m,l))/dx2(0);
      normals[i](1) = -(FastIndex3(phi,n,m+ss(1),l)-FastIndex3(phi,n,m-ss(1),l))/dx2(1);
      normals[i](2) = -(FastIndex3(phi,n,m,l+ss(2))-FastIndex3(phi,n,m,l-ss(2)))/dx2(2);
      normals[i] /= normals[i].abs();
    }
    EndFastIndex3(phi);
  }

  virtual void CellDistances(GridHierarchy& GH, grid_data_type& gdphi, const int& nc, 
                             const int* idx, DataType* distance) {
    BeginFastIndex3(phi, gdphi.bbox(), gdphi.data(), DataType);
    for (register int n = 0; n<nc; n++)
      distance[n] = FastIndex3(phi,idx[3*n],idx[3*n+1],idx[3*n+2]);
    EndFastIndex3(phi);
  }
};

#endif

---