source: src/Box.cpp@ f3be87

/*
 * Box.cpp
 *
 *  Created on: Jun 30, 2010
 *      Author: crueger
 */

//#include "Helpers/MemDebug.hpp"

#include "Box.hpp"

#include <cmath>

#include "Matrix.hpp"
#include "vector.hpp"
#include "Plane.hpp"

#include "Helpers/Assert.hpp"

Box::Box()
{
  M= new Matrix();
  M->one();
  Minv = new Matrix();
  Minv->one();
  conditions.resize(3);
  conditions[0] = conditions[1] = conditions[2] = Wrap;
}

Box::Box(const Box& src){
  M=new Matrix(*src.M);
  Minv = new Matrix(*src.Minv);
  conditions = src.conditions;
}

Box::~Box()
{
  delete M;
  delete Minv;
}

const Matrix &Box::getM() const{
  return *M;
}
const Matrix &Box::getMinv() const{
  return *Minv;
}

void Box::setM(Matrix _M){
  ASSERT(_M.determinant()!=0,"Matrix in Box construction was not invertible");
  *M    =_M;
  *Minv = M->invert();
}

Vector Box::translateIn(const Vector &point) const{
  return (*M) * point;
}

Vector Box::translateOut(const Vector &point) const{
  return (*Minv) * point;
}

Vector Box::WrapPeriodically(const Vector &point) const{
  Vector helper = translateOut(point);
  for(int i=NDIM;i--;){

    switch(conditions[i]){
    case Wrap:
      helper.at(i)=fmod(helper.at(i),1);
      helper.at(i)+=(helper.at(i)>0)?0:1;
      break;
    case Bounce:
      {
        // there probably is a better way to handle this...
        // all the fabs and fmod modf probably makes it very slow
        double intpart,fracpart;
        fracpart = modf(fabs(helper.at(i)),&intpart);
        helper.at(i) = fabs(fracpart-fmod(intpart,2));
      }
      break;
    case Ignore:
      break;
    default:
      ASSERT(0,"No default case for this");
    }

  }
  return translateIn(helper);
}

bool Box::isInside(const Vector &point) const
{
  bool result = true;
  Vector tester = translateOut(point);

  for(int i=0;i<NDIM;i++)
    result = result &&
              ((conditions[i] == Ignore) ||
               ((tester[i] >= -MYEPSILON) &&
               ((tester[i] - 1.) < MYEPSILON)));

  return result;
}


VECTORSET(std::list) Box::explode(const Vector &point,int n) const{
  VECTORSET(std::list) res;

  // translate the Vector into each of the 27 neighbourhoods

  // first create all translation Vectors
  // there are (n*2+1)^3 such vectors
  int max_dim = (n*2+1);
  int max_dim2 = max_dim*max_dim;
  int max = max_dim2*max_dim;
  // only one loop to avoid unneccessary jumps
  for(int i = 0;i<max;++i){
    // get all coordinates for this iteration
    int n1 = (i%max_dim)-n;
    int n2 = ((i/max_dim)%max_dim)-n;
    int n3 = ((i/max_dim2))-n;
    Vector translation = translateIn(Vector(n1,n2,n3));
    res.push_back(translation);
  }
  // translate all the translation vector by the offset defined by point
  res.translate(point);
  return res;
}

VECTORSET(std::list) Box::explode(const Vector &point) const{
  VECTORSET(std::list) res;

  // translate the Vector into each of the 27 neighbourhoods

  // first create all 27 translation Vectors
  // these loops depend on fixed parameters and can easily be expanded
  // by the compiler to allow code without jumps
  for(int n1 = -1;n1<=1;++n1){
    for(int n2 = -1;n2<=1;++n2){
      for(int n3 = -1;n3<=1;++n3){
        // get all coordinates for this iteration
        Vector translation = translateIn(Vector(n1,n2,n3));
        res.push_back(translation);
      }
    }
  }
  // translate all the translation vector by the offset defined by point
  res.translate(point);
  return res;
}

double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
  Vector helper1 = WrapPeriodically(point1);
  Vector helper2 = WrapPeriodically(point2);
  VECTORSET(std::list) expansion = explode(helper1);
  double res = expansion.minDistSquared(helper2);
  return res;
}

double Box::periodicDistance(const Vector &point1,const Vector &point2) const{
  double res;
  res = sqrt(periodicDistanceSquared(point1,point2));
  return res;
}

const Box::Conditions_t Box::getConditions(){
  return conditions;
}

void Box::setCondition(int i,Box::BoundaryCondition_t condition){
  conditions[i]=condition;
}

const vector<pair<Plane,Plane> >  Box::getBoundingPlanes(){
  vector<pair<Plane,Plane> > res;
  for(int i=0;i<NDIM;++i){
    Vector base1,base2,base3;
    base2[(i+1)%NDIM] = 1.;
    base3[(i+2)%NDIM] = 1.;
    Plane p1(translateIn(base1),
             translateIn(base2),
             translateIn(base3));
    Vector offset;
    offset[i]=1;
    Plane p2(translateIn(base1+offset),
             translateIn(base2+offset),
             translateIn(base3+offset));
    res.push_back(make_pair(p1,p2));
  }
  return res;
}

Box &Box::operator=(const Box &src){
  if(&src!=this){
    delete M;
    delete Minv;
    M    = new Matrix(*src.M);
    Minv = new Matrix(*src.Minv);
    conditions = src.conditions;
  }
  return *this;
}

Box &Box::operator=(const Matrix &mat){
  setM(mat);
  return *this;
}