source: src/Box.cpp@ c0e28c

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 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 "CodePatterns/MemDebug.hpp"

#include "Box.hpp"

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

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Observer/Channels.hpp"
#include "CodePatterns/Observer/Notification.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "Shapes/BaseShapes.hpp"
#include "Shapes/ShapeOps.hpp"


Box::Box() :
    Observable("Box"),
    M(new RealSpaceMatrix()),
    Minv(new RealSpaceMatrix())
{
  internal_list.reserve(pow(3,3));
  coords.reserve(NDIM);
  index.reserve(NDIM);

  // observable stuff
  Channels *OurChannel = new Channels;
  NotificationChannels.insert( std::make_pair(this, OurChannel) );
  // add instance for each notification type
  for (size_t type = 0; type < NotificationType_MAX; ++type)
    OurChannel->addChannel(type);

  M->setIdentity();
  Minv->setIdentity();
  conditions.resize(3);
  conditions[0] = conditions[1] = conditions[2] = Wrap;
}

Box::Box(const Box& src) :
    Observable("Box"),
    conditions(src.conditions),
    M(new RealSpaceMatrix(*src.M)),
    Minv(new RealSpaceMatrix(*src.Minv))
{
  internal_list.reserve(pow(3,3));
  coords.reserve(NDIM);
  index.reserve(NDIM);

  // observable stuff
  Channels *OurChannel = new Channels;
  NotificationChannels.insert( std::make_pair(this, OurChannel) );
  // add instance for each notification type
  for (size_t type = 0; type < NotificationType_MAX; ++type)
    OurChannel->addChannel(type);
}

Box::Box(RealSpaceMatrix _M) :
    Observable("Box"),
    M(new RealSpaceMatrix(_M)),
    Minv(new RealSpaceMatrix())
{
  internal_list.reserve(pow(3,3));
  coords.reserve(NDIM);
  index.reserve(NDIM);

  // observable stuff
  Channels *OurChannel = new Channels;
  NotificationChannels.insert( std::make_pair(this, OurChannel) );
  // add instance for each notification type
  for (size_t type = 0; type < NotificationType_MAX; ++type)
    OurChannel->addChannel(type);

  ASSERT(M->determinant()!=0,"Matrix in Box construction was not invertible");
  *Minv = M->invert();
  conditions.resize(3);
  conditions[0] = conditions[1] = conditions[2] = Wrap;
}

Box::~Box()
{
  // observable stuff
  std::map<Observable *, Channels*>::iterator iter = NotificationChannels.find(this);
  delete iter->second;
  NotificationChannels.erase(iter);

  delete M;
  delete Minv;
}

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

void Box::setM(RealSpaceMatrix _M){
  ASSERT(_M.determinant()!=0,"Matrix in Box construction was not invertible");
  OBSERVE;
  NOTIFY(MatrixChanged);
  *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");
      break;
    }

  }
  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::vector) Box::explode(const Vector &point,int n) const{
  ASSERT(isInside(point),"Exploded point not inside Box");
  internal_explode(point, n);
  VECTORSET(std::vector) res(internal_list);
  return res;
}

void Box::internal_explode(const Vector &point,int n) const{
  internal_list.clear();
  size_t list_index = 0;

  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
  coords.clear();
  index.clear();
  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
  internal_list.resize(pow(2*n+1, dims));

  if(!dims){
    // all boundaries are ignored
    internal_list[list_index++] = point;
    return;
  }

  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");
          break;
        default:
          ASSERT(0,"No default case for this switch-case");
          break;
      }

    }
    // add back all ignored coordinates (not handled in above loop)
    helper+=mask;
    ASSERT(list_index < internal_list.size(),
        "Box::internal_explode() - we have estimated the number of vectors wrong: "
        +toString(list_index) +" >= "+toString(internal_list.size())+".");
    internal_list[list_index++] = 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];
    }
  }
}

VECTORSET(std::vector) 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(!isInside(point1) ? WrapPeriodically(point1) : point1);
  Vector helper2(!isInside(point2) ? WrapPeriodically(point2) : point2);
  internal_explode(helper1,1);
  double res = internal_list.minDistSquared(helper2);
  return res;
}

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

double Box::DistanceToBoundary(const Vector &point) const
{
  std::map<double, Plane> DistanceSet;
  std::vector<std::pair<Plane,Plane> > Boundaries = getBoundingPlanes();
  for (int i=0;i<NDIM;++i) {
    const double tempres1 = Boundaries[i].first.distance(point);
    const double tempres2 = Boundaries[i].second.distance(point);
    DistanceSet.insert( make_pair(tempres1, Boundaries[i].first) );
    LOG(1, "Inserting distance " << tempres1 << " and " << tempres2 << ".");
    DistanceSet.insert( make_pair(tempres2, Boundaries[i].second) );
  }
  ASSERT(!DistanceSet.empty(), "Box::DistanceToBoundary() - no distances in map!");
  return (DistanceSet.begin())->first;
}

Shape Box::getShape() const{
  return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
}

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

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

const std::vector<std::pair<Plane,Plane> >  Box::getBoundingPlanes() const
{
  std::vector<std::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));
  }
  ASSERT(res.size() == 3, "Box::getBoundingPlanes() - does not have three plane pairs!");
  return res;
}

void Box::setCuboid(const Vector &endpoint)
{
  OBSERVE;
  NOTIFY(MatrixChanged);
  ASSERT(endpoint[0]>0 && endpoint[1]>0 && endpoint[2]>0,"Vector does not define a full cuboid");
  M->setIdentity();
  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){
    OBSERVE;
    // new matrix
    NOTIFY(MatrixChanged);
    delete M;
    delete Minv;
    M    = new RealSpaceMatrix(*src.M);
    Minv = new RealSpaceMatrix(*src.Minv);
    // new boundary conditions
    NOTIFY(BoundaryConditionsChanged);
    conditions = src.conditions;
  }
  return *this;
}

Box &Box::operator=(const RealSpaceMatrix &mat)
{
  OBSERVE;
  NOTIFY(MatrixChanged);
  setM(mat);
  return *this;
}

std::ostream & operator << (std::ostream& ost, const Box &m)
{
  ost << m.getM();
  return ost;
}