source: src/Box.cpp@ bcf653

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

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

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "Helpers/MemDebug.hpp"

#include "Box.hpp"

#include <cmath>
#include <iostream>
#include <cstdlib>

#include "LinearAlgebra/Matrix.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/Plane.hpp"

#include "Helpers/Assert.hpp"

using namespace std;

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{
  ASSERT(isInside(point),"Exploded point not inside Box");
  VECTORSET(std::list) res;

  Vector translater = translateOut(point);
  Vector mask; // contains the ignored coordinates

  // count the number of coordinates we need to do
  int dims = 0; // number of dimensions that are not ignored
  vector<int> coords;
  vector<int> index;
  for(int i=0;i<NDIM;++i){
    if(conditions[i]==Ignore){
      mask[i]=translater[i];
      continue;
    }
    coords.push_back(i);
    index.push_back(-n);
    dims++;
  } // there are max vectors in total we need to create

  if(!dims){
    // all boundaries are ignored
    res.push_back(point);
    return res;
  }

  bool done = false;
  while(!done){
    // create this vector
    Vector helper;
    for(int i=0;i<dims;++i){
      switch(conditions[coords[i]]){
        case Wrap:
          helper[coords[i]] = index[i]+translater[coords[i]];
          break;
        case Bounce:
          {
            // Bouncing the coordinate x produces the series:
            // 0 -> x
            // 1 -> 2-x
            // 2 -> 2+x
            // 3 -> 4-x
            // 4 -> 4+x
            // the first number is the next bigger even number (n+n%2)
            // the next number is the value with alternating sign (x-2*(n%2)*x)
            // the negative numbers produce the same sequence reversed and shifted
            int n = abs(index[i]) + ((index[i]<0)?-1:0);
            int sign = (index[i]<0)?-1:+1;
            int even = n%2;
            helper[coords[i]]=n+even+translater[coords[i]]-2*even*translater[coords[i]];
            helper[coords[i]]*=sign;
          }
          break;
        case Ignore:
          ASSERT(0,"Ignored coordinate handled in generation loop");
        default:
          ASSERT(0,"No default case for this switch-case");
      }

    }
    // add back all ignored coordinates (not handled in above loop)
    helper+=mask;
    res.push_back(translateIn(helper));
    // set the new indexes
    int pos=0;
    ++index[pos];
    while(index[pos]>n){
      index[pos++]=-n;
      if(pos>=dims) { // it's trying to increase one beyond array... all vectors generated
        done = true;
        break;
      }
      ++index[pos];
    }
  }
  return res;
}

VECTORSET(std::list) Box::explode(const Vector &point) const{
  ASSERT(isInside(point),"Exploded point not inside Box");
  return explode(point,1);
}

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

void Box::setCuboid(const Vector &endpoint){
  ASSERT(endpoint[0]>0 && endpoint[1]>0 && endpoint[2]>0,"Vector does not define a full cuboid");
  M->one();
  M->diagonal()=endpoint;
  Vector &dinv = Minv->diagonal();
  for(int i=NDIM;i--;)
    dinv[i]=1/endpoint[i];
}

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