hexahedron.h

#pragma once
 
#include "base/src/cube.h"
#include "base/src/vect.h"
 
// Assumptions: 1) Faces are planar; 2) face normals are oriented inward; Hexahedron is convex
// Class used to create a grid. This grid is generated by a user and used to import data (see function readDepletionGridFile
// Convex hexahedron. Expecting points to be added in the following order: clockwise from top left closest corner (top-southwest corner)
class Hexahedron
{
public:
  Hexahedron():ext_(DCube::nothing()),vSum_(0.0) {}
  void addPoint(const DVect3 &v) 
  { 
    points_.push_back(v); 
    ext_ = ext_.expandToInclude(v);
    vSum_ +=v;
  }
  const std::vector<DVect3> &points() const { return points_; }
  const DCube               &extent() const { return ext_; }
 
  bool isInConvex(DVect3 const& point) const
  {
    static const int faceGpIndex[6][4]  = {{4,7,3,0} ,  /* Hexahefron Face front */
                                           {5,1,2,6} ,  /*                 back */
                                           {4,0,1,5} ,  /*                 left */
                                           {7,6,2,3} ,  /*                 rigth */
                                           {0,3,2,1} ,  /*                 top */
                                           {4,5,6,7}};  /*                 bottom */
 
    assert(points_.size() == 8);
 
    if (ext_.isIn(point)==false) return false;
 
    for (int i = 0; i < 6; ++i) 
    {
      const DVect3 &v0 = points_[faceGpIndex[i][0]];
      const DVect3 &v1 = points_[faceGpIndex[i][1]];
      const DVect3 &v2 = points_[faceGpIndex[i][2]];
 
      DVect3 point2face = v0 - point;         // use an arbitrary point on face
      // calculate face normal vector
      const DVect3 &dir1 = v1 - v0;
      const DVect3 &dir2 = v2 - v0;
      const DVect3 &nV   = (dir1&dir2).unit(); // cross prodcut
 
      const double d = point2face|nV;        // dot product of the normal face vector, face distance to the point 
      /*
      d /= point2face.mag();           // for numeric stability
      //static const double bound = -1e-16; // use 1e16 to exclude boundaries
      if (d < bound)
      */
      if (d < 1e-16) // negative dot prodcut  indicates point outside 
        return false;
    }
    return true;
  }
 
  bool isRegular() const
  {
    assert(points_.size()==8);
    DVect3 v1 = (points_[0]-points_[4]);//.unit();
    DVect3 v2 = (points_[0]-points_[1]);//.unit();
    DVect3 v3 = (points_[0]-points_[3]);//.unit();
 
    double dotP  = v1|v2;
    if (fabs(dotP) > 1e-14) return false; 
    dotP = v1|v3;
    if (fabs(dotP) > 1e-14) return false; 
    dotP = v2|v3;
    if (fabs(dotP) > 1e-14) return false; 
    return true;
  }
 
  DVect3 geometricCenter() const 
  {
    int nPoints = to<int>(points_.size());
    if (nPoints == 0) return DVect3(0.0);
    return vSum_/nPoints;
  }
 
  DVect3 maxEdges()
  {
    assert(points_.size() == 8);
    DVect3 maxEdges(ext_.width(),ext_.height(),ext_.depth());
    return maxEdges;
  }
 
private:
  std::vector<DVect3> points_;
  DCube ext_;
  DVect3 vSum_;
};