direct_product_mesh.h

Go to the documentation of this file.

/*
This file is part of the Ristra Wonton project.
Please see the license file at the root of this repository, or at:
    https://github.com/laristra/wonton/blob/master/LICENSE
*/

#ifndef WONTON_DIRECT_PRODUCT_MESH_H_
#define WONTON_DIRECT_PRODUCT_MESH_H_

#include <array>
#include <cassert>
#include <vector>

#include "wonton/support/wonton.h"
#include "wonton/support/CoordinateSystem.h"
#include "wonton/support/structured_partitioner.h"

namespace Wonton {

template<int D, class CoordSys=Wonton::DefaultCoordSys>
class Direct_Product_Mesh {

 public:

  // ==========================================================================
  // Constructors and destructors

  Direct_Product_Mesh() = delete;

  Direct_Product_Mesh(const std::array<std::vector<double>,D> &axis_points_in,
                      const Wonton::Executor_type *executor = nullptr,
                      const int num_ghost_layers = 0,
                      const int dpart = D);

  Direct_Product_Mesh & operator=(
      const Direct_Product_Mesh<D,CoordSys> &) = delete;

  // ==========================================================================
  // Accessors

  int space_dimension() const;

  bool distributed() const;
  
  int num_axis_points(const int dim,
                      const Entity_type etype = Entity_type::ALL) const;

  int num_ghost_layers() const;

  Entity_type axis_point_type(const int dim, const int pointid) const;
  
  GID_t point_global_index(const int dim, const int pointid) const;
  
  int num_axis_cells(const int dim,
                     const Entity_type etype = Entity_type::ALL) const;

  double get_axis_point(const int dim, const int pointid) const;

  bool point_on_external_boundary(const int dim, const int pointid) const;

  void get_global_coord_bounds(const int dim, double *lo, double *hi) const;
    
  void get_local_coord_bounds(const int dim, double *lo, double *hi) const;
  
 private:

  // ==========================================================================
  // Class data

  std::array<std::array<double,2>,D> coord_bounds_global_;

  std::array<GID_t,D> num_axis_points_global_;

  int num_ghost_layers_;

  std::array<std::array<GID_t,2>,D> gid_bounds_local_;

  std::array<std::vector<double>,D> axis_points_;

