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_;
};