Geometry.H

Go to the documentation of this file.

 
#ifndef _GEOMETRY_H_
#define _GEOMETRY_H_
 
#include <algorithm>
#include <cmath>
#include <vector>
#include "CRPacket.H"
#include "../IO/ParmParser.H"
#include "../Utils/RealBox.H"
#include "../Utils/Types.H"
#include "../Utils/Vec3.H"
 
namespace criptic {
 
  class Geometry {
 
  public:
 
    typedef enum {
      openBC = 0,
      absorbingBC,
      reflectingBC,
      periodicBC
#ifdef TRACK_PITCH_ANGLE
      , reflectingPitchAngleBC
#endif
    } bcType;
 
    Geometry(const ParmParser &pp);
 
    const RealBox& xBox() const { return box; }
 
    const RealVec& xLo() const { return box.lo; }
 
    const RealVec& xHi() const { return box.hi; }
 
    const RealVec& xSize() const { return probSize; }
 
    const Vec3<bcType>& bLo() const { return bcLo; }
    
    const Vec3<bcType>& bHi() const { return bcHi; }
 
    const BoolVec& isPeriodic() const { return periodicity; }
 
    bool isAnyPeriodic() const { return hasPeriodicDim; }
 
    bool isAnyOpen() const {
      BoolVec lo = bcLo == openBC;
      BoolVec hi = bcHi == openBC;
      return max(lo, hi).max();
    }
 
    void applyBC(RealVec& x
#ifdef TRACK_PITCH_ANGLE
      , CRPacket& pkt
#endif
    ) const {
      for (int i = 0; i < 3; i++) {
        if (periodicity[i]) {
          while (x[i] < box.lo[i]) x[i] += probSize[i];
          while (x[i] > box.hi[i]) x[i] += probSize[i];
        } else if (bcLo[i] == reflectingBC && x[i] < box.lo[i]) {
          x[i] += 2 * (box.lo[i] - x[i]);
        } else if (bcHi[i] == reflectingBC && x[i] > box.hi[i]) {
          x[i] -= 2 * (x[i] - box.hi[i]);
        }
#ifdef TRACK_PITCH_ANGLE
        else if (bcLo[i] == reflectingPitchAngleBC && x[i] < box.lo[i]) {
          x[i] += 2 * (box.lo[i] - x[i]);
          pkt.mu = -pkt.mu;
        } else if (bcHi[i] == reflectingPitchAngleBC && x[i] > box.hi[i]) {
          x[i] -= 2 * (x[i] - box.hi[i]);
          pkt.mu = -pkt.mu;
        }
#endif
      }
    }
 
    Real periodicLoop(const Real x,
              const int i,
              const bool recenter = false) const {
      Real x_ = x;
      if (periodicity[i]) {
    while (x_ < box.lo[i]) x_ += probSize[i];
    while (x_ > box.hi[i]) x_ -= probSize[i];
    if (recenter) x_ -= (box.lo[i] + box.hi[i]) / 2.0;
      }
      return x_;
    }
 
    RealVec periodicLoop(const RealVec& x,
             const bool recenter = false) const {
      if (!hasPeriodicDim) return x;
      RealVec x_ = x;
      for (int i=0; i<3; i++) x_[i] = periodicLoop(x_[i], i, recenter);
      return x_;
    }
 
    std::vector<RealBox> periodicLoopBox(const RealBox& b,
                     const bool recenter = false) const {
      std::vector<RealBox> vrb;
      RealBox rb;
      rb.lo = periodicLoop(b.lo, recenter);
      rb.hi = periodicLoop(b.hi, recenter);
      for (int i=0; i<3; i++) {
    if (rb.lo[i] > rb.hi[i]) {
      RealBox rb1 = rb;
      RealBox rb2 = rb;
      if (recenter) {
        rb1.hi[i] = probSize[i]/2;
        rb2.lo[i] = -probSize[i]/2;
      } else {
        rb1.hi[i] = box.hi[i];
        rb2.lo[i] = box.lo[i];
      }
      std::vector<RealBox> vrb1 = periodicLoopBox(rb1, recenter);
      std::vector<RealBox> vrb2 = periodicLoopBox(rb2, recenter);
      for (auto r : vrb1) vrb.push_back(r);
      for (auto r : vrb2) vrb.push_back(r);
    }
      }
      if (vrb.size() == 0) vrb.push_back(rb);
      return vrb;
    }
    
