source: src/Box.cpp@ 014475

/*
 * 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 "Helpers/Assert.hpp"

Box::Box()
{
  M= new Matrix();
  M->one();
  Minv = new Matrix();
  Minv->one();
}

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

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.at(1,0)==_M.at(0,1),"Matrix used as cell_size was not symmetric");
  ASSERT(_M.at(2,0)==_M.at(0,2),"Matrix used as cell_size was not symmetric");
  ASSERT(_M.at(1,2)==_M.at(2,1),"Matrix used as cell_size was not symmetric");
  *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--;){
    double intpart,fracpart;
    fracpart = modf(helper.at(i),&intpart);
    if(fracpart<0.)
      fracpart+=1.;
    helper.at(i)=fracpart;
  }
  return translateIn(helper);
}

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

  // translate the Vector into each of the 27 neighbourhoods

  // first create 27 translation Vectors
  int N[NDIM];
  for (N[0]=-1;N[0]<=1;N[0]++){
    for (N[1]=-1;N[1]<=1;N[1]++){
      for (N[2]=-1;N[2]<=1;N[2]++){
        Vector translation = translateIn(Vector(N[0],N[1],N[2]));
        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{
  VECTORSET(std::list) expansion = explode(point1);
  double res = expansion.minDistSquared(point2);
  return res;
}

double Box::periodicDistance(const Vector &point1,const Vector &point2) const{
  double res;
  res = sqrt(periodicDistanceSquared(point1,point2));
  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);
  }
  return *this;
}

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