  bool distributed_=false;

#ifdef WONTON_ENABLE_MPI
  MPI_Comm mycomm_;
  int nprocs_=1;
  int rank_=0;
#endif
  
};  // class Direct_Product_Mesh


// ============================================================================
// Constructors and destructors

// ____________________________________________________________________________
// Constructor

template<int D, class CoordSys>
Direct_Product_Mesh<D,CoordSys>::Direct_Product_Mesh(
    const std::array<std::vector<double>,D> &axis_points_in,
    Wonton::Executor_type const *executor, int num_ghost_layers, int dpart) {

  for (int d = 0; d < D; d++) {
    coord_bounds_global_[d][0] = axis_points_in[d][0];
    coord_bounds_global_[d][1] = axis_points_in[d][axis_points_in[d].size()-1];
    num_axis_points_global_[d] = axis_points_in[d].size();
  }
  
#ifdef WONTON_ENABLE_MPI
  
  mycomm_ = MPI_COMM_NULL;
  auto mpiexec = dynamic_cast<MPIExecutor_type const *>(executor);
  if (mpiexec && mpiexec->mpicomm != MPI_COMM_NULL) {
    mycomm_ = mpiexec->mpicomm;
    MPI_Comm_size(mycomm_, &nprocs_);
    MPI_Comm_rank(mycomm_, &rank_);
    if (nprocs_ > 1) distributed_ = true;
  }
  if (distributed_) {
    // compute partition of the global domain on each processor

    int seed = 42;  // seed for randomization to ensure we get
                    // consistent results on all processors

    std::array<GID_t, D> ncells;
    for (int d = 0; d < D; d++) ncells[d] = axis_points_in[d].size()-1;
    auto partition_bounds =
        structured_partitioner<D>(nprocs_, ncells, dpart, seed);
    
    // integer bounds of points of OWNED cells in each direction on
    // local partition. The 'partition_bounds' are returned as
    // ((xlo,ylo,zlo),(xhi,yhi,zhi)) while bounds in this class are
    // stored as ((xlo,xhi),(ylo,yhi),(zlo,zhi))
    
    for (int d = 0; d < D; d++) {
      gid_bounds_local_[d][0] = partition_bounds[rank_][0][d];
      gid_bounds_local_[d][1] = partition_bounds[rank_][1][d];
    }
    
    // Build list of axis points including ghost points
    
    num_ghost_layers_ = num_ghost_layers;
    int NG = num_ghost_layers_;

    double inf  =  std::numeric_limits<double>::infinity();
    for (int d = 0; d < D; d++) {
      GID_t lo = gid_bounds_local_[d][0];
      GID_t hi = gid_bounds_local_[d][1];
      
      // owned + ghost points
      int npall = static_cast<int>(hi-lo+1+2*NG);
      axis_points_[d].resize(npall);

      if (lo == 0)
        for (int k = 0; k < NG; k++)
          axis_points_[d][k] = -inf;
      else
        for (int k = 0; k < NG; k++)
          axis_points_[d][k] = (lo-NG+k < 0) ? -inf :
              axis_points_[d][k] = axis_points_in[d][lo-NG+k];
        
      std::copy(axis_points_in[d].begin()+lo, axis_points_in[d].begin()+hi+1,
                axis_points_[d].begin()+NG);

      if (hi == num_axis_points_global_[d]-1)
        for (int k = 0; k < NG; k++)
          axis_points_[d][npall-NG+k] = inf;
      else
        for (int k = 0; k < NG; k++)
          axis_points_[d][npall-NG+k] = (hi+1+k >= num_axis_points_global_[d]) ?
              inf : axis_points_in[d][hi+1+k];
    }
  } else {
    num_ghost_layers_ = 0;
    for (int d = 0; d < D; ++d)
      axis_points_[d] = axis_points_in[d];
    for (int d = 0; d < D; ++d) {
      gid_bounds_local_[d][0] = 0;
      gid_bounds_local_[d][1] = num_axis_points_global_[d];
    }
  }
    
#else
  num_ghost_layers_ = 0;
  for (int d = 0; d < D; ++d)
    axis_points_[d] = axis_points_in[d];
  for (int d = 0; d < D; ++d) {
    gid_bounds_local_[d][0] = 0;
    gid_bounds_local_[d][1] = num_axis_points_global_[d];
  }
#endif

}


// ============================================================================
// Accessors

// ____________________________________________________________________________
// Get the dimensionality of the mesh.
template<int D, class CoordSys>
int Direct_Product_Mesh<D,CoordSys>::space_dimension() const {
  return(D);
}

// ____________________________________________________________________________
// Is the mesh distributed?
template<int D, class CoordSys>
bool Direct_Product_Mesh<D,CoordSys>::distributed() const {
  return(distributed_);
}

// ____________________________________________________________________________
// Get the number of axis points (cell boundary coordinates) of particular type
template<int D, class CoordSys>
int Direct_Product_Mesh<D,CoordSys>::num_axis_points(const int dim,
                                                     Entity_type etype) const {
  assert(dim >= 0);
  assert(dim < D);
  switch(etype) {
    case PARALLEL_OWNED: {
      // Eliminate the guaranteed ghost points (i.e. lower points of
      // ghost cells before the owned cells and the upper points of
      // ghost cells after the owned cells). For instance, if there
      // are 2 ghost layers, then 4 points are guaranteed to be
      // ghosts as shown here (g is ghost, o is owned):
      //
      // axis pts  g     g    g/o    o     o     o     g     g
      //           *-----*-----*-----*-----*-----*-----*-----*
      // cells        g     g     o     o     o     o     o
      
      int n = axis_points_[dim].size()-2*num_ghost_layers_;
      
      // the upper coordinate of the partition in any direction is
      // ALWAYS owned by the current partition while the lower
      // coordinate of the partition in any direction is ghost (owned
      // by another partition) _unless_ it is on the outer low
      // boundary of global domain
      
      if (gid_bounds_local_[dim][0] != 0)
        n--;
      return n;
    }
    case PARALLEL_GHOST: {
      int n = 2*num_ghost_layers_;
      if (gid_bounds_local_[dim][0] != 0)  // see figure above
        n++;
      return n;
    }
    case ALL:
      return axis_points_[dim].size();
    default:
      return 0;
  }
}


// ____________________________________________________________________________
// Get the number of ghost layers on this mesh
template<int D, class CoordSys>
int Direct_Product_Mesh<D,CoordSys>::num_ghost_layers() const {
  return num_ghost_layers_;
}

// ____________________________________________________________________________
// Get the specified axis point (cell boundary coordinate).
template<int D, class CoordSys>
double Direct_Product_Mesh<D,CoordSys>::get_axis_point(
    const int dim, const int pointid) const {
  assert(dim >= 0);
  assert(dim < D);

#if DEBUG
  int npall = axis_points_[dim].size();
  assert(pointid >= -num_ghost_layers_);
  assert(pointid <= npall);
#endif

  // pointid starts at -num_ghost_layers for lowest ghost point  
  return axis_points_[dim][pointid+num_ghost_layers_];
}

// ____________________________________________________________________________
// is point on external (global) boundary along given axis
// Both the owned points on the global boundary AND their boundary
// ghost counterparts will return true in this call
template<int D, class CoordSys>
bool Direct_Product_Mesh<D,CoordSys>::point_on_external_boundary(
    const int dim, const int pointid) const {
  assert(dim >= 0);
  assert(dim < D);

  GID_t global_point_index = point_global_index(dim, pointid);
  if (global_point_index <= 0)
    return true;               // coord on or outside of lower global bounds
  else return global_point_index >= num_axis_points_global_[dim] - 1;
}


// ____________________________________________________________________________
// Entity type of a point along the axis
template<int D, class CoordSys>
Entity_type Direct_Product_Mesh<D,CoordSys>::axis_point_type(
    const int dim, const int pointid) const {
  int npall = num_axis_points(dim, ALL);
  if (point_on_external_boundary(dim, pointid)) {  // global boundary
    if (pointid < 0 || pointid > npall-2*num_ghost_layers_-1)
      return BOUNDARY_GHOST;
    else if (pointid == 0)
      return PARALLEL_OWNED;  // lo bndry of owned cells in partition
    else if (pointid == npall-2*num_ghost_layers_-1)
      return PARALLEL_GHOST;  // hi bndry of owned cells in partition
  } else {
    if (pointid < 0 || pointid >= npall-2*num_ghost_layers_-1)
      return PARALLEL_GHOST;
    else
      return PARALLEL_OWNED;
  }
  return TYPE_UNKNOWN;
}
  
// ____________________________________________________________________________
// get global index of point along given axis

// Note: we expect the number of points along any one axis to be well
// within the range representable by an int
template<int D, class CoordSys>
GID_t Direct_Product_Mesh<D,CoordSys>::point_global_index(
    const int dim, const int pointid) const {
  assert(dim >= 0);
  assert(dim < D);

  return gid_bounds_local_[dim][0] + pointid;
}

// ____________________________________________________________________________
// number of cells along given axis
template<int D, class CoordSys>
int Direct_Product_Mesh<D,CoordSys>::num_axis_cells(
    const int dim, const Entity_type etype) const {
  assert(dim >= 0);
  assert(dim < D);
  switch (etype) {
    case ALL: 
      return num_axis_points(dim, ALL) - 1;
    case PARALLEL_OWNED:
      return num_axis_points(dim, ALL) - 1 - 2*num_ghost_layers_;
    case PARALLEL_GHOST:
      return 2*num_ghost_layers_;
    default:
      return 0;
  }
}


// ____________________________________________________________________________
// Coordinate bounds of global domain
template<int D, class CoordSys>
void Direct_Product_Mesh<D,CoordSys>::get_global_coord_bounds(const int dim,
                                                              double *plo,
                                                              double *phi) const {
    *plo = coord_bounds_global_[dim][0];
    *phi = coord_bounds_global_[dim][1];
}
    
// ____________________________________________________________________________
// Coordinate bounds of owned cells in local domain (ghost layers EXCLUDED)
template<int D, class CoordSys>
void Direct_Product_Mesh<D,CoordSys>::get_local_coord_bounds(const int dim,
                                                             double *plo,
                                                             double *phi) const {
  int npall = axis_points_[dim].size();
  *plo = axis_points_[dim][num_ghost_layers_];
  *phi = axis_points_[dim][npall-1-num_ghost_layers_];
}
    

}  // namespace Wonton

#endif  // WONTON_DIRECT_PRODUCT_MESH_H_