    bool contains(const RealVec& x) const {
      for (int i=0; i<3; i++) {
    if (!periodicity[i] && (x[i] < box.lo[i] || x[i] > box.hi[i]))
      return false;
      }
      return true;
    }
 
    RealVec disp(const RealVec& x1, const RealVec& x2) const {
      return periodicLoop(x1 - x2, true);
    }
    
    Real dist(const RealVec& x1, const RealVec& x2) const {
      return disp(x1, x2).mag();
    }
    
    Real minBoxDist(const RealVec& x,
            const RealVec& boxLo,
            const RealVec& boxHi) const {
      Real d2 = 0.0;
      for (int i=0; i<3; i++) {
    if (boxLo[i] <= x[i] && x[i] <= boxHi[i]) continue;
    Real dxLo = x[i] - boxLo[i];
    if (periodicity[i]) {
      while (dxLo < 0) dxLo += probSize[i];
      while (dxLo > probSize[i]) dxLo -= probSize[i];
    }
    Real dxHi = x[i] - boxHi[i];
    if (periodicity[i]) {
      while (dxHi < 0) dxHi += probSize[i];
      while (dxHi > probSize[i]) dxHi -= probSize[i];
    }
    Real dx2Lo = dxLo * dxLo;
    Real dx2Hi = dxHi * dxHi;
    d2 += std::max(dx2Lo, dx2Hi);
      }
      return std::sqrt(d2);
    }
 
    Real maxBoxDist(const RealVec& x,
            const RealVec& boxLo,
            const RealVec& boxHi) const {
      Real d2 = 0.0;
      for (int i=0; i<3; i++) {
    Real dxLo = x[i] - boxLo[i];
    if (periodicity[i]) {
      while (dxLo < 0) dxLo += probSize[i];
      while (dxLo > probSize[i]) dxLo -= probSize[i];
    }
    Real dxHi = x[i] - boxHi[i];
    if (periodicity[i]) {
      while (dxHi < 0) dxHi += probSize[i];
      while (dxHi > probSize[i]) dxHi -= probSize[i];
    }
    Real dx2Lo = dxLo * dxLo;
    Real dx2Hi = dxHi * dxHi;
    d2 += std::max(dx2Lo, dx2Hi);
      }
      return std::sqrt(d2);
    }
 
    void minMaxBoxDist2(const RealVec& x,
            const RealVec& boxLo,
            const RealVec& boxHi,
            RealVec& d2Min,
            RealVec& d2Max) const {
      d2Min = 0.0;
      d2Max = 0.0;
      for (int i=0; i<3; i++) {
    Real dxLo = x[i] - boxLo[i];
    if (periodicity[i]) {
      while (dxLo < 0) dxLo += probSize[i];
      while (dxLo > probSize[i]) dxLo -= probSize[i];
    }
    Real dxHi = x[i] - boxHi[i];
    if (periodicity[i]) {
      while (dxHi < 0) dxHi += probSize[i];
      while (dxHi > probSize[i]) dxHi -= probSize[i];
    }
    Real dx2Lo = dxLo * dxLo;
    Real dx2Hi = dxHi * dxHi;
    d2Max[i] = std::max(dx2Lo, dx2Hi);
    if (x[i] < boxLo[i] || x[i] > boxHi[i])
      d2Min[i] = std::min(dx2Lo, dx2Hi);
    else
      d2Min[i] = 0.0;
      }
    }
 
    void minMaxBoxDist2(const RealVec& x,
            const RealVec& boxLo,
            const RealVec& boxHi,
            Real& d2Min,
            Real& d2Max) const {
      RealVec d2MinVec, d2MaxVec;
      minMaxBoxDist2(x, boxLo, boxHi, d2MinVec, d2MaxVec);
      d2Min = d2MinVec.sum();
      d2Max = d2MaxVec.sum();
    }
    
    void minMaxBoxDist(const RealVec& x,
               const RealVec& boxLo,
               const RealVec& boxHi,
               Real& dMin,
               Real& dMax) const {
      Real d2Min, d2Max;
      minMaxBoxDist2(x, boxLo, boxHi, d2Min, d2Max);
      dMin = std::sqrt(d2Min);
      dMax = std::sqrt(d2Max);
    }
    
  private:
 
    Vec3<bcType> bcLo;         
    Vec3<bcType> bcHi;         
    RealBox box;               
    RealVec probSize;          
    BoolVec periodicity;       
    bool hasPeriodicDim;       
  };
 
}
 
#endif
// _GEOMETRY_H_