diff --git a/src/coupling/quasi_continuum.cc b/src/coupling/quasi_continuum.cc index 7e760fd..b1fd10a 100644 --- a/src/coupling/quasi_continuum.cc +++ b/src/coupling/quasi_continuum.cc @@ -1,247 +1,246 @@ /** * @file quasi_continuum.cc * * @author Till Junge * * @date Mon Jul 07 17:09:15 2014 * * @brief Quasicontinuum point to point coupling * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ #include "quasi_continuum.hh" #include "communicator.hh" #include "lib_bridging.hh" #include "lm_common.hh" #include "lm_log.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ QuasiContinuum::QuasiContinuum(const std::string &name) : LMObject(name), CouplingAtomicContinuum(name), interface_association("InterfaceAssociation:" + name), interface_path(""), pad_association("PADAssociation:" + name), pad_path(""), weight_path("") { } /* -------------------------------------------------------------------------- */ QuasiContinuum::~QuasiContinuum() {} /* -------------------------------------------------------------------------- */ template void QuasiContinuum::init(DomainA &domA, DomainC &domC) { const UInt Dim = spatial_dimension; std::vector pbc_flags(spatial_dimension); for (UInt i = 0; i < Dim; ++i) { if (this->is_in_atomic) pbc_flags[i] = domA.isPeriodic(i); if (this->is_in_continuum) pbc_flags[i] = domC.isPeriodic(i); } this->interface_association.setParam("PBC", pbc_flags); this->pad_association.setParam("PBC", pbc_flags); if (!this->arlequin_flag && this->reset_position_freq != 1) { LM_FATAL("Without arlequin, reset_position_freq must be 1 (which is the " "default)"); } STARTTIMER("Interface Association"); this->interface_association.setParam("PATH", this->interface_path); this->interface_association.init(domA, domC); this->setBoundary(); STOPTIMER("Interface Association"); STARTTIMER("Pad Association"); this->pad_association.setParam("WEIGHT_PATH", this->weight_path); this->pad_association.setParam("PATH", this->pad_path); this->pad_association.init(domA, domC); this->clearPadForces(); this->clearInterfaceForces(); STOPTIMER("Pad Association"); } /* -------------------------------------------------------------------------- */ void QuasiContinuum::setBoundary() { if (this->is_in_continuum) { LM_TOIMPLEMENT; // typedef FilterCompatibility PointsBFilter; // auto *ptr = // static_cast(FilterManager::getManager().getObject( // "pointsBlist:InterfaceAssociation:" + this->name)); // UInt index = 0; // for (auto &&point : *ptr) { // for (UInt dir = 0; dir < Dim; ++dir) { // point.boundary(dir) = true; // } // ++index; // } // DUMP("SetBoundary fully blocked " << index << " nodes.", DBG_DETAIL); } } /* -------------------------------------------------------------------------- */ void QuasiContinuum::updateForMigration() { LM_TOIMPLEMENT; // Communicator &comm = Communicator::getCommunicator(); // comm.synchronize(Communicator::group_all); STARTTIMER("updatesForMigrations"); this->interface_association.updateForMigration(); STOPTIMER("updatesForMigrations"); LM_TOIMPLEMENT; // comm.synchronize(group_all); STARTTIMER("updatesForMigrations"); this->pad_association.updateForMigration(); STOPTIMER("updatesForMigrations"); LM_TOIMPLEMENT; // comm.synchronize(group_all); } /* -------------------------------------------------------------------------- */ template void QuasiContinuum::coupling(CouplingStage stage, DomainA &domA, DomainC &domC) { if (stage == COUPLING_STEP1) { updateForMigration(); if (current_step % reset_position_freq == 0) { this->enforceInterface<_displacement>(); this->enforcePad<_position>(); } } else if (stage == COUPLING_STEP2) { this->clearPadForces(); } // } else if (stage == COUPLING_STEP2) { // STARTTIMER("clearPadForces"); // this->clearPadForces(); // //this->pad_association.template mixField<_force>(); // STOPTIMER("clearPadForces"); // STARTTIMER("clearInterfaceForces"); // this->clearInterfaceForces(); // //this->interface_association.template mixField<_force>(); // STOPTIMER("clearInterfaceForces"); // } else if (arlequin_flag && stage == COUPLING_STEP3) { this->pad_association.template mixField<_velocity>(domA.getContainer(), domC.getContainer()); this->enforceInterface<_velocity>(); } } /* -------------------------------------------------------------------------- */ template void QuasiContinuum::enforceInterface() { this->interface_association.template enforceField(); } /* -------------------------------------------------------------------------- */ template void QuasiContinuum::enforcePad() { this->pad_association.template enforceField(); } /* -------------------------------------------------------------------------- */ void QuasiContinuum::clearPadForces() { this->pad_association.zeroFields<_force>(); this->pad_association.zeroFields<_velocity>(); } /* -------------------------------------------------------------------------- */ void QuasiContinuum::clearInterfaceForces() { this->interface_association.zeroFields<_force>(); this->interface_association.zeroFields<_velocity>(); } /* -------------------------------------------------------------------------- */ void QuasiContinuum::clearAll() { this->interface_association.clearAll(); this->pad_association.clearAll(); } /* -------------------------------------------------------------------------- */ /* LMDESC QUASICONTINUUM This coupler implements a CADD/QC style sharp interface coupling between atomic and continuum domains. The interface and pad atom-nodes are synchronised using PointAssociation couplers. Based on these, the coupler creates dumpable containers for: - 1) the interface nodes - (pointsBlist:InterfaceAssociation:) - 2) the interface atoms - (pointsAlist:InterfaceAssociation:) - 3) the pad nodes (pointsAlist:PadAssociation:) - 4) the pad atoms (pointsBlist:PadAssociation:) + + 1) the interface nodes (pointsBlist:InterfaceAssociation:) + 2) the interface atoms (pointsAlist:InterfaceAssociation:) + 3) the pad nodes (pointsAlist:PadAssociation:) + 4) the pad atoms (pointsBlist:PadAssociation:) */ /* LMEXAMPLE COUPLING_CODE ID md fe QUASICONTINUUM PAD_PATH ./pad_atoms.txt INTERFACE_PATH ./interface_atoms.txt GEOMETRY pad GRID_DIVISION 8 8 4 */ /* LMHERITANCE point_association */ void QuasiContinuum::declareParams() { // so they can parse their respective tolerances this->addSubParsableObject(pad_association); this->addSubParsableObject(interface_association); /* LMKEYWORD PAD_PATH String file path of the file containing pad atoms */ this->parseKeyword("PAD_PATH", pad_path); /* LMKEYWORD INTERFACE_PATH String file path of the file containing pad atoms */ this->parseKeyword("INTERFACE_PATH", interface_path); /* LMKEYWORD WEIGHT_SCRIPT String file path of the script to use to generate the weights */ this->parseKeyword("WEIGHT_SCRIPT", weight_path, ""); /* LMKEYWORD POS_FREQ The frequency at which the positions should be reset */ this->parseKeyword("POS_FREQ", reset_position_freq, 1u); /* LMKEYWORD ARLEQUIN Set the arlequin style flavor */ this->parseTag("ARLEQUIN", arlequin_flag, false); } /* -------------------------------------------------------------------------- */ DECLARE_COUPLER_INIT_MAKE_CALL(QuasiContinuum, domMD, domCONTINUUM) DECLARE_COUPLING_MAKE_CALL(QuasiContinuum, domMD, domCONTINUUM) /* ------------------------------------------------------------------------ */ __END_LIBMULTISCALE__ diff --git a/src/dumper/dumper_grid.cc b/src/dumper/dumper_grid.cc index 29c05ff..3420cf3 100644 --- a/src/dumper/dumper_grid.cc +++ b/src/dumper/dumper_grid.cc @@ -1,231 +1,230 @@ /** * @file dumper_grid.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This dumper allows to use a grid so as to reduce on a per-cell basis * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dumper_grid.hh" #include "cube.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "lm_file.hh" #include "ref_point_data.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template void DumperGrid::dump(Cont &cont) { switch (field) { case _displacement: FillGrid<_displacement>(cont); break; case _velocity: FillGrid<_velocity>(cont); break; case _force: FillGrid<_force>(cont); break; default: LM_FATAL("cannot treat this field"); } std::stringstream fname; fname << this->getBaseName() << "_grid." << std::setfill('0') << std::setw(4) << this->action_step << ".txt"; LMFile file(fname.str()); // std::vector *> &blocks = grid->getBlocks(); LM_TOIMPLEMENT; // Real x = 0, y = 0, z = 0; // std::vector force_x; // for (UInt i = 0; i < grid->getSize(); ++i) // if (blocks[i]) { // grid->index2Coordinates(i, x, y, z); // Tensor t; // t = grid->extractBoxValue(i, op); // file.printf( // "%d\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\n", // i, x, y, z, t[0], t[1], t[2], t[3], t[4], t[5]); // } // file.close(); // delete (grid); // grid = NULL; } /* -------------------------------------------------------------------------- */ template void DumperGrid::FillGrid(Cont &cont) { if (cont.size() == 0) return; LM_TOIMPLEMENT; // Cube &cube = cont.getBoundingBox(); // constexpr UInt Dim = Cont::Dim; // for (UInt i = 0; i < Dim; ++i) { // if (gridspace[i] == -1 && griddivision[i] == UINT_MAX) // LM_FATAL("GRIDSPACE or GRIDDIVISION keyword should be used for // direction " // << i); // } // Vector gridspace2 = gridspace.block(0,0); // Vector griddivision2 = griddivision.block(0,0); // grid = // new SpatialGridPoint(cube, gridspace2, // griddivision2); // for (auto &&at : cont) { // Vector x = at.position0(); // if (current_flag) // x = at.position(); // auto tmp = at.template field(); // Tensor t; // for (UInt i = 0; i < Dim; ++i) // t[i] = tmp[i]; // grid->addElement(t, x); // } // DUMP("we have as average " << grid->getAverageEltByBlock() // << " points per block", // DBG_MESSAGE); } /* -------------------------------------------------------------------------- */ /* LMDESC GRID This dumper allows to use a grid so as to reduce on a per-cell basis. A point/nodes based sorting into a spatial grid is used. It admits several options that could be model dependent. */ /* LMEXAMPLE DUMPER GRID GRIDSPACE 10 10 10 PREFIX /home/titeuf */ /* LMHERITANCE dumper_interface */ inline void DumperGrid::declareParams() { DumperInterface::declareParams(); /* LMKEYWORD GRID_DIVISIONX Set the number of grid cells in x direction */ this->parseKeyword("GRID_DIVISIONX", griddivision[0]); /* LMKEYWORD GRID_DIVISIONY Set the number of grid cells in y direction */ this->parseKeyword("GRID_DIVISIONY", griddivision[1]); /* LMKEYWORD GRID_DIVISIONZ Set the number of grid cells in x direction */ this->parseKeyword("GRID_DIVISIONZ", griddivision[2]); /* LMKEYWORD GRID_DIVISION Allows to set the number of grid cells in all directions of space, see GRID\_DIVISION\_X/Y/Z */ this->parseVectorKeyword("GRID_DIVISION", 3, griddivision); /* LMKEYWORD GRIDSPACEX Set the cell size in the X direction for the grid generation */ this->parseKeyword("GRIDSPACEX", gridspace[0]); /* LMKEYWORD GRIDSPACEY Set the cell size in the Y direction for the grid generation */ this->parseKeyword("GRIDSPACEY", gridspace[1]); /* LMKEYWORD GRIDSPACEZ Set the cell size in the Z direction for the grid generation */ this->parseKeyword("GRIDSPACEZ", gridspace[2]); /* LMKEYWORD GRIDSPACE The implementation of the Bridging domain method needs a match between elements and atoms to be made at initialisation of the bridging zones. This could be a very costful operation. \\ In order to break the complexity of such a stage, the elements are first stored in a grid for which number of cells (in each direction) are specified by this keyword. */ this->parseVectorKeyword("GRIDSPACE", 3, gridspace); /* LMKEYWORD OPERATOR Specify the operator to reduce the information per cell \\ The possible operators are for now: - \begin{itemize} - \item SUM: perform a sum per cell $\sum_{i \in cell} v_i$ - \item AVERAGE: perform an average per cell $\frac{1}{\#cell}\sum_{i \in - cell} v_i$ - \item DENSITY: perform the volumetric average $\overline{v} = - \frac{1}{vol(cell)}\sum_{i \in cell} v_i$ - \item DEVIATION: compute the standard deviation $\frac{1}{\#cell}\sum_{i - \in cell} (v_i - \overline{v} )^2$ - \item MAX: compute the max value - \item MIN: compute the min value - \end{itemize} + + - SUM: perform a sum per cell :math:`\sum_{i \in cell} v_i` + - AVERAGE: perform an average per cell :math:`\frac{1}{\#cell}\sum_{i \in + cell} v_i` + - DENSITY: perform the volumetric average :math:`\overline{v} = + \frac{1}{vol(cell)}\sum_{i \in cell} v_i` + - DEVIATION: compute the standard deviation :math:`\frac{1}{\#cell}\sum_{i + \in cell} (v_i - \overline{v} )^2` + - MAX: compute the max value + - MIN: compute the min value */ this->parseKeyword("OPERATOR", op); /* LMKEYWORD CURRENT_POSITION All the points/atoms/nodes are sorted in a spatial grid. By default the initial position of these points/atoms/nodes are considered. This flag allows to consider current positions instead. */ this->parseTag("CURRENT_POSITION", current_flag); /* LMKEYWORD FIELD Request to perform the grid reduction on the provided field */ this->parseKeyword("FIELD", field); } /* -------------------------------------------------------------------------- */ DumperGrid::DumperGrid(const std::string &name) : LMObject(name) { current_flag = 0; for (UInt i = 0; i < 3; ++i) gridspace[i] = -1.; for (UInt i = 0; i < 3; ++i) griddivision[i] = UINT_MAX; op = OP_AVERAGE; // grid = NULL; } /* -------------------------------------------------------------------------- */ DumperGrid::~DumperGrid() {} /* -------------------------------------------------------------------------- */ DECLARE_DUMPER_MAKE_CALL(DumperGrid) __END_LIBMULTISCALE__ diff --git a/src/factory/domain_multiscale.cc b/src/factory/domain_multiscale.cc index 4341f0d..a111c65 100644 --- a/src/factory/domain_multiscale.cc +++ b/src/factory/domain_multiscale.cc @@ -1,454 +1,451 @@ /** * @file domain_multiscale.cc * * @author Guillaume Anciaux * * @date Wed Jul 02 16:25:38 2014 * * @brief This is a meta domain whose purpose is to handle other domains * * @section LICENSE * * Copyright INRIA and CEA * * The LibMultiScale is a C++ parallel framework for the multiscale * coupling methods dedicated to material simulations. This framework * provides an API which makes it possible to program coupled simulations * and integration of already existing codes. * * This Project was initiated in a collaboration between INRIA Futurs Bordeaux * within ScAlApplix team and CEA/DPTA Ile de France. * The project is now continued at the Ecole Polytechnique Fédérale de Lausanne * within the LSMS/ENAC laboratory. * * This software is governed by the CeCILL-C license under French law and * abiding by the rules of distribution of free software. You can use, * modify and/ or redistribute the software under the terms of the CeCILL-C * license as circulated by CEA, CNRS and INRIA at the following URL * "http://www.cecill.info". * * As a counterpart to the access to the source code and rights to copy, * modify and redistribute granted by the license, users are provided only * with a limited warranty and the software's author, the holder of the * economic rights, and the successive licensors have only limited * liability. * * In this respect, the user's attention is drawn to the risks associated * with loading, using, modifying and/or developing or reproducing the * software by the user in light of its specific status of free software, * that may mean that it is complicated to manipulate, and that also * therefore means that it is reserved for developers and experienced * professionals having in-depth computer knowledge. Users are therefore * encouraged to load and test the software's suitability as regards their * requirements in conditions enabling the security of their systems and/or * data to be ensured and, more generally, to use and operate it in the * same conditions as regards security. * * The fact that you are presently reading this means that you have had * knowledge of the CeCILL-C license and that you accept its terms. * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "action_interface.hh" #include "coupler_manager.hh" #include "domain_multiscale.hh" #include "factory_multiscale.hh" #include "filter_interface.hh" #include "geometry_manager.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "lm_communicator.hh" #include "lm_parser.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ DomainMultiScale::~DomainMultiScale() {} /* -------------------------------------------------------------------------- */ DomainMultiScale::DomainMultiScale() : LMObject("DomainMultiscale") {} /* -------------------------------------------------------------------------- */ void DomainMultiScale::coupling(CouplingStage stage) { #ifdef LIBMULTISCALE_USE_EPSN if (epsn_itfc->isNodeInitialized()) epsn_itfc->beginTask("updatecoupling"); #endif CouplerManager::getManager().coupling(stage); #ifdef LIBMULTISCALE_USE_EPSN if (epsn_itfc->isNodeInitialized()) epsn_itfc->endTask("updatecoupling"); #endif } /* -------------------------------------------------------------------------- */ UInt DomainMultiScale::build(const std::string &config) { Parser::parseConfigFile(config, "MultiScale", "", (*this)); // if (action_manager) action_manager->init(); // if (filter) filter->init(); Communicator::getGroup("all").synchronize(); #ifdef LIBMULTISCALE_USE_EPSN UInt epsn_status = 0; // EPSN Ready DUMP("EPSN is going to be READY..."); epsn_status = epsn_itfc->ready(); if (!epsn_status) { FATAL("error on EPSN Node ready operation\n"); } DUMP("EPSN READY"); // the main loop DUMP("EPSN Begin HTG"); epsn_itfc->beginHTM(); epsn_itfc->beginLoop("main"); DUMP("simulation ready to start for EPSN"); #endif return Dim; } /* -------------------------------------------------------------------------- */ #define BOOST_PARSE_MODEL(keyword, class, dim) \ DUMPFILE(Parser::fout, "trying to create class " \ << stringify_macro(class) \ << " while dimension is " << Dim << " dim is " \ << dim << " key is " << key << " keyword is " \ << keyword); \ if (key == keyword && dim == Dim) { \ if (Communicator::getCommunicator().getNbGroups() == 0) \ LM_FATAL("PROCESSOR KEYWORD SHOULD BE USED ON ID " << ID); \ auto &dom_group = Communicator::getCommunicator().getGroup(ID); \ model = std::make_shared(ID, dom_group); \ Parser::parseConfigFile(config_file, keyword, ID, *model); \ DUMP("Loading of config file ... done", DBG_INFO_STARTUP); \ if (model == NULL) \ LM_FATAL("Could not instanciate the model: " \ << " conf line was " << line.str()); \ model->checkAllKeywordsAreParsed(); \ model->init(); \ addObject(model); \ DUMPFILE(Parser::fout, \ "creation of " << stringify_macro(class) << " successful"); \ \ return true; \ } #define BOOST_PARSE_MODEL_KEY(n, data, class) \ BOOST_PARSE_MODEL(stringify_macro(BOOST_PP_TUPLE_ELEM(3, 2, class)), \ BOOST_PP_TUPLE_ELEM(3, 0, class), \ BOOST_PP_TUPLE_ELEM(3, 1, class)) /* -------------------------------------------------------------------------- */ UInt DomainMultiScale::parseLineModel(std::stringstream &line) { // first of all check if the communicator is already set correctly try { Communicator &myComm __attribute__((unused)) = Communicator::getCommunicator(); } catch (LibMultiScaleException &e) { if (lm_world_size == 1) Communicator::createCommunicator(); else LM_THROW( "Communicator was not set: " << " before MODEL keyword the COM keyword should have been used"); } std::string key, config_file; std::shared_ptr model; LMID ID = invalidDomain; // parse the keyword of the model Parser::parse(key, line); // parse the ID of the model Parser::parse(ID, line); // parse the config file of the model try { Parser::parse(config_file, line); } catch (LibMultiScaleException &e) { config_file = Parser::getCurrentConfigFile(); } BOOST_PP_SEQ_FOR_EACH(BOOST_PARSE_MODEL_KEY, f, LIST_ATOM_MODEL) BOOST_PP_SEQ_FOR_EACH(BOOST_PARSE_MODEL_KEY, f, LIST_CONTINUUM_MODEL) BOOST_PP_SEQ_FOR_EACH(BOOST_PARSE_MODEL_KEY, f, LIST_DD_MODEL) LM_FATAL("Unknown Model type " << key << " maybe you should compile support " "for this plugin ? or check the " "dimension required ?"); return false; } /* -------------------------------------------------------------------------- */ /* LMDESC MultiScale Section (toplevel domain) This section describe the most general keyword that can be used. It allows to specify general values of the simulation to be run. */ /* LMSYNTAX TODO */ /* LMEXAMPLE TODO */ ParseResult DomainMultiScale::setParam(const std::string &key, std::stringstream &line) { if (key == "MD_CODE" or key == "ELAST_CODE" or key == "DD_CODE") return parseLineModel(line); /* LMKEYWORD MODEL String Integer config_file(current_file) Declare a new domain to be created based on the plugins compiled. Up to now the possible models are: - \begin{itemize} - \item LAMMPS - \item MD1D - \item MECA1D - \item SIMULPACK - \item LIBMESH - \end{itemize} - - More to come in future (Akantu). \\ + + - LAMMPS + - MD1D + - MECA1D + - AKANTU + - PARADIS */ else if (key == "MODEL") return parseLineModel(line); /* LMKEYWORD COUPLING_CODE ARG TODO */ else if (key == "COUPLING_CODE") { return CouplerManager::getManager().parseCouplingLine(line, Dim); } /* LMKEYWORD GEOMETRY ARG TODO */ else if (key == "GEOMETRY") { if (!GeometryManager::allocated()) LM_FATAL("Please use DIMENSION keyword before" << Parser::getParserState()); return GeometryManager::getManager().parseLine(line, Dim); } /* LMKEYWORD DIMENSION Integer TODO */ else if (key == "DIMENSION") { Parser::parse((UInt &)Dim, line); spatial_dimension = Dim; DUMP("parsing DIMENSION commande = " << Dim, DBG_INFO_STARTUP); ActionManager::getManager(); FilterManager::getManager(); CouplerManager::getManager(); GeometryManager::getManager(); return true; } /* LMKEYWORD DUMPER String TODO \ref{dumpersSection} */ else if (key == "DUMPER") { if (!ActionManager::allocated()) LM_FATAL("Please use DIMENSION keyword before" << Parser::getParserState()); return ActionManager::getManager().parseActionLine(line, at_dumper); } /* LMKEYWORD FILTER String TODO \ref{filterSection} */ else if (key == "FILTER") { if (!FilterManager::allocated()) LM_FATAL("Please use DIMENSION keyword before" << Parser::getParserState()); return FilterManager::getManager().parseActionLine(line, at_filter); } else if (key == "COMPUTE") { if (!FilterManager::allocated()) LM_FATAL("Please use DIMENSION keyword before" << Parser::getParserState()); return FilterManager::getManager().parseActionLine(line, at_compute); } /* LMKEYWORD STIMULATION String TODO \ref{stimulatorsSection} */ else if (key == "STIMULATION") { if (!ActionManager::allocated()) LM_FATAL("Please use DIMENSION keyword before" << Parser::getParserState()); return ActionManager::getManager().parseActionLine(line, at_stimulator); } /* LMKEYWORD UNITCODE String TODO */ else if (key == "UNITCODE") { DUMP("parsing UNITCODE commande = " << line.str(), DBG_INFO_STARTUP); Parser::parse(code_unit_system, line); Parser::units_converter.setOutUnits(code_unit_system); Parser::units_converter.computeConversions(); if (lm_my_proc_id == 0) code_unit_system.print(); return true; } /* LMKEYWORD PROCESSORS Integer Integer TODO */ else if (key == "PROCESSORS") { LMID id; Parser::parse(id, line); UInt nb_procs = 0; Parser::parse(nb_procs, line); if (nb_procs == 0) LM_FATAL("a group cannot be constituted of zero processors"); DUMPFILE(Parser::fout, "parsing " << key << "( Integer , Integer ) with value " << id << " " << nb_procs); try { Communicator::getCommunicator().addGroup(id, nb_procs); } catch (LibMultiScaleException &e) { LM_FATAL_RE(e, "COM keyword should be used before the PROCESSOR keyword"); } return true; } /* LMKEYWORD COM String TODO */ else if (key == "COM") { std::string mot; Parser::parse(mot, line); Communicator::createCommunicator(); return true; } // /* LMKEYWORD GLOBAL_TIMESTEP Time // TODO // */ // else if (key == "GLOBAL_TIMESTEP") { // Real dt; // Parser::parse(dt, line); // this->setTimeStep(dt); // DUMPBYPROC("changing global timestep to " << dt, DBG_MESSAGE, 0); // return true; // } // /* LMKEYWORD SET_CURRENT_STEP Integer // TODO // */ // else if (key == "SET_CURRENT_STEP") { // UInt ts; // Parser::parse(ts, line); // this->setCurrentStep(ts); // DUMPBYPROC("changing current timestep to " << ts, DBG_MESSAGE, 0); // return true; // } /* LMKEYWORD RESTART_FILE String TODO */ else if (key == "RESTART_FILE") { LMID domainID; Parser::parse(domainID, line); std::string mot; Parser::parse(mot, line); DUMP("reloading file " << mot << " for domain " << domainID, DBG_MESSAGE); if (getSharedObject(domainID)) { if (Communicator::getGroup(domainID).amIinGroup()) getSharedObject(domainID)->readXMLFile(mot); } else LM_FATAL("there is no domain with ID " << domainID); return true; } Communicator::getGroup("all").synchronize(); return false; } /* -------------------------------------------------------------------------- */ std::shared_ptr DomainMultiScale::getModelByNumber(UInt i) { auto iest_domain = std::next(insert_order.begin(), i); return shared_pointer_objects[*iest_domain]; } /* -------------------------------------------------------------------------- */ UInt DomainMultiScale::getNbModels() { return objects.size(); } /* -------------------------------------------------------------------------- */ #define DOMAINLOOP(inloop) \ \ for (auto &obj : objects) { \ if (Communicator::getGroup(obj.first).amIinGroup()) \ inloop; \ } #define DEFINE_RETURN_REAL_FUNCTION_BYSUM(name) \ Real DomainMultiScale::name() { \ Real res = 0; \ DOMAINLOOP(res += obj.second->name()) \ return res; \ } /* -------------------------------------------------------------------------- */ Component &getRegisteredComponent(const std::string &component_name) { auto managers = std::forward_as_tuple( DomainMultiScale::getManager(), FilterManager::getManager(), CouplerManager::getManager(), ActionManager::getManager()); Component *ptr = nullptr; int found = 0; auto find_comp = [&](auto &manager) { try { auto &component = manager.getObject(component_name); ptr = &component; ++found; } catch (UnknownID &e) { } }; std::apply([&](auto &&... manager) { (find_comp(manager), ...); }, managers); if (found > 1) { LM_FATAL("several components are registered with the name: " << component_name << " " << ptr); } if (not ptr) { LM_THROW("component " << component_name << " cannot be found" << std::endl); } return *ptr; } /* -------------------------------------------------------------------------- */ OutputContainer getRegisteredOutput(const std::string &name) { std::string component_name = name; std::string output_name = ""; auto pos = component_name.find("."); if (pos != std::string::npos) { output_name = component_name.substr(pos + 1); component_name = component_name.substr(0, pos); } auto &comp = getRegisteredComponent(component_name); if (output_name != "") return comp.evalOutput(output_name); OutputContainer out; out = comp; return out; } /* -------------------------------------------------------------------------- */ template <> std::unique_ptr> IDManager::static_pointer = NULL; /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ diff --git a/src/filter/compute_reduce.cc b/src/filter/compute_reduce.cc index dd83b8a..994edc0 100644 --- a/src/filter/compute_reduce.cc +++ b/src/filter/compute_reduce.cc @@ -1,127 +1,124 @@ /** * @file compute_reduce.cc * * @author Guillaume Anciaux * * @date Tue Dec 10 16:00:49 2013 * * @brief This is a compute which reduce to scalar quantities * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "compute_reduce.hh" #include "lib_md.hh" #include "lm_common.hh" #include "math_tools.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ ComputeReduce::ComputeReduce(const std::string &name) : LMObject(name) { this->createInput("input"); this->createArrayOutput("reduced"); } /* -------------------------------------------------------------------------- */ ComputeReduce::~ComputeReduce() {} /* -------------------------------------------------------------------------- */ template std::enable_if_t> reduce(Cont &cont, V &reduced_data) { MathTools::ReduceOperator ops; ops.init(reduced_data.rows(), reduced_data.cols()); for (auto &&val : cont.rowwise()) { ops.account(val); } reduced_data = ops.result(cont.getCommGroup()); } /* -------------------------------------------------------------------------- */ template void ComputeReduce::build(Cont &cont) { const UInt Dim = cont.getDim(); Array reduced_data(1, Dim); auto &array = this->getOutputAsArray(); array.resize(array.rows(), Dim); switch (op) { case OP_MAX: reduce(cont, reduced_data); break; case OP_MIN: reduce(cont, reduced_data); break; case OP_SUM: reduce(cont, reduced_data); break; case OP_AVERAGE: reduce(cont, reduced_data); break; default: LM_FATAL("Operator " << op << " not handled"); } this->clear(); this->getOutputAsArray().push_back(reduced_data); } /* -------------------------------------------------------------------------- */ /* LMDESC REDUCE This computes a reduced version of a set of reals */ /* LMEXAMPLE COMPUTE disp EXTRACT INPUT md FIELD displacement\\ COMPUTE sum REDUCE INPUT disp OPERATOR MAX */ inline void ComputeReduce::declareParams() { ComputeInterface::declareParams(); /* LMKEYWORD OPERATOR Specify the operator to use for the reduction. \\ At present time the available operations are: - \begin{itemize} - \item MIN - \item MAX - \item SUM - \item AVERAGE - \end{itemize} - + - MIN + - MAX + - SUM + - AVERAGE */ this->parseKeyword("OPERATOR", op); } /* -------------------------------------------------------------------------- */ DECLARE_COMPUTE_MAKE_CALL(ComputeReduce) __END_LIBMULTISCALE__ diff --git a/src/filter/filter_interface.cc b/src/filter/filter_interface.cc index eb6ec89..6f5ba4e 100644 --- a/src/filter/filter_interface.cc +++ b/src/filter/filter_interface.cc @@ -1,65 +1,63 @@ /** * @file filter.cc * * @author Guillaume Anciaux * * @date Tue Jul 29 17:03:26 2014 * * @brief This is the mother class of all filters * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "filter_interface.hh" #include "lib_continuum.hh" #include "lm_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ /* LMDESC FILTER This is a generic internal filter */ void FilterInterface::declareParams() { ActionInterface::declareParams(); /* LMKEYWORD DOF_TYPE Allow to specify the dof type to iterate on: useful only for parallelism. The possible choices are: - \begin{itemize} - \item dt\_master : only the master dofs - \item dt\_slave : only the slave dofs - \item dt\_pbc_master : only pbc master - \item dt\_pbc_slave : only pbc slaves - \item dt\_all : all dofs - \end{itemize} + - dt_master : only the master dofs + - dt_slave : only the slave dofs + - dt_pbc_master : only pbc master + - dt_pbc_slave : only pbc slaves + - dt_all : all dofs */ this->parseKeyword("DOF_TYPE", dt, dt_master); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_langevin.cc b/src/stimulation/stimulation_langevin.cc index ee63f3f..67b05e6 100644 --- a/src/stimulation/stimulation_langevin.cc +++ b/src/stimulation/stimulation_langevin.cc @@ -1,390 +1,395 @@ /** * @file stimulation_langevin.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This stimulator applies a Langevin thermostat * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ #include "stimulation_langevin.hh" #include "container_array.hh" #include "factory_multiscale.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "math_tools.hh" #include "reference_manager.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ StimulationLangevin::StimulationLangevin(const std::string &name) : LMObject(name) { damp = 100; temperature = 100; timestep = 1.; seed = 0; compute_work_flag = false; // initialize the compute auto output_names = std::list{ "temperature", "friction_work", "restitution_work", "cumulated_frictional_work", "cumulated_restitution_work"}; this->createArrayOutputs(output_names); for (auto f : output_names) { this->getOutputAsArray(f).assign(1, 0); } // initialize the per atom computes output_names = std::list{ "old_friction_per_atom", "old_restitution_per_atom", "friction_work_per_atom", "restitution_work_per_atom" }; this->createArrayOutputs(output_names); cumulated_friction_work = 0; cumulated_restitution_work = 0; } /* -------------------------------------------------------------------------- */ StimulationLangevin::~StimulationLangevin() {} /* -------------------------------------------------------------------------- */ template void StimulationLangevin::stimulate(Cont &cont) { if (compute_work_flag && compute_work_peratom_flag) this->addForce(cont); else if ((!compute_work_flag) && compute_work_peratom_flag) this->addForce(cont); else if ((!compute_work_flag) && (!compute_work_peratom_flag)) this->addForce(cont); else if (compute_work_flag && (!compute_work_peratom_flag)) this->addForce(cont); else LM_FATAL("internal error: should not happend"); } /* -------------------------------------------------------------------------- */ template void StimulationLangevin::resetOldWorkPerAtom(Cont &cont) { constexpr UInt Dim = Cont::Dim; if (this->compute_work_peratom_flag) { getOutputAsArray("friction_work_per_atom").resize(cont.size(), Dim); getOutputAsArray("restitution_work_per_atom").resize(cont.size(), Dim); } auto &old_friction_per_atom = getOutputAsArray("old_friction_per_atom"); auto &old_restitution_per_atom = getOutputAsArray("old_restitution_per_atom"); DUMP("AAAAA " << cont.size() << " " << old_friction_per_atom.size() << " " << old_friction_per_atom.rows() << " " << old_friction_per_atom.cols(), DBG_DETAIL); if (old_friction_per_atom.size() == 0) { old_friction_per_atom.resize(cont.size(), Dim); old_restitution_per_atom.resize(cont.size(), Dim); cont.attachObject(old_friction_per_atom); cont.attachObject(old_restitution_per_atom); } else { if (old_friction_per_atom.size() != cont.size() * Dim) LM_FATAL("inconsistent previous storage of force vector (friction) " << old_friction_per_atom.size() << "(" << &old_friction_per_atom << ")" << " != " << cont.size() * Dim << "(" << &cont << ")"); if (old_restitution_per_atom.size() != cont.size() * Dim) LM_FATAL("inconsistent previous storage of force vector (restitution) " << old_restitution_per_atom.size() << "(" << &old_restitution_per_atom << ")" << " != " << cont.size() * Dim << "(" << &cont << ")"); } } /* -------------------------------------------------------------------------- */ template void StimulationLangevin::gatherWork(Cont &cont) { // Communicator &comm = Communicator::getCommunicator(); CommGroup &group = cont.getCommGroup(); if (compute_work_flag) { group.allReduce(&friction_work, 1, "frictional work", OP_SUM); group.allReduce(&restitution_work, 1, "restitution work", OP_SUM); friction_work *= 0.5 * timestep * code_unit_system.fd2e; restitution_work *= 0.5 * timestep * code_unit_system.fd2e; cumulated_friction_work += friction_work; cumulated_restitution_work += restitution_work; std::map toto{ {"temperature", temperature}, {"friction_work", friction_work}, {"restitution_work", restitution_work}, {"cumulated_frictional_work", cumulated_friction_work}, {"cumulated_restitution_work", cumulated_restitution_work}}; for (auto &&kv : toto) { auto &name = kv.first; auto &value = kv.second; auto &f = getOutputAsArray(name); f.clear(); f.push_back(value); } } } /* -------------------------------------------------------------------------- */ template void StimulationLangevin::addForce(Cont &cont) { if (computeWork) resetOldWorkPerAtom(cont); this->frictionCorrection(cont); this->restitutionCorrection(cont); if (computeWork) gatherWork(cont); } /* -------------------------------------------------------------------------- */ template void __attribute__((noinline)) StimulationLangevin::frictionCorrection(Cont &cont) { // constexpr UInt Dim = Cont::Dim; Real conv1 = code_unit_system.mv_t2f / damp; UInt index = 0; friction_work = 0; auto &old_friction_per_atom = getOutputAsArray("old_friction_per_atom"); auto &friction_work_per_atom = getOutputAsArray("friction_work_per_atom"); for (auto &&at : cont) { Real gamma1 = -conv1 * at.mass(); auto &&vel = at.velocity(); Vector friction_term = gamma1 * vel; if (computeWork) { Vector old_friction_term = old_friction_per_atom.row(index); friction_work += (friction_term + old_friction_term).dot(vel); old_friction_per_atom.row(index) = friction_term; if (computeWorkPerAtom) friction_work_per_atom.row(index) += friction_work; // accountFrictionWork( // friction_term, at.velocity()[i], index); } at.force() += friction_term; ++index; } } /* -------------------------------------------------------------------------- */ template void __attribute__((noinline)) StimulationLangevin::restitutionCorrection(Cont &cont) { // constexpr UInt Dim = Cont::Dim; Real conv2 = sqrt(24. * temperature * boltzmann * code_unit_system.kT2fd * code_unit_system.m_tt2f_d / damp / timestep); restitution_work = 0; auto &old_restitution_per_atom = getOutputAsArray("old_restitution_per_atom"); auto &restitution_work_per_atom = getOutputAsArray("restitution_work_per_atom"); UInt index = 0; for (auto &&at : cont) { Real gamma2 = conv2 * sqrt(at.mass()); Vector restitution_term = gamma2 * 0.5 * Vector::Random(); if (computeWork) { auto &&vel = at.velocity(); Vector old_restitution_term = old_restitution_per_atom.row(index); restitution_work += (restitution_term + old_restitution_term).dot(vel); old_restitution_term = restitution_term; if (computeWorkPerAtom) restitution_work_per_atom[index] += restitution_work; // accountRestitutionWork( // restitution_term, at.velocity()[i], index); } at.force() += restitution_term; ++index; } } /* -------------------------------------------------------------------------- */ // template // inline void StimulationLangevin::accountFrictionWork(Real friction_term, // Real velocity, // UInt index) { // Real old_friction_term = old_friction_per_atom[index]; // friction_work += (friction_term + old_friction_term) * velocity; // old_friction_per_atom[index] = friction_term; // if (computeWorkPerAtom) // friction_work_per_atom[index] += friction_work; // } // /* -------------------------------------------------------------------------- // */ // template // inline void StimulationLangevin::accountRestitutionWork(Real // restitution_term, // Real velocity, // UInt index) { // Real old_restitution_term = old_restitution_per_atom[index]; // restitution_work += (restitution_term + old_restitution_term) * velocity; // old_restitution_per_atom[index] = restitution_term; // if (computeWorkPerAtom) // restitution_work_per_atom[index] += restitution_work; // } /* -------------------------------------------------------------------------- */ /* LMDESC LANGEVIN This stimulator applies a Langevin thermostat to domain. The equation of motion becomes: .. math:: m\ddot{x} = - \nabla \phi \underbrace{- \frac{m}{damp}{v}}_{f_{damp}} + \underbrace{\sqrt{\frac{k_{b}\,T\,m}{dt\,damp}}R}_{f_{rand}} with :math:`\phi` the classical potential for interatomic forces, :math:`T` the temperature of the heat bath, :math:`damp` the relaxation time, :math:`v` the velocity field, :math:`k_b` the boltzmann constant, :math:`dt` the used timestep, :math:`m` the particle mass and :math:`R` a vector of random numbers uniformly distributed between -1/2 and 1/2. Explicitely, this stimulator append :math:`f_{damp}` and :math:`f_{rand}` to the force field during stage PRE_STEP3. */ /* LMEXAMPLE STIMULATION thermostat LANGEVIN INPUT md TEMP 100 SEED 32 */ /* LMHERITANCE action_interface */ void StimulationLangevin::declareParams() { StimulationInterface::declareParams(); this->changeDefault("STAGE", PRE_STEP3); /* LMKEYWORD DAMP Set the damping parameter expressed in the units of time */ this->parseKeyword("DAMP", damp); /* LMKEYWORD TEMP Set the temperature desired to be set */ this->parseKeyword("TEMP", temperature); /* LMKEYWORD SEED Set the seed for the random generator */ this->parseKeyword("SEED", seed); /* LMKEYWORD TIMESTEP Set the time step needed for random noise contruction */ this->parseKeyword("TIMESTEP", timestep); /* LMKEYWORD WORK_PERATOM TODO */ this->parseTag("WORK_PERATOM", compute_work_peratom_flag, false); /* LMKEYWORD WORK Activates the computation of the work done by ther thermostat. A compute named friction\_restitution:ID is created and allows the visualization of the work done by the thermostat (ID should be the name given to this compute). More details available - The work is computed like:\ \ - $$dW^n = \int_{x^n}^{x^{n+1}} f \cdot dx$$ + The work is computed like: + .. math:: + + dW^n = \int_{x^n}^{x^{n+1}} f \cdot dx + which for the case of the velocity verlet scheme is equivalent to - $$dW^n = \frac{1}{2}(f^n + f^{n+1}) v^{n+1/2}$$ + + .. math:: + dW^n = \frac{1}{2}(f^n + f^{n+1}) v^{n+1/2} + Thus when this flag is activated a compute with - with the name "friction\_restitution:STIMULATOR\_NAME" + with the name "friction_restitution:STIMULATOR_NAME" and 5 entries is registered. These entries are: - \begin{enumerate} - \item Requested temperature - \item frictional\_work: $dW_{damp}^n = \frac{1}{2}(f_{damp}^n + - f_{damp}^{n+1}) v^{n+1/2}$ - \item restitution\_work: $dW_{restitution}^n = \frac{1}{2}(f_{rand}^n + - f_{rand}^{n+1}) v^{n+1/2}$ - \item cumulated\_frictional\_work: $W_{damp}^n = \sum_0^n dW_{damp}^n$ - \item cumulated\_restitution\_work: $W_{restitution}^n = \sum_0^n - dW_{restitution}^n$ - \end{enumerate} + + - Requested temperature + - frictional_work: :math:`dW_{damp}^n = \frac{1}{2}(f_{damp}^n + + f_{damp}^{n+1}) v^{n+1/2}` + - restitution_work: :math:`dW_{restitution}^n = \frac{1}{2}(f_{rand}^n + + f_{rand}^{n+1}) v^{n+1/2}` + - cumulated_frictional_work: :math:`W_{damp}^n = \sum_0^n dW_{damp}^n` + - cumulated_restitution_work: :math:`W_{restitution}^n = \sum_0^n + dW_{restitution}^n` */ this->parseTag("WORK", compute_work_flag, false); } /* -------------------------------------------------------------------------- */ void StimulationLangevin::init() { MathTools::setSeed(seed); } /* -------------------------------------------------------------------------- */ DECLARE_STIMULATION_MAKE_CALL(StimulationLangevin) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_zero.cc b/src/stimulation/stimulation_zero.cc index acbddc5..771e5e8 100644 --- a/src/stimulation/stimulation_zero.cc +++ b/src/stimulation/stimulation_zero.cc @@ -1,153 +1,151 @@ /** * @file stimulation_raz.cc * * @author Guillaume Anciaux * @author Till Junge * * @date Wed Jul 09 21:59:47 2014 * * @brief Remise a zero (reset to zero) stimulator * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "stimulation_zero.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "lm_communicator.hh" #include "ref_point_data.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ StimulationZero::StimulationZero(const std::string &name) : LMObject(name), COEF(0), average_flag(false), direction(-1) {} /* -------------------------------------------------------------------------- */ StimulationZero::~StimulationZero() {} /* -------------------------------------------------------------------------- */ template void StimulationZero::stimulate(Cont &cont) { constexpr UInt Dim = Cont::Dim; std::vector directions; if (this->direction < 0) { for (UInt i = 0; i < Dim; ++i) { directions.push_back(i); } } else { directions.push_back(this->direction); } typedef std::vector::iterator iter; DUMP("Stimulation ZERO fields " << " COEF : " << COEF << " pointer this " << this, DBG_INFO); for (UInt i = 0; i < field.size(); ++i) { FieldType f = field[i]; for (auto &&at : cont) { if (f == _velocity) { for (iter i = directions.begin(); i != directions.end(); ++i) { at.velocity()[*i] += (COEF - 1) * at.velocity()[*i]; } } else if (f == _displacement) { for (iter i = directions.begin(); i != directions.end(); ++i) { throw; // at.displacement()[*i] += (COEF-1) * at.displacement()[*i]; } } else if (f == _force) { for (iter i = directions.begin(); i != directions.end(); ++i) { at.force()[*i] += (COEF - 1) * at.force()[*i]; } } else if (f == _position0) { for (iter i = directions.begin(); i != directions.end(); ++i) { at.position0()[*i] = at.position()[*i]; } } } } } /* -------------------------------------------------------------------------- */ /* LMDESC ZERO Remise a zero (reset to zero) stimulator\\ This stimulator reset some fields to zero or even rescale them according to the options provided. */ /* LMEXAMPLE STIMULATION reset ZERO INPUT md FIELD displacement FIELD velocity */ /* LMHERITANCE action_interface */ void StimulationZero::declareParams() { StimulationInterface::declareParams(); /* LMKEYWORD FIELD Set the fields that should be altered. The value selected can be: - \begin{itemize} - \item displacement : selects the displacement field - \item velocity : selects the velocity field - \item force : selects the force field - \item position0 : selects the initial position field. Actually it is - relevant for - atomic systems assigning the initial positions to the current one. \\ - This is useful when a restart file needs to be generated with current - positions - and a zero displacement field. - \end{itemize} + + - displacement : selects the displacement field + - velocity : selects the velocity field + - force : selects the force field + - position0 : selects the initial position field. Actually it is + relevant for atomic systems assigning the initial positions to the + current one. This is useful when a restart file needs to be generated + with current positions and a zero displacement field. + This command can be called several times for each field that need to be altered. */ this->parseKeyword("FIELD", field); /* LMKEYWORD COEF This is effective only for displacement, velocity or force field. Instead of reset to zero it multiplies by the provided coefficient. It enables, for instance to perform velocity damping. */ this->parseKeyword("COEF", COEF, 0.); /* LMKEYWORD COMPONENT Just apply stimulation in direction COMPONENT. By default all directions are stimulated. */ this->parseKeyword("COMPONENT", direction, -1); } /* -------------------------------------------------------------------------- */ DECLARE_STIMULATION_MAKE_CALL(StimulationZero) __END_LIBMULTISCALE__