source: src/Actions/FillAction/FillRegularGridAction.cpp@ 345eda

/*
 * 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.
 */

/*
 * FillRegularGridAction.cpp
 *
 *  Created on: Jan 12, 2012
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "Atom/atom.hpp"
#include "Atom/CopyAtoms/CopyAtoms_withBonds.hpp"
#include "CodePatterns/Log.hpp"
#include "Descriptors/MoleculeOrderDescriptor.hpp"
#include "Filling/Cluster.hpp"
#include "Filling/Filler.hpp"
#include "Filling/Mesh/CubeMesh.hpp"
#include "Filling/Predicates/IsInsideSurface_FillPredicate.hpp"
#include "Filling/Predicates/IsVoidNode_FillPredicate.hpp"
#include "Filling/Predicates/Ops_FillPredicate.hpp"
#include "LinkedCell/linkedcell.hpp"
#include "LinkedCell/PointCloudAdaptor.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "Parser/FormatParserInterface.hpp"
#include "Parser/FormatParserStorage.hpp"
#include "Shapes/BaseShapes.hpp"
#include "Tesselation/tesselation.hpp"
#include "Tesselation/BoundaryLineSet.hpp"
#include "Tesselation/BoundaryTriangleSet.hpp"
#include "Tesselation/CandidateForTesselation.hpp"
#include "World.hpp"


#include <algorithm>
#include <iostream>
#include <string>
#include <vector>

#include "Actions/FillAction/FillRegularGridAction.hpp"

using namespace MoleCuilder;

// and construct the stuff
#include "FillRegularGridAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */
Action::state_ptr FillRegularGridAction::performCall() {
  typedef std::vector<atom*> AtomVector;

  // obtain information
  getParametersfromValueStorage();

  // get all present atoms for UndoState
  AtomVector presentatoms = World::getInstance().getAllAtoms();

  // get selected atoms and find containing sphere
  World::AtomComposite atoms = World::getInstance().getSelectedAtoms();
  if (atoms.size() == 0) {
    ELOG(2, "No atoms selected, aborting,");
    return Action::failure;
  }

  Vector center;
  double radius = 0.;
  {
    center.Zero();
    for(World::AtomComposite::const_iterator iter = atoms.begin();
        iter != atoms.end(); ++iter)
      center += (*iter)->getPosition();
    center *= 1./(double)atoms.size();
    for(World::AtomComposite::const_iterator iter = atoms.begin();
        iter != atoms.end(); ++iter) {
      const Vector &position = (*iter)->getPosition();
      for (size_t i=0;i<NDIM;++i) {
        const double temp_distance = fabs(position[i] - center[i]);
        if (temp_distance > radius)
          radius = temp_distance;
      }
    }
  }
  // always move cluster to be centered at origin
  for(World::AtomComposite::iterator iter = atoms.begin();
      iter != atoms.end(); ++iter)
    (*iter)->setPosition((*iter)->getPosition() - center);
  if (radius == 0.) { // never have zero radius: set to one or less
    ASSERT( (params.counts[0] != 0.) && (params.counts[1] != 0.) && (params.counts[2] != 0.),
        "FillRegularGridAction::performCall() - counts must not be zero.");
    const Vector diag(1./(double)(params.counts[0]+1), 1./(double)(params.counts[1]+1), 1./(double)(params.counts[2]+1));
    radius = std::min(1.,World::getInstance().getDomain().translateIn(diag).Norm()/2.);
  }
  // add some small boundary
  radius += 1e+6*std::numeric_limits<double>::epsilon();
  radius = std::max(params.mindistance, radius);
  LOG(2, "INFO: The " << atoms.size() << " are contained in a sphere at "
      << center << " with radius " << radius << ".");
  Shape *s = new Shape(Sphere(zeroVec, radius));
  LOG(2, "INFO: Created sphere at " << s->getCenter() << " and radius " << s->getRadius() << ".");

  // create set of atomic ids
  std::vector<atomId_t> atomIds(atoms.size(), (size_t)-1);
  std::transform(atoms.begin(), atoms.end(), atomIds.begin(),
      boost::bind(&atom::getId, _1) );
  Cluster::atomIdSet atomset(atomIds.begin(), atomIds.end());
  LOG(2, "INFO: Created atomset of size " << atomset.size() << ".");

  // create cluster
  ClusterInterface::Cluster_impl cluster(new Cluster(atomset, *s));
  LOG(2, "INFO: Created cluster of size " << atomset.size() << ".");

  // check for selected molecules and create surfaces from them
  const std::vector<molecule *> molecules = World::getInstance().getSelectedMolecules();
  typedef std::list<FillPredicate *> predicate_list_t;
  typedef std::list<LinkedCell_deprecated *> LC_list_t;
  typedef std::list<Tesselation *> Tesselation_list_t;
  predicate_list_t surface_predicate_list;
  LC_list_t LinkedCell_deprecated_list;
  Tesselation_list_t Tesselation_list;
  if (params.SphereRadius != 0.) {
    if ( molecules.size() == 0) {
      ELOG(1, "You have given a sphere radius " << params.SphereRadius
          << " != 0, but have not select any molecules.");
    }
    for (std::vector<molecule *>::const_iterator iter = molecules.begin();
        iter != molecules.end(); ++iter) {
      // create adaptor for the current molecule
      PointCloudAdaptor< molecule > cloud(*iter, (*iter)->name);
      LinkedCell_deprecated_list.push_back(
          new LinkedCell_deprecated(cloud, 2.*params.SphereRadius)
      );

      // create tesselation
      Tesselation_list.push_back(
          new Tesselation
      );
      (*Tesselation_list.back())(cloud, params.SphereRadius);

      // and create predicate
      surface_predicate_list.push_back(
          new FillPredicate( IsInsideSurface_FillPredicate(
              *Tesselation_list.back(),
              *LinkedCell_deprecated_list.back()) )
      );
    }
  }

  // create predicate, mesh, and filler
  {
    FillPredicate *voidnode_predicate = new FillPredicate( IsVoidNode_FillPredicate(*s) );
    FillPredicate Andpredicate = (*voidnode_predicate);
    for (predicate_list_t::iterator iter = surface_predicate_list.begin();
        iter != surface_predicate_list.end(); ++iter)
      Andpredicate = (Andpredicate) && !(**iter);
    Mesh *mesh = new CubeMesh(params.counts, params.offset, World::getInstance().getDomain().getM());
    Filler *fillerFunction = new Filler(*mesh, Andpredicate);

    // fill
    CopyAtoms_withBonds copyMethod;
    (*fillerFunction)(copyMethod, cluster);

    // remove
    delete fillerFunction;
    delete mesh;
    delete voidnode_predicate;
    for (LC_list_t::iterator iter = LinkedCell_deprecated_list.begin();
            iter != LinkedCell_deprecated_list.end(); ++iter)
      delete (*iter);
    LinkedCell_deprecated_list.clear();
    for (Tesselation_list_t::iterator iter = Tesselation_list.begin();
            iter != Tesselation_list.end(); ++iter)
      delete (*iter);
    Tesselation_list.clear();
    delete s;
  }

  // generate list of newly created atoms
  // (we can in general remove more quickly from a list than a vector)
  AtomVector filleratoms = World::getInstance().getAllAtoms();
//  LOG(0, filleratoms.size() << " atoms are present.");
  std::list<atom *> filleratoms_list;
  std::copy( filleratoms.begin(),  filleratoms.end(), std::back_inserter( filleratoms_list ));
//  LOG(0, filleratoms_list.size() << " atoms have been copied.");
  for (AtomVector::const_iterator iter = presentatoms.begin();
      iter != presentatoms.end();
      ++iter) {
    filleratoms_list.remove(*iter);
  }
//  LOG(0, filleratoms_list.size() << " atoms left after removal.");
  filleratoms.clear();
  std::copy(filleratoms_list.begin(), filleratoms_list.end(), std::back_inserter( filleratoms ));

//  LOG(0, filleratoms.size() << " atoms have been inserted.");

  return Action::state_ptr(new FillRegularGridState(filleratoms,params));
}

Action::state_ptr FillRegularGridAction::performUndo(Action::state_ptr _state) {
//  FillRegularGridState *state = assert_cast<FillRegularGridState*>(_state.get());
//
//  BOOST_FOREACH(atom *_atom, state->filleratoms) {
//    World::getInstance().destroyAtom(Walker);
//  }
//
//  // as atoms and atoms from state are removed, we have to create a new one
//  std::vector<atom*> filleratoms;
//  return Action::state_ptr(new FillRegularGridState(filleratoms,state->params));
  return Action::failure;
}

Action::state_ptr FillRegularGridAction::performRedo(Action::state_ptr _state){
  //FillRegularGridState *state = assert_cast<FillRegularGridState*>(_state.get());

  return Action::failure;
  //return Action::state_ptr(_state);
}

bool FillRegularGridAction::canUndo() {
  return false;
}

bool FillRegularGridAction::shouldUndo() {
  return false;
}
/** =========== end of function ====================== */