/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 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 .
*/
/*
* LinkedCell_Controller.cpp
*
* Created on: Nov 15, 2011
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
//#include "CodePatterns/MemDebug.hpp"
#include
#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Observer/Notification.hpp"
#include "CodePatterns/Range.hpp"
#include "LinkedCell_Controller.hpp"
#include "LinkedCell_Model.hpp"
#include "LinkedCell_View.hpp"
#include "LinkedCell_View_ModelWrapper.hpp"
#include "IPointCloud.hpp"
#include "WorldTime.hpp"
using namespace LinkedCell;
double LinkedCell_Controller::lower_threshold = 1.;
double LinkedCell_Controller::upper_threshold = 20.;
/** Constructor of class LinkedCell_Controller.
*
*/
LinkedCell_Controller::LinkedCell_Controller(const Box &_domain) :
Observer("LinkedCell_Controller"),
domain(_domain)
{
// sign on to specific notifications
domain.signOn(this, Box::MatrixChanged);
WorldTime::getInstance().signOn(this, WorldTime::TimeChanged);
/// Check that upper_threshold fits within half the box.
Vector diagonal(1.,1.,1.);
diagonal.Scale(upper_threshold);
Vector diagonal_transformed = domain.getMinv() * diagonal;
double max_factor = 1.;
for (size_t i=0; i 1./max_factor)
max_factor = 1./diagonal_transformed.at(i);
upper_threshold *= max_factor;
/// Check that lower_threshold is still lower, if not set to half times upper_threshold.
if (lower_threshold > upper_threshold)
lower_threshold = 0.5*upper_threshold;
}
/** Destructor of class LinkedCell_Controller.
*
* Here, we free all LinkedCell_Model instances again.
*
*/
LinkedCell_Controller::~LinkedCell_Controller()
{
// sign off
domain.signOff(this, Box::MatrixChanged);
WorldTime::getInstance().signOff(this, WorldTime::TimeChanged);
/// we free all LinkedCell_Model instances again.
for(MapEdgelengthModel::iterator iter = ModelsMap.begin();
!ModelsMap.empty(); iter = ModelsMap.begin()) {
delete iter->second;
ModelsMap.erase(iter);
}
}
/** Internal function to obtain the range within which an model is suitable.
*
* \note We use statics lower_threshold and upper_threshold as min and max
* boundaries.
*
* @param distance desired egde length
* @return range within which model edge length is acceptable
*/
const range LinkedCell_Controller::getHeuristicRange(const double distance) const
{
const double lower = 0.5*distance < lower_threshold ? lower_threshold : 0.5*distance;
const double upper = 2.*distance > upper_threshold ? upper_threshold : 2.*distance;
range HeuristicInterval(lower, upper);
return HeuristicInterval;
}
static size_t checkLengthsLess(
const Vector &Lengths,
const double min_distance)
{
size_t counter=0;
for (size_t i=0;i0; --dim) {
const double fraction = pow(MAX_LINKEDCELLNODES+1., 1./(double)dim);
for (size_t i=0;i upper_threshold)
distance = upper_threshold - std::numeric_limits::round_error();
/// Look for all models within [0.5 distance, 2. distance).
MapEdgelengthModel::const_iterator bestmatch = ModelsMap.end();
if (!ModelsMap.empty()) {
for(MapEdgelengthModel::const_iterator iter = ModelsMap.begin();
iter != ModelsMap.end(); ++iter) {
// check that we are truely within range
range HeuristicInterval(getHeuristicRange(iter->first));
if (HeuristicInterval.isInRange(distance)) {
// if it's the first match or a closer one, pick
if ((bestmatch == ModelsMap.end())
|| (fabs(bestmatch->first - distance) > fabs(iter->first - distance)))
bestmatch = iter;
}
}
}
/// Return best match or NULL if none found.
if (bestmatch != ModelsMap.end())
return bestmatch->second;
else
return NULL;
}
/** Internal function to insert a new model and check for valid insertion.
*
* @param distance edge length of new model
* @param instance pointer to model
*/
void LinkedCell_Controller::insertNewModel(const double edgelength, const LinkedCell_Model* instance)
{
std::pair< MapEdgelengthModel::iterator, bool> inserter =
ModelsMap.insert( std::make_pair(edgelength, instance) );
ASSERT(inserter.second,
"LinkedCell_Controller::getView() - LinkedCell_Model instance with distance "
+toString(edgelength)+" already present.");
}
/** Returns the a suitable LinkedCell_Model contained in a LinkedCell_View
* for the requested \a distance.
*
* \sa getBestModel()
*
* @param distance edge length of the requested linked cell structure
* @param set of initial points to insert when new model is created (not always), should be World's
* @return LinkedCell_View wrapping the best LinkedCell_Model
*/
LinkedCell_View LinkedCell_Controller::getView(const double distance, IPointCloud &set)
{
// distance should be a given minimum length dependent on domain at least
const double ConstraintDistance = getMinimumSensibleLength(distance);
/// Look for best instance.
const LinkedCell_Model * const LCModel_best = getBestModel(ConstraintDistance);
/// Construct new instance if none found,
if (LCModel_best == NULL) {
LinkedCell_Model * const LCModel_new = new LinkedCell_Model(ConstraintDistance, domain);
LCModel_new->insertPointCloud(set);
insertNewModel(ConstraintDistance, LCModel_new);
LinkedCell_View view(*LCModel_new);
return view;
} else {
/// else construct interface and return.
LinkedCell_View view(*LCModel_best);
return view;
}
}
/** Internal function to re-create all present and used models for the new Box.
*
* This is necessary in the following cases:
* -# the Box is changed
* -# we step on to a different time step, i.e. all atomic positions change
*
* The main problem are the views currently in use.
*
* We make use of LinkedCell:LinkedCell_View::RAIIMap as there all present are
* listed. We go through the list, create a map with old model ref as keys to
* just newly created ones, and finally go again through each view and exchange
* the model against the new ones via a simple map lookup.
*
*/
void LinkedCell_Controller::updateModels()
{
LOG(1, "INFO: Updating all models.");
typedef std::map ModelLookup;
ModelLookup models;
// set up map, for now with NULL pointers
for (LinkedCell_View::ModelInstanceMap::const_iterator iter = LinkedCell_View::RAIIMap.begin();
iter != LinkedCell_View::RAIIMap.end(); ++iter) {
#ifndef NDEBUG
std::pair< ModelLookup::iterator, bool > inserter =
#endif
models.insert( std::pair( (*iter)->LC->getModel(), NULL) );
LOG(2, "INFO: Added " << (*iter)->LC->getModel() << " to list of models to replace.");
ASSERT( inserter.second,
"LinkedCell_Controller::updateModelsForNewBoxMatrix() - failed to insert old model "
+toString( (*iter)->LC->getModel() )+",NULL into models, is already present");
}
// invert MapEdgelengthModel
LOG(2, "INFO: ModelsMap is " << ModelsMap << ".");
typedef std::map MapEdgelengthModel_inverted;
MapEdgelengthModel_inverted ModelsMap_inverted;
for (MapEdgelengthModel::const_iterator iter = ModelsMap.begin();
iter != ModelsMap.end(); ++iter) {
#ifndef NDEBUG
MapEdgelengthModel_inverted::const_iterator assertiter = ModelsMap_inverted.find(iter->second);
ASSERT( assertiter == ModelsMap_inverted.end(),
"LinkedCell_Controller::updateModelsForNewBoxMatrix() - ModelsMap is not invertible, value "
+toString(iter->second)+" is already present.");
#endif
ModelsMap_inverted.insert( std::make_pair(iter->second, iter->first) );
}
LOG(2, "INFO: Inverted ModelsMap is " << ModelsMap_inverted << ".");
// go through map and re-create models
for (ModelLookup::iterator iter = models.begin(); iter != models.end(); ++iter) {
// delete old model
const LinkedCell_Model * const oldref = iter->first;
#ifndef NDEBUG
MapEdgelengthModel_inverted::const_iterator assertiter = ModelsMap_inverted.find(oldref);
ASSERT( assertiter != ModelsMap_inverted.end(),
"LinkedCell_Controller::updateModelsForNewBoxMatrix() - ModelsMap_inverted does not contain old model "
+toString(oldref)+".");
#endif
const double distance = ModelsMap_inverted[oldref];
// create new one, afterwards erase old model (this is for unit test to get different memory addresses)
LinkedCell_Model * const newref = new LinkedCell_Model(distance, domain);
MapEdgelengthModel::iterator oldmodeliter = ModelsMap.find(distance);
delete oldmodeliter->second;
ModelsMap.erase(oldmodeliter);
LOG(2, "INFO: oldref is " << oldref << ", newref is " << newref << ".");
iter->second = newref;
// replace in ModelsMap
#ifndef NDEBUG
std::pair< MapEdgelengthModel::iterator, bool > inserter =
#endif
ModelsMap.insert( std::make_pair(distance, newref) );
ASSERT( inserter.second,
"LinkedCell_Controller::updateModelsForNewBoxMatrix() - failed to insert new model "
+toString(distance)+","+toString(newref)+" into ModelsMap, is already present");
}
// remove all remaining active models (also those that don't have an active View on them)
for (MapEdgelengthModel::iterator iter = ModelsMap.begin();
!ModelsMap.empty();
iter = ModelsMap.begin()) {
delete iter->second;
ModelsMap.erase(iter);
}
// delete inverted map for safety (values are gone)
ModelsMap_inverted.clear();
// go through views and exchange the models
for (LinkedCell_View::ModelInstanceMap::const_iterator iter = LinkedCell_View::RAIIMap.begin();
iter != LinkedCell_View::RAIIMap.end(); ++iter) {
ModelLookup::const_iterator modeliter = models.find((*iter)->LC->getModel());
ASSERT( modeliter != models.end(),
"LinkedCell_Controller::updateModelsForNewBoxMatrix() - we miss a model "
+toString((*iter)->LC->getModel())+" in ModelLookup.");
// this is ugly but the only place where we have to set ourselves over the constness of the member variable
if (modeliter != models.end()) {
LOG(2, "INFO: Setting model to " << modeliter->second << " in view " << *iter << ".");
(*iter)->LC->setModel(modeliter->second);
}
}
}
/** Callback function for Observer mechanism.
*
* @param publisher reference to the Observable that calls
*/
void LinkedCell_Controller::update(Observable *publisher)
{
ELOG(2, "LinkedCell_Model received inconclusive general update from "
<< publisher << ".");
}
/** Callback function for the Notifications mechanism.
*
* @param publisher reference to the Observable that calls
* @param notification specific notification as cause of the call
*/
void LinkedCell_Controller::recieveNotification(Observable *publisher, Notification_ptr notification)
{
if (publisher == &domain) {
switch(notification->getChannelNo()) {
case Box::MatrixChanged:
LOG(1, "INFO: LinkedCell_Controller got update from Box.");
updateModels();
break;
default:
ASSERT(0,
"LinkedCell_Controller::recieveNotification() - unwanted notification from Box "
+toString(notification->getChannelNo())+" received.");
break;
}
} else if (publisher == WorldTime::getPointer()) {
switch(notification->getChannelNo()) {
case WorldTime::TimeChanged:
LOG(1, "INFO: LinkedCell_Controller got update from WorldTime.");
updateModels();
break;
default:
ASSERT(0,
"LinkedCell_Controller::recieveNotification() - unwanted notification from WorldTime "
+toString(notification->getChannelNo())+" received.");
break;
}
} else {
ELOG(1, "Notification " << notification->getChannelNo()
<< " from unknown publisher " << publisher << ".");
}
}
/** Callback function when an Observer dies.
*
* @param publisher reference to the Observable that calls
*/
void LinkedCell_Controller::subjectKilled(Observable *publisher)
{}