/* * 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 . */ /* * 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; } /** Internal function to decide whether a suitable model is present or not. * * Here, the heuristic for deciding whether a new linked cell structure has to * be constructed or not is implemented. The current heuristic is as follows: * -# the best model should have at least half the desired length (such * that most we have to look two neighbor shells wide and not one). * -# the best model should have at most twice the desired length but * no less than 1 angstroem. * * \note Dealing out a pointer is here (hopefully) safe because the function is * internal and we - inside this class - know what we are doing. * * @param distance edge length of the requested linked cell structure * @return NULL - there is no fitting LinkedCell_Model, else - pointer to instance */ const LinkedCell_Model *LinkedCell_Controller::getBestModel(double distance) const { /// Bound distance to be within [lower_threshold, upper_threshold). /// Note that we need to stay away from upper boundary a bit, /// otherwise the distance will end up outside of the interval. if (distance < lower_threshold) distance = lower_threshold; if (distance > 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) { /// Look for best instance. const LinkedCell_Model * const LCModel_best = getBestModel(distance); /// Construct new instance if none found, if (LCModel_best == NULL) { LinkedCell_Model * const LCModel_new = new LinkedCell_Model(distance, domain); LCModel_new->insertPointCloud(set); insertNewModel(distance, 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() )+","+toString(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) {}