/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. 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/>.
 */

/*
 * CubeMesh.cpp
 *
 *  Created on: Jan 20, 2012
 *      Author: heber
 */

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

//#include "CodePatterns/MemDebug.hpp"

#include <boost/lambda/lambda.hpp>
#include <cmath>

#include "CodePatterns/Log.hpp"
#include "CubeMesh.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/Vector.hpp"

/** Constructor for class CubeMesh.
 *
 * Here, we generate nodes homogeneously distributed over a cuboid.
 *
 * \a offset shifts the coordinates, e.g. if counts = 2, we would set nodes at
 * 0 and 0.5 with offset = 0 and 0.4999... and 0.9999.. with offset = 0.999...
 *
 * @param counts number of points per axis
 * @param offset offset Vector (coordinates in [0,1))
 * @param M matrix to transform the default cuboid from (0,0,0) to (1,1,1).
 */
CubeMesh::CubeMesh(const Vector &counts, const Vector &offset, const RealSpaceMatrix &M)
{
#ifndef NDEBUG
  for (size_t i=0;i<NDIM;++i) {
    ASSERT((offset[i] >= 0.) && (offset[i] < 1.),
        "CubeMesh::CubeMesh() - offset coordinates must be in [0,1) but offset["
        +toString(i)+"] is "+toString(offset[i])+".");
    ASSERT(counts[i] == (int)counts[i],
        "CubeMesh::CubeMesh() - counts["+toString(i)+"] must be integer: != "+toString(counts[i])+".");
  }
#endif

  std::vector< unsigned int> size_counts;
  for (size_t i=0;i<NDIM;++i) {
    size_counts.push_back( (int)counts[i] );
  }
  init(size_counts, offset, M);
}

/** Constructor for class CubeMesh.
 *
 * Here, we generate nodes homogeneously distributed over a cuboid.
 *
 * \a offset shifts the coordinates, e.g. if counts = 2, we would set nodes at
 * 0 and 0.5 with offset = 0 and 0.4999... and 0.9999.. with offset = 0.999...
 *
 * @param counts number of points per axis
 * @param offset offset Vector (coordinates in [0,1))
 * @param M matrix to transform the default cuboid from (0,0,0) to (1,1,1).
 */
CubeMesh::CubeMesh(const std::vector< unsigned int > &counts, const Vector &offset, const RealSpaceMatrix &M)
{
  ASSERT(counts.size() == 3,
      "CubeMesh::CubeMesh() - counts does not have three but "
      +toString(counts.size())+" entries.");

  init(counts, offset, M);
}

void CubeMesh::init(const std::vector< unsigned int > &counts, const Vector &offset, const RealSpaceMatrix &M)
{
  RealSpaceMatrix partition;

  partition.setZero();
  for(int i=0;i<NDIM;i++)
    partition.at(i,i) =  1./(double)counts[i];
  LOG(1, "INFO: partition is " << partition << ".");

  // go over [0,1]^3 filler grid
  size_t n[NDIM];
  for (n[0] = 0; n[0] < counts[0]; n[0]++)
    for (n[1] = 0; n[1] < counts[1]; n[1]++)
      for (n[2] = 0; n[2] < counts[2]; n[2]++) {
        const Vector CurrentPosition((double)n[0], (double)n[1], (double)n[2]);
        nodes.push_back(M * partition * (CurrentPosition + offset));
      }
}

/** Destructor for class CubeMesh.
 *
 */
CubeMesh::~CubeMesh()
{}