source: src/Box.cpp@ a42fee

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * 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 <cstdlib>
#include <iostream>
#include <sstream>

#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(static_cast<Observable *>(this), OurChannel) );
  // add instance for each notification type
  for (size_t type = 0; type < NotificationType_MAX; ++type)
    OurChannel->addChannel(type);

  M->setIdentity();
  Minv->setIdentity();
}

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(static_cast<Observable *>(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(static_cast<Observable *>(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();
}

Box::~Box()
{
  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;
  if (_M != *M)
    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::enforceBoundaryConditions(const Vector &point) const{
  Vector helper = translateOut(point);
  for(int i=NDIM;i--;){

    switch(conditions[i]){
    case BoundaryConditions::Wrap:
      helper.at(i)=fmod(helper.at(i),1);
      helper.at(i)+=(helper.at(i)>=0)?0:1;
      break;
    case BoundaryConditions::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 BoundaryConditions::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] == BoundaryConditions::Ignore) ||
               ((tester[i] >= -MYEPSILON) &&
               ((tester[i] - 1.) < MYEPSILON)));

  return result;
}

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

  for(int i=0;i<NDIM;i++)
    result = result &&
              ((conditions[i] != BoundaryConditions::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]==BoundaryConditions::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

  {
    const size_t new_size = pow(2*n+1, dims);
    if (internal_list.size() != new_size)
      internal_list.resize(new_size);
  }

  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 BoundaryConditions::Wrap:
          helper[coords[i]] = index[i]+translater[coords[i]];
          break;
        case BoundaryConditions::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 BoundaryConditions::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);
}

const Vector Box::periodicDistanceVector(const Vector &point1,const Vector &point2) const{
  Vector helper1(enforceBoundaryConditions(point1));
  Vector helper2(enforceBoundaryConditions(point2));
  internal_explode(helper1,1);
  const Vector res = internal_list.minDistance(helper2);
  return res;
}

double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
  const Vector res = periodicDistanceVector(point1, point2);
  return res.NormSquared();
}

double Box::periodicDistance(const Vector &point1,const Vector &point2) const{
  double 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 std::string Box::getConditionNames() const
{
  std::stringstream outputstream;
  outputstream << conditions;
  return outputstream.str();
}

const BoundaryConditions::Conditions_t & Box::getConditions() const
{
  return conditions.get();
}

const BoundaryConditions::BoundaryCondition_t Box::getCondition(size_t i) const
{
  return conditions.get(i);
}

void Box::setCondition(size_t i, const BoundaryConditions::BoundaryCondition_t _condition)
{
  OBSERVE;
  if (conditions.get(i) != _condition)
    NOTIFY(BoundaryConditionsChanged);
  conditions.set(i, _condition);
}

void Box::setConditions(const BoundaryConditions::Conditions_t & _conditions)
{
  OBSERVE;
  if (conditions.get() != _conditions)
  NOTIFY(BoundaryConditionsChanged);
  conditions.set(_conditions);
}

void Box::setConditions(const std::string & _conditions)
{
  OBSERVE;
  NOTIFY(BoundaryConditionsChanged);
  std::stringstream inputstream(_conditions);
  inputstream >> conditions;
}

void Box::setConditions(const std::vector< std::string >& _conditions)
{
  OBSERVE;
  NOTIFY(BoundaryConditionsChanged);
  conditions.set(_conditions);
}

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