diff --git a/.gitignore b/.gitignore index dbce7f0..4306c4a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,11 @@ *~ build* *.orig third-party/pybind11 +third-party/eigen +third-party/lammps +third-party/akantu __pycache__ TAGS -.ipynb_checkpoints \ No newline at end of file +.ipynb_checkpoints +#* diff --git a/doc/source/development.rst b/doc/source/development.rst new file mode 100644 index 0000000..2f79ad3 --- /dev/null +++ b/doc/source/development.rst @@ -0,0 +1,2 @@ +Developer Guide +======================================= diff --git a/doc/source/index.rst b/doc/source/index.rst index d7bd113..422efd4 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -1,41 +1,42 @@ .. image:: images/logo.svg :width: 200 :target: https://gitlab.com/libmultiscale/libmultiscale LibMultiScale manual ==================== In simulations, particle approaches can be helpful when the discreteness of matter needs to be taken into account. Multiscale coupling methods allow to reduce the prohibitive computational costs of discrete methods. For instance, with these approaches one can couple an atomic description with a macroscopic model of continuum mechanics. The LibMultiScale is a C++ parallel framework for the multiscale coupling methods dedicated to material simulations. This framework is designed as a library providing an API which makes it possible to program coupled simulations. The coupled parts can be provided by existing projects. In such a manner, the API gives C++ templated interfaces to reduce to the maximum the cost of integration taking the form of plugins or alike. LAMMPS (Sandia laboratories) and Akantu (LSMS) have been integrated to provide a functional framework. The LibMultiScale is now distributed with a joint :ref:`CECILL-C` and :ref:`LGPL` open-source licence with a shared copyrights between INRIA, CEA and EPFL. Table of content ================ .. toctree:: :maxdepth: 2 licences Installation Conffile Python keywords api + dev diff --git a/examples/atc_xiao_Al/material.dat b/examples/atc_xiao_Al/material.dat index 1f49155..7855d4d 100644 --- a/examples/atc_xiao_Al/material.dat +++ b/examples/atc_xiao_Al/material.dat @@ -1,9 +1,14 @@ material elastic_orthotropic [ name = aluminum-film-0K-LJ-orthotropic rho = 3.8328782866680342684 # density in g/mol/A^2 E1 = 9.7874197348334701168 # young's in kcal/mol/A^3 E2 = 9.7874197348334701168 # young's in kcal/mol/A^3 + E3 = 9.7874197348334701168 nu12 = 0.33333333333333333333 # poisson's ratio + nu13 = 0.33333333333333333333 + nu23 = 0.33333333333333333333 G12 = 3.6702824005625512938 # shear modulus in kcal/mol/A^3 + G13 = 3.6702824005625512938 + G23 = 3.6702824005625512938 ] diff --git a/src/common/lm_python_types.hh b/src/common/lm_python_types.hh index 9ceaa4e..a0471d3 100644 --- a/src/common/lm_python_types.hh +++ b/src/common/lm_python_types.hh @@ -1,90 +1,91 @@ #ifndef __LM_PYTHON_TYPES_HH__ #define __LM_PYTHON_TYPES_HH__ /* -------------------------------------------------------------------------- */ #include "lm_communicator.hh" #include "lm_types.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ inline void declare_types(py::module &m) { py::enum_(m, "IntegrationSchemeStage") .value("NONE_STEP", NONE_STEP) .value("PRE_DUMP", PRE_DUMP) .value("PRE_STEP1", PRE_STEP1) .value("PRE_STEP2", PRE_STEP2) .value("PRE_STEP3", PRE_STEP3) .value("PRE_STEP4", PRE_STEP4) .value("PRE_STEP5", PRE_STEP5) .value("PRE_STEP6", PRE_STEP6) .value("PRE_STEP7", PRE_STEP7) .value("PRE_FATAL", PRE_FATAL) .value("ALL_STEP", ALL_STEP) .export_values(); py::enum_(m, "CouplingStage") .value("COUPLING_STEP1", COUPLING_STEP1) .value("COUPLING_STEP2", COUPLING_STEP2) .value("COUPLING_STEP3", COUPLING_STEP3) .value("COUPLING_STEP4", COUPLING_STEP4) .value("COUPLING_STEP5", COUPLING_STEP5) .value("COUPLING_STEP6", COUPLING_STEP6) .value("COUPLING_STEP7", COUPLING_STEP7) .export_values(); py::enum_(m, "DOFType") .value("dt_master", dt_master) .value("dt_slave", dt_slave) .value("dt_pbc_master", dt_pbc_master) .value("dt_pbc_slave", dt_pbc_slave) .value("dt_all", dt_all) .export_values(); py::enum_(m, "FieldType") .value("position0", _position0) .value("position", _position) .value("displacement", _displacement) .value("velocity", _velocity) .value("force", _force) .value("stress", _stress) .value("temperature", _temperature) .value("grouprank", _grouprank) .value("strain", _strain) .value("epot", _epot) .value("applied_force", _applied_force) .value("mass", _mass) .value("tag", _tag) .value("id", _id) .value("proc", _proc) .value("charge", _charge) .value("burgers", _burgers) .value("normal", _normal) .value("boundary", _boundary) + .value("torque", _torque) .value("angular_velocity", _angular_velocity) .value("angular_acceleration", _angular_acceleration) .export_values(); py::class_(m, "CommGroup") .def(py::init()); py::class_(m, "Communicator") .def("getCommunicator", &Communicator::getCommunicator, py::return_value_policy::reference) .def("getNumberFreeProcs", &Communicator::getNumberFreeProcs) .def("addGroup", &Communicator::addGroup) .def("getGroup", &Communicator::getGroup) .def("reset", &Communicator::reset); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif //__LM_PYTHON_TYPES_HH__ diff --git a/src/common/possible_types.hh b/src/common/possible_types.hh index 191578d..34bda2f 100644 --- a/src/common/possible_types.hh +++ b/src/common/possible_types.hh @@ -1,255 +1,257 @@ #include "auto_arguments.hh" #include "lm_common.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ template using enable_if_type = AutoDispatch::enable_if_type; template using enable_if_not_type = AutoDispatch::enable_if_not_type; template using _or = AutoDispatch::_or; template using _params = AutoDispatch::_params; /* -------------------------------------------------------------------------- */ template class ContainerArray; template class ContainerMesh; /* -------------------------------------------------------------------------- */ // TODO : need to boost these declarations /* -------------------------------------------------------------------------- */ // MD /* -------------------------------------------------------------------------- */ // possible containers of MD #ifdef LIBMULTISCALE_MD1D class ContainerMD1D; class RefMD1D; class DomainMD1D; #endif #ifdef LIBMULTISCALE_LAMMPS_PLUGIN template class ContainerLammps; template class RefLammps; template class DomainLammps; #endif using md_sub_containers = _or< #ifdef LIBMULTISCALE_MD1D ContainerArray #endif #if defined(LIBMULTISCALE_LAMMPS_PLUGIN) && defined(LIBMULTISCALE_MD1D) , #endif #ifdef LIBMULTISCALE_LAMMPS_PLUGIN ContainerArray>, ContainerArray> #endif >; using md_containers = _or<_or< #ifdef LIBMULTISCALE_MD1D ContainerMD1D #endif #if defined(LIBMULTISCALE_LAMMPS_PLUGIN) && defined(LIBMULTISCALE_MD1D) , #endif #ifdef LIBMULTISCALE_LAMMPS_PLUGIN ContainerLammps<2>, ContainerLammps<3> #endif >, md_sub_containers>; using md_domains = _or< #ifdef LIBMULTISCALE_MD1D DomainMD1D #endif #if defined(LIBMULTISCALE_LAMMPS_PLUGIN) && defined(LIBMULTISCALE_MD1D) , #endif #ifdef LIBMULTISCALE_LAMMPS_PLUGIN DomainLammps<2>, DomainLammps<3> #endif >; using subMD = md_sub_containers; using MD = md_containers; using domMD = md_domains; /* -------------------------------------------------------------------------- */ // DD /* -------------------------------------------------------------------------- */ #ifdef LIBMULTISCALE_PARADIS_PLUGIN class ContainerNodesParaDiS; class ContainerElemsParaDiS; class RefNodeParaDiS; class RefElemParaDiS; class DomainParaDiS; #endif // possible containers of DD using dd_sub_containers = _or< #ifdef LIBMULTISCALE_PARADIS_PLUGIN ContainerMesh, ContainerArray> #endif >; using dd_containers = _or< #ifdef LIBMULTISCALE_PARADIS_PLUGIN ContainerMesh #endif >; using dd_domains = _or< #ifdef LIBMULTISCALE_PARADIS_PLUGIN DomainParaDiS #endif >; using subDD = dd_sub_containers; using DD = dd_containers; using domDD = dd_domains; /* -------------------------------------------------------------------------- */ // CONTINUUM /* -------------------------------------------------------------------------- */ // possible containers of CONTINUUM #ifdef LIBMULTISCALE_MECA1D class ContainerElemsMeca1D; class ContainerNodesMeca1D; class RefNodeMeca1D; class RefElemMeca1D; class DomainMeca1D; #endif #ifdef LIBMULTISCALE_AKANTU_PLUGIN template class ContainerElemsAkantu; template class ContainerNodesAkantu; template class RefNodeAkantu; template class RefElemAkantu; template class DomainAkantu; #endif using continuum_sub_containers = std::tuple< #ifdef LIBMULTISCALE_MECA1D ContainerMesh, ContainerArray> #endif #if defined(LIBMULTISCALE_MECA1D) && defined(LIBMULTISCALE_AKANTU_PLUGIN) , #endif #ifdef LIBMULTISCALE_AKANTU_PLUGIN ContainerMesh>, ContainerArray>>, ContainerMesh>, ContainerArray>>, ContainerMesh>, ContainerArray>> #endif >; using continuum_containers = _or #endif #if defined(LIBMULTISCALE_MECA1D) && defined(LIBMULTISCALE_AKANTU_PLUGIN) , #endif #ifdef LIBMULTISCALE_AKANTU_PLUGIN ContainerMesh, ContainerElemsAkantu<1>>, ContainerMesh, ContainerElemsAkantu<2>>, ContainerMesh, ContainerElemsAkantu<3>> #endif >, continuum_sub_containers>; using continuum_domains = std::tuple< #ifdef LIBMULTISCALE_MECA1D DomainMeca1D #endif #if defined(LIBMULTISCALE_MECA1D) && defined(LIBMULTISCALE_AKANTU_PLUGIN) , #endif #ifdef LIBMULTISCALE_AKANTU_PLUGIN DomainAkantu<1>, DomainAkantu<2>, DomainAkantu<3> #endif >; using subCONTINUUM = continuum_sub_containers; using CONTINUUM = continuum_containers; using domCONTINUUM = continuum_domains; /* -------------------------------------------------------------------------- */ // GenericContainers /* -------------------------------------------------------------------------- */ template class RefPoint; template class RefPointData; template class RefGenericElem; template using ContainerGenericMesh = ContainerMesh>, ContainerArray>>; template class RefGenericDDElem; template class RefGenericDDNode; template using ContainerGenericDDMesh = ContainerMesh>, ContainerArray>>; // possible containers of points using point_containers = _or>, ContainerArray>, ContainerArray>>; using POINT = point_containers; // possible containers of meshes using mesh_containers = _or, ContainerGenericMesh<2u>, ContainerGenericMesh<3u>>; using MESH = mesh_containers; // possible containers of DD using dd_mesh_containers = _or>; using DDMESH = dd_mesh_containers; // possible computes as input class Component; using COMPONENT = _or; using ARRAY = _or>; /* -------------------------------------------------------------------------- */ template using enable_if_md = enable_if_type<_params, MD>; template +using enable_if_not_md = enable_if_not_type<_params, MD>; +template using enable_if_mesh = enable_if_type<_params, _or>; template using enable_if_continuum = enable_if_type<_params, CONTINUUM>; template using enable_if_point = enable_if_type<_params, _or>; template using enable_if_dd = enable_if_type<_params, DD>; template using enable_if_component = enable_if_type<_params, COMPONENT>; template using enable_if_not_component = enable_if_not_type<_params, COMPONENT>; /* --------------------------------------------------------------------- */ __END_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ diff --git a/src/dumper/dumper_restart.cc b/src/dumper/dumper_restart.cc index e8660ed..07ac239 100644 --- a/src/dumper/dumper_restart.cc +++ b/src/dumper/dumper_restart.cc @@ -1,348 +1,348 @@ /** * @file dumper_restart.cc * * @author Guillaume Anciaux * * @date Wed Aug 20 16:58:08 2014 * * @brief This dumper saves the state of the simulation * * @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 "dumper_restart.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 #include /* ------------------------------------------------------------------------ */ __BEGIN_LIBMULTISCALE__ /* ------------------------------------------------------------------------ */ //! local number of Dofs per proc. Only valid for root (rank = 0) processor static std::vector nb_local_per_proc; //! buffer data used to gather information on the one that will write to disk. //! Only valid for root static std::vector data_per_proc; /* ------------------------------------------------------------------------ */ DumperRestart::DumperRestart(const std::string &name) : LMObject(name) { nb_local_dofs = 0; allocated_size = 0; data = NULL; text_flag = false; } /* ------------------------------------------------------------------------ */ DumperRestart::~DumperRestart() { // if (lm_my_proc_id == 0){ // //liberation de la memoire que je veux pas maintenir en parallel // DUMP("freeing data_per_proc",DBG_INFO); // for (UInt i = 1 ; i < lm_world_size ; ++i){ // DUMP("freeing ptr[" << i << "] " << data_per_proc[i],DBG_INFO); // if (data_per_proc[i] != NULL){ // delete (data_per_proc[i]); // data_per_proc[i] = NULL; // } // } // nb_local_per_proc.clear(); // data_per_proc.clear(); // requests.clear(); // } // if (data) delete(data); } /* ------------------------------------------------------------------------ */ void DumperRestart::initVectors(CommGroup &group) { UInt my_rank = group.getMyRank(); UInt group_size = group.size(); if (my_rank == 0) { nb_local_per_proc.resize(group_size); data_per_proc.resize(group_size); nb_local_per_proc.assign(group_size, 0); data_per_proc.assign(group_size, 0); for (UInt i = 0; i < group_size; ++i) data_per_proc[i] = NULL; } } /* ------------------------------------------------------------------------ */ template void DumperRestart::dump(Cont &cont) { CommGroup &group = cont.getCommGroup(); this->initVectors(group); UInt my_rank = group.getMyRank(); UInt group_size = group.size(); constexpr UInt Dim = Cont::Dim; LMFile file; UInt nbTotalDofs = 0; if (my_rank == 0) { std::stringstream temp; temp << this->getBaseName() << "_restart-" << std::setfill('0') << std::setw(4) << this->action_step << ".xml"; file.open(temp.str(), "wb", true); b64.setOutputFile(file); } bool need_reallocate = false; if (allocated_size < cont.size() && cont.size() != 0) need_reallocate = true; nb_local_dofs = cont.size(); group.gather(&nb_local_dofs, 1, &nb_local_per_proc[0], 0, "gather nb local dofs"); DUMP("Gather done", DBG_INFO); if (my_rank == 0) { for (UInt i = 0; i < group_size; ++i) { DUMP("nb_local[" << i << "]=" << nb_local_per_proc[i] << " " << nb_local_dofs, DBG_MESSAGE); nbTotalDofs += nb_local_per_proc[i]; } printHeaders(nbTotalDofs, Dim, file); } // allocation des array de stockage if (need_reallocate) { allocated_size = nb_local_dofs; DUMP("(re)allocating for size = " << nb_local_dofs, DBG_INFO); if (data != NULL) data = (Real *)realloc(data, nb_local_dofs * sizeof(Real) * Dim); else { data = (Real *)malloc(nb_local_dofs * sizeof(Real) * Dim); memset(data, 0, sizeof(Real) * nb_local_dofs * Dim); } } if (my_rank == 0) { // je switch le pointeur pour les donnees locales data_per_proc[0] = data; for (UInt i = 1; i < group_size; ++i) { if (nb_local_per_proc[i] == 0) continue; if (data_per_proc[i] != NULL) data_per_proc[i] = (Real *)realloc( data_per_proc[i], sizeof(Real) * nb_local_per_proc[i] * Dim); else { data_per_proc[i] = (Real *)malloc(sizeof(Real) * nb_local_per_proc[i] * Dim); memset(data_per_proc[i], 0, sizeof(Real) * nb_local_per_proc[i] * Dim); } } } // chaqun construit en parallele les data qu'il veut envoyer UInt cpt = 0; for (auto &&at : cont) { auto X = at.position0(); for (UInt i = 0; i < Dim; ++i) data[cpt * Dim + i] = X[i]; ++cpt; } DUMP("a la fin cpt = " << cpt, DBG_INFO); DUMP("local P0 construction done", DBG_INFO); DumpField("P0", file, group); cpt = 0; for (auto &&at : cont) { auto U = at.displacement(); for (UInt i = 0; i < Dim; ++i) data[cpt * Dim + i] = U[i]; ++cpt; } DUMP("local U construction done", DBG_INFO); DumpField("U", file, group); cpt = 0; for (auto &&at : cont) { auto V = at.velocity(); for (UInt i = 0; i < Dim; ++i) data[cpt * Dim + i] = V[i]; ++cpt; } DUMP("local V construction done", DBG_INFO); DumpField("V", file, group); // on fait le champ de position P0 // je rassemble le tout en local if (my_rank == 0) { printTail(file); file.close(); } group.synchronize(); } /* ------------------------------------------------------------------------ */ template inline void DumperRestart::DumpField(const std::string &fieldName, LMFile &file, CommGroup &group) { UInt group_size = group.size(); UInt my_rank = group.getMyRank(); if (my_rank == 0) { // I switch the pointer for local data for (UInt i = 1; i < group_size; ++i) { if (nb_local_per_proc[i] == 0) continue; DUMP("pending recv from " << i << " of size " << nb_local_per_proc[i] << " ...", DBG_INFO); group.receive(data_per_proc[i], nb_local_per_proc[i] * Dim, i, "restart data per proc"); DUMP("pending recv from " << i << " of size " << nb_local_per_proc[i] << " ... done", DBG_INFO); } } else if (nb_local_dofs != 0) { DUMP("sending to root " << nb_local_dofs << " ...", DBG_INFO); group.send(data, nb_local_dofs * Dim, 0, "restart data per proc"); DUMP("sending to root " << nb_local_dofs << " ... done", DBG_INFO); } if (my_rank == 0) { if (!text_flag) b64.clearBuffer(); file.printf("<%s>\n", fieldName.c_str()); for (int i = group_size - 1; i > -1; --i) { DUMP("waiting com from " << i, DBG_INFO); // attend les coms avant d'ecrire sur disque if (i > 0 && nb_local_per_proc[i] > 0) { LM_TOIMPLEMENT; } DUMP("reception from " << i << " complete", DBG_INFO); for (UInt j = 0; j < nb_local_per_proc[i]; ++j) { if (text_flag) { for (UInt k = 0; k < Dim; ++k) { file.printf("%.15e", data_per_proc[i][j * Dim + k]); if (k != Dim - 1) file.printf("\t"); } file.printf("\n"); } else { for (UInt k = 0; k < Dim; ++k) { b64.pushRealInBase64(data_per_proc[i][j * Dim + k]); } } } } if (!text_flag) { b64.finish(); b64.dumpToFile(); } file.printf("\n", fieldName.c_str()); } } /* ------------------------------------------------------------------------ */ void DumperRestart::printHeaders(UInt nbDofs, const UInt Dim, LMFile &file) { file.printf("\n"); file.printf("\n"); else file.printf(" dataType=\"BINARY\">\n"); } /* ------------------------------------------------------------------------ */ void DumperRestart::printTail(LMFile &file) { file.printf("\n"); } /* ------------------------------------------------------------------------ */ /* LMDESC RESTART This dumper saves the state of the simulation in order to allow - restarting simulations from + restarting simulations from this state. */ /* LMEXAMPLE DUMPER res RESTART INPUT md PREFIX /home/titeuf */ /* LMHERITANCE dumper_interface */ void DumperRestart::declareParams() { DumperInterface::declareParams(); /* LMKEYWORD TEXT Flag to request human-readable text output instead of binary. */ this->parseTag("TEXT", text_flag, false); } /* ------------------------------------------------------------------------ */ DECLARE_DUMPER_MAKE_CALL(DumperRestart) __END_LIBMULTISCALE__ diff --git a/src/dumper/dumper_restart.hh b/src/dumper/dumper_restart.hh index c0cddd8..2b8c439 100644 --- a/src/dumper/dumper_restart.hh +++ b/src/dumper/dumper_restart.hh @@ -1,99 +1,99 @@ /** * @file dumper_restart.hh * * @author Guillaume Anciaux * * @date Tue Oct 29 22:09:05 2013 * * @brief This dumper saves the state of the simulation * * @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. * */ #ifndef __LIBMULTISCALE_DUMPER_RESTART_HH__ #define __LIBMULTISCALE_DUMPER_RESTART_HH__ /* -------------------------------------------------------------------------- */ #include "base64_writer.hh" #include "dumper_interface.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ class DumperRestart : public DumperInterface { public: - DECLARE_DUMPER(DumperRestart, _or); + DECLARE_DUMPER(DumperRestart, _or); //! most generic dump function template void dump(_Input &cont); void initVectors(CommGroup &group); protected: //! prUInt the headers in a given file for the xml restart format void printHeaders(UInt nbDofs, const UInt Dim, LMFile &file); //! prUInt the headers in a given file for the xml restart format void printTail(LMFile &file); //! make coms and dump to file with little communication overlapping template void DumpField(const std::string &fieldName, LMFile &file, CommGroup &group); protected: //! local number of dofs UInt nb_local_dofs; //! maintain info of allocated buffer sizes UInt allocated_size; //! local data Real *data; //! flag to know if a binary or text dump will be done bool text_flag; //! helper for base64 Base64Writer b64; }; /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_DUMPER_RESTART_HH__ */ diff --git a/src/filter/compute_periodic_position.cc b/src/filter/compute_periodic_position.cc index 89b7c47..2a80127 100644 --- a/src/filter/compute_periodic_position.cc +++ b/src/filter/compute_periodic_position.cc @@ -1,136 +1,136 @@ /** * @file compute_periodic_position.cc * * @author Guillaume Anciaux * @author Till Junge * * @date Wed Jul 09 21:59:47 2014 * * @brief This computes periodic positions by correcting for PBC * * @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 . * */ #define TIMER /* -------------------------------------------------------------------------- */ #include "compute_periodic_position.hh" #include "compute_extract.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_point_data.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ ComputePeriodicPosition::ComputePeriodicPosition(const std::string &name) : LMObject(name), ComputeTrueDisplacement(name) {} /* -------------------------------------------------------------------------- */ ComputePeriodicPosition::~ComputePeriodicPosition() {} /* -------------------------------------------------------------------------- */ template void ComputePeriodicPosition::build(Cont &cont) { constexpr UInt Dim = Cont::Dim; if (this->pbc[0] == 0 && this->pbc[1] == 0 && this->pbc[2] == 0) buildWithPBC(cont); if (this->pbc[0] == 1 && this->pbc[1] == 0 && this->pbc[2] == 0) buildWithPBC(cont); if (this->pbc[0] == 1 && this->pbc[1] == 1 && this->pbc[2] == 0) buildWithPBC(cont); if (this->pbc[0] == 1 && this->pbc[1] == 1 && this->pbc[2] == 1) buildWithPBC(cont); if (this->pbc[0] == 1 && this->pbc[1] == 0 && this->pbc[2] == 1) buildWithPBC(cont); if (this->pbc[0] == 0 && this->pbc[1] == 1 && this->pbc[2] == 0) buildWithPBC(cont); if (this->pbc[0] == 0 && this->pbc[1] == 1 && this->pbc[2] == 1) buildWithPBC(cont); if (this->pbc[0] == 0 && this->pbc[1] == 0 && this->pbc[2] == 1) buildWithPBC(cont); } /* -------------------------------------------------------------------------- */ template void ComputePeriodicPosition::buildWithPBC(Cont &cont) { this->clear(); Cube &bbox = cont.getBoundingBox(); auto d_size = bbox.getSize(); const Vector<3> &xmin = bbox.getXmin(); DUMP("here is my geom " << bbox, DBG_DETAIL); DUMP("here is interesting info " << d_size[0] << " " << d_size[1] << " " << d_size[2] << " " << xmin[0] << " " << xmin[1] << " " << xmin[2], DBG_DETAIL); UInt cpt = 0; for (auto &&at : cont) { Vector position = at.position(); for (UInt i = 0; i < Dim; ++i) { if ((pbc_x && i == 0) || (pbc_y && i == 1) || (pbc_z && i == 2)) { int offset = floor((position[i] - xmin[i]) / d_size[i]); position[i] -= offset * d_size[i]; } this->getOutputAsArray().push_back(position[i]); } ++cpt; } } /* -------------------------------------------------------------------------- */ /* LMDESC PERIODICPOSITION TODO */ /* LMEXAMPLE - COMPUTE periodicDisp PERIODICDISP INPUT md PBC 1 1 0 + COMPUTE periodicPosition PERIODICPOSITION INPUT md PBC 1 1 0 */ /* LMHERITANCE compute_true_displacement */ void ComputePeriodicPosition::declareParams() { ComputeTrueDisplacement::declareParams(); } /* -------------------------------------------------------------------------- */ DECLARE_COMPUTE_MAKE_CALL(ComputePeriodicPosition) __END_LIBMULTISCALE__ diff --git a/src/md/1d/ref_atom_md1d.hh b/src/md/1d/ref_atom_md1d.hh index 012a70c..b9250aa 100644 --- a/src/md/1d/ref_atom_md1d.hh +++ b/src/md/1d/ref_atom_md1d.hh @@ -1,208 +1,216 @@ /** * @file ref_atom_md1d.hh * * @author Guillaume Anciaux * @author Srinivasa Babu Ramisetti * * @date Fri Jan 10 20:47:45 2014 * * @brief This is the reference over atoms handled by the domainMD1D * * @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. * */ #ifndef __LIBMULTISCALE_REF_ATOM_MD1D_HH__ #define __LIBMULTISCALE_REF_ATOM_MD1D_HH__ /* -------------------------------------------------------------------------- */ #include "ref_atom.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ class ContainerMD1D; class ComparatorRefMD1D; class DomainMD1D; /* -------------------------------------------------------------------------- */ /** * Class RefMD1D * */ class RefMD1D : public RefAtom<1, RefMD1D> { public: typedef ComparatorRefMD1D RefComparator; typedef DomainMD1D Domain; friend class IteratorMD1D; friend class ContainerMD1D; friend class ComparatorRefMD1D; static constexpr UInt Dim = 1; using fields = field_list<_position0, _position, _displacement, _velocity, _force, _stress, _mass, _acceleration>; RefMD1D(); RefMD1D(ContainerMD1D &dom, UInt index = 0); bool operator==(RefMD1D &ref) { return (ref.index_atom == index_atom); }; VectorView position0(); VectorView position(); AccessorAtomDof displacement(); VectorView velocity(); VectorView force(); VectorView acceleration(); + VectorView torque(); + VectorView angular_velocity(); Real NormContrainte() { return 0; }; Real &mass(); void setMass(Real mass); UInt typeID() { return 0; }; Real &charge(); UInt id(); UInt tag(); void setIndex(UInt n); UInt getIndex(); Real getEPot(); private: Real *_m; Real *_charge; Real *_p0; Real *_p; Real *_v; Real *_f; Real *_a; Real *_pe; - // currently pointed atom UInt index_atom; }; /* -------------------------------------------------------------------------- */ inline RefMD1D::RefMD1D() : _m(nullptr), _charge(nullptr), _p0(nullptr), _p(nullptr), _v(nullptr), _f(nullptr), _a(nullptr), _pe(nullptr), index_atom(0) {} /* -------------------------------------------------------------------------- */ inline void RefMD1D::setIndex(UInt n) { index_atom = n; } /* -------------------------------------------------------------------------- */ inline UInt RefMD1D::getIndex() { return index_atom; } /* -------------------------------------------------------------------------- */ inline UInt RefMD1D::id() { return index_atom; } /* -------------------------------------------------------------------------- */ inline void RefMD1D::setMass(Real mass) { _m[index_atom] = mass; } /* -------------------------------------------------------------------------- */ inline Real &RefMD1D::charge() { return _charge[index_atom]; } /* -------------------------------------------------------------------------- */ inline UInt RefMD1D::tag() { return index_atom; } /* -------------------------------------------------------------------------- */ inline VectorView<1> RefMD1D::position0() { return VectorView<1>(_p0 + index_atom); } /* -------------------------------------------------------------------------- */ inline VectorView<1> RefMD1D::position() { return VectorView<1>(_p + index_atom); } /* -------------------------------------------------------------------------- */ inline AccessorAtomDof RefMD1D::displacement() { return AccessorAtomDof(*this); } /* -------------------------------------------------------------------------- */ inline VectorView<1> RefMD1D::velocity() { return VectorView<1>(_v + index_atom); } /* -------------------------------------------------------------------------- */ +inline VectorView<1> RefMD1D::torque() { LM_TOIMPLEMENT; } +/* -------------------------------------------------------------------------- */ + +inline VectorView<1> RefMD1D::angular_velocity() { LM_TOIMPLEMENT; } + +/* -------------------------------------------------------------------------- */ + inline VectorView<1> RefMD1D::force() { return VectorView<1>(_f + index_atom); } /* -------------------------------------------------------------------------- */ inline VectorView<1> RefMD1D::acceleration() { return VectorView<1>(_a + index_atom); } /* -------------------------------------------------------------------------- */ inline Real &RefMD1D::mass() { return _m[index_atom]; } /* -------------------------------------------------------------------------- */ inline Real RefMD1D::getEPot() { return _pe[index_atom]; } /* -------------------------------------------------------------------------- */ class ComparatorRefMD1D { public: bool operator()(RefMD1D &r1, RefMD1D &r2) const { return r1.index_atom < r2.index_atom; } bool operator()(const RefMD1D &r1, const RefMD1D &r2) const { return r1.index_atom < r2.index_atom; } }; /* -------------------------------------------------------------------------- */ inline std::ostream &operator<<(std::ostream &os, RefMD1D &ref) { UInt n = ref.getIndex(); os << "MD1D reference : " << n; return os; } __END_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ #endif /* __LIBMULTISCALE_REF_ATOM_MD1D_HH__ */ diff --git a/src/parser/lm_parser_inline_impl.hh b/src/parser/lm_parser_inline_impl.hh index a26f30b..a101c4a 100644 --- a/src/parser/lm_parser_inline_impl.hh +++ b/src/parser/lm_parser_inline_impl.hh @@ -1,386 +1,392 @@ /** * @file lm_parser_inline_impl.hh * * @author Guillaume Anciaux * * @date Thu Jul 24 14:21:58 2014 * * @brief This is the central parser for LM * * @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 "lm_parser.hh" #include // backward declaration namespace pybind11 { class object; class dict; } // namespace pybind11 namespace py = pybind11; __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ template inline UInt parseNumber(T &d, std::stringstream &line) { std::string buffer; UInt nb = Parser::strNext(buffer, line); Real res; bool alg_parsed = false; try { alg_parsed = Parser::math_parser.parse(buffer, res); } catch (std::string mess) { LM_FATAL(Parser::getParserState() << std::endl << "\"" << line.str() << "\"" << std::endl << mess); } if (!alg_parsed) { std::stringstream sbuf(buffer); sbuf >> res; if (sbuf.fail()) LM_FATAL(Parser::getParserState() << std::endl << "\"" << line.str() << "\"" << std::endl << "could not parse " << line.str().substr(line.tellg()) << " or " << buffer); } d = static_cast(res); return nb; } /* -------------------------------------------------------------------------- */ template <> inline UInt parseNumber(bool &d, std::stringstream &line) { try { std::string buffer; UInt nb = Parser::strNext(buffer, line); std::transform(buffer.begin(), buffer.end(), buffer.begin(), ::tolower); if (buffer == "true") d = true; else if (buffer == "false") d = false; else if (buffer == "1") d = true; else if (buffer == "0") d = false; else { d = !d; return 0; } return nb; } catch (...) { d = !d; } return 0; } /* -------------------------------------------------------------------------- */ inline UInt _parse(std::string &val, std::stringstream &line, UInt = 0) { return Parser::strNext(val, line); } /* -------------------------------------------------------------------------- */ template std::enable_if_t::value, UInt> _parse(T &val, std::stringstream &line, UInt n_entries = 0); /* -------------------------------------------------------------------------- */ template inline std::enable_if_t< std::is_arithmetic::value and not std::is_same::value, UInt> _parse(T &val, std::stringstream &line, UInt = 0) { return parseNumber(val, line); } /* -------------------------------------------------------------------------- */ class Component; UInt _parse(Component &comp, std::stringstream &line, UInt n_entries = 0); /* -------------------------------------------------------------------------- */ UInt _parse(std::map> val, std::stringstream &line, UInt n_entries = 0); /* -------------------------------------------------------------------------- */ inline UInt _parse(char &val, std::stringstream &line, UInt = 0) { std::string buffer; UInt nb = Parser::strNext(buffer, line); std::stringstream sbuf(buffer); sbuf >> val; if (sbuf.fail()) LM_FATAL(Parser::getParserState() << std::endl << "\"" << line.str() << "\"" << std::endl << "could not parse " << line.str().substr(line.tellg()) << " or " << buffer); return nb; } /* -------------------------------------------------------------------------- */ template inline UInt _parse(std::vector &vec, std::stringstream &line, UInt = 0) { T tmp; UInt nb_word = _parse(tmp, line); vec.push_back(tmp); return nb_word; } /* -------------------------------------------------------------------------- */ template inline UInt _parse_vector(T &val, std::stringstream &line, UInt n_entries) { UInt nb_word = 0; using type = typename std::remove_reference::type; type tmp = type(); if (n_entries > int(nb)) LM_FATAL("invalid vector size: " << n_entries << " > " << nb); UInt size = nb; if (n_entries > 0) size = n_entries; for (UInt i = 0; i < size; ++i) { nb_word += _parse(tmp, line); val[i] = tmp; } return nb_word; } /* -------------------------------------------------------------------------- */ template inline UInt _parse(Vector &val, std::stringstream &line, UInt n_entries = 0) { return _parse_vector(val, line, n_entries); } /* -------------------------------------------------------------------------- */ template inline UInt _parse(Quantity &val, std::stringstream &line, UInt n_entries = 0) { return _parse_vector(val, line, n_entries); } /* -------------------------------------------------------------------------- */ template inline UInt _parse(PhysicalScalar &val, std::stringstream &line, UInt = 0) { Real tmp; UInt nb = _parse(tmp, line); val = tmp * Parser::units_converter.getConversion(); return nb; } /* -------------------------------------------------------------------------- */ inline UInt _parse(IntegrationSchemeMask &mask, std::stringstream &line, UInt = 0) { IntegrationSchemeMask value; std::string read_stage; UInt nb = Parser::strNext(read_stage, line); if (read_stage == "PRE_DUMP") value = PRE_DUMP; else if (read_stage == "PRE_STEP1") value |= PRE_STEP1; else if (read_stage == "PRE_STEP2") value |= PRE_STEP2; else if (read_stage == "PRE_STEP3") value |= PRE_STEP3; else if (read_stage == "PRE_STEP4") value |= PRE_STEP4; else if (read_stage == "PRE_STEP5") value |= PRE_STEP5; else if (read_stage == "PRE_STEP6") value |= PRE_STEP6; else if (read_stage == "PRE_STEP7") value |= PRE_STEP7; else if (read_stage == "PRE_FATAL") value |= PRE_FATAL; else { std::stringstream sstr; sstr << read_stage; nb += parseNumber(value, sstr); } mask |= value; return nb; } /* -------------------------------------------------------------------------- */ inline UInt _parse(UnitSystem &usys, std::stringstream &line, UInt = 0) { std::string name; UInt nb = Parser::strNext(name, line); if (name == "SIUnits") usys = si_unit_system; else if (name == "RealUnits") - usys = real_unit_system; + usys = real_unit_system; else if (name == "MetalUnits") usys = metal_unit_system; else LM_FATAL("unknown Units system " << name); return nb; } /* -------------------------------------------------------------------------- */ template inline UInt _parse(T (&val)[nb], std::stringstream &line, UInt n_entries = 0) { return _parse_vector(val, line, n_entries); } /* -------------------------------------------------------------------------- */ inline UInt _parse(FieldType &ftype, std::stringstream &line, UInt = 0) { std::string name; UInt nb = Parser::strNext(name, line); if (name == "position0") ftype = _position0; else if (name == "position") ftype = _position; else if (name == "displacement") ftype = _displacement; else if (name == "velocity") ftype = _velocity; else if (name == "force") ftype = _force; else if (name == "temperature") ftype = _temperature; else if (name == "grouprank") ftype = _grouprank; else if (name == "stress") ftype = _stress; else if (name == "strain") ftype = _strain; else if (name == "epot") ftype = _epot; else if (name == "burgers") ftype = _burgers; else if (name == "normal") ftype = _normal; + else if (name == "torque") + ftype = _torque; + else if (name == "angular_velocity") + ftype = _angular_velocity; + else if (name == "angular_acceleration") + ftype = _angular_acceleration; else LM_FATAL("unknown field type " << name); return nb; } /* -------------------------------------------------------------------------- */ inline UInt _parse(Operator &op, std::stringstream &line, UInt = 0) { std::string name; UInt nb = Parser::strNext(name, line); if (name == "NONE") op = OP_NONE; else if (name == "AVERAGE") op = OP_AVERAGE; else if (name == "DENSITY") op = OP_DENSITY; else if (name == "DEVIATION") op = OP_DEVIATION; else if (name == "MAX") op = OP_MAX; else if (name == "MIN") op = OP_MIN; else if (name == "SUM") op = OP_SUM; else LM_FATAL("unknown operator " << name); return nb; } /* -------------------------------------------------------------------------- */ inline UInt _parse(LatticeType <ype, std::stringstream &line, UInt = 0) { std::string name; UInt nb = Parser::strNext(name, line); if (name == "mono_atom_1d") ltype = mono_atom_1d; else if (name == "hex_2d") ltype = hex_2d; else if (name == "hex_2d_radial") ltype = hex_2d_radial; else if (name == "fcc_3d") ltype = fcc_3d; else LM_FATAL("unknown lattice type " << name); return nb; } /* -------------------------------------------------------------------------- */ inline UInt _parse(DOFType &dt, std::stringstream &line, UInt = 0) { std::string name; UInt nb = Parser::strNext(name, line); if (name == "dt_master") dt = dt_master; else if (name == "dt_slave") dt = dt_slave; else if (name == "dt_pbc_master") dt = dt_pbc_master; else if (name == "dt_pbc_slave") dt = dt_pbc_slave; else if (name == "dt_all") dt = dt_all; else LM_FATAL("unknown DOFType " << name); return nb; } /* -------------------------------------------------------------------------- */ // inline UInt _parse(char &val, std::stringstream &line) { // std::string val_str; // UInt nb = Parser::strNext(val_str, line); // if (val_str.size() != 1) { // LM_THROW("takes only one character"); // } // val = val_str[0]; // return nb; // } /* -------------------------------------------------------------------------- */ template inline UInt Parser::parse(T &var, std::stringstream &line, UInt n_entries) { return _parse(var, line, n_entries); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_addforce.cc b/src/stimulation/stimulation_addforce.cc index d0fad46..8f612db 100644 --- a/src/stimulation/stimulation_addforce.cc +++ b/src/stimulation/stimulation_addforce.cc @@ -1,124 +1,124 @@ /** * @file stimulation_addforce.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief Add a force to a set of atoms/nodes/points * * @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_addforce.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" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ StimulationAddForce::StimulationAddForce(const std::string &name) : LMObject(name) { for (UInt i = 0; i < 3; ++i) { total_force[i] = 0.0; per_atom_force[i] = 0.0; } } /* -------------------------------------------------------------------------- */ StimulationAddForce::~StimulationAddForce() {} /* -------------------------------------------------------------------------- */ template void StimulationAddForce::stimulate(Cont &cont) { constexpr UInt Dim = Cont::Dim; UInt cpt = cont.size(); CommGroup &group = cont.getCommGroup(); UInt rank = group.getMyRank(); group.allReduce(&cpt, 1, "reduce total number of points concerned", OP_SUM); if (rank == 0) { for (UInt i = 0; i < Dim; ++i) DUMP("In direction " << i << " add an additional force " << total_force[i] / (Real)cpt + per_atom_force[i], DBG_INFO); } for (auto &&at : cont) { for (UInt i = 0; i < Dim; ++i) { at.force()[i] += total_force[i] / Real(cpt) + per_atom_force[i]; } } } /* -------------------------------------------------------------------------- */ /* LMDESC ADDFORCE Add a force to a set of atoms/nodes/points The stimulator is provided a force :math:`f^{req}` and for each atoms/nodes/points in the input the stimulator do: .. math:: f += \frac{1}{\#points} f^{total} + f^{per-atom-force} */ /* LMEXAMPLE .. code-block:: STIMULATION force ADDFORCE TOTAL_FORCE 1 STIMULATION force ADDFORCE PER_ATOM_FORCE 0.001 */ void StimulationAddForce::declareParams() { StimulationInterface::declareParams(); this->changeDefault("STAGE", PRE_STEP3); /* LMKEYWORD TOTAL_FORCE Give the total force vector requested to be added. This quantity will be divided by the number of atoms in the container. */ this->parseVectorKeyword("TOTAL_FORCE", spatial_dimension, total_force, VEC_DEFAULTS(0., 0., 0.)); /* LMKEYWORD PER_ATOM_FORCE Give the force vector requested to be added on each atom. */ this->parseVectorKeyword("PER_ATOM_FORCE", spatial_dimension, per_atom_force, VEC_DEFAULTS(0., 0., 0.)); } /* -------------------------------------------------------------------------- */ -DECLARE_STIMULATION_MAKE_CALL(StimulationAddForce) +DECLARE_STIMULATION_MAKE_CALL(StimulationAddForce, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_block.cc b/src/stimulation/stimulation_block.cc index 1f52044..aa85001 100644 --- a/src/stimulation/stimulation_block.cc +++ b/src/stimulation/stimulation_block.cc @@ -1,105 +1,105 @@ /** * @file stimulation_block.cc * * @author Guillaume Anciaux * @author Till Junge * * @date Tue Jul 29 17:03:26 2014 * * @brief blocks nodes of a fem mesh, thus eliminating equations * * @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_block.hh" #include "container_mesh.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_point_data.hh" __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ StimulationBlock::StimulationBlock(const std::string &name) : LMObject(name) { block = Vector<3, bool>::Zero(); } /* -------------------------------------------------------------------------- */ StimulationBlock::~StimulationBlock() {} /* -------------------------------------------------------------------------- */ template void StimulationBlock::stimulate(_Input &input) { if (input.size() == 0) return; constexpr UInt Dim = _Input::Dim; std::vector directions; for (UInt i = 0; i < Dim; ++i) { if (this->block[i]) directions.push_back(i); } typedef std::vector::iterator dir_iterator; DUMP("Stimulation BLOCK " << " Dim = " << Dim << " block_x = " << this->block[0] << " block_y = " << this->block[1] << " block_z = " << this->block[2], DBG_INFO); for (__attribute__((unused)) auto &&nd : input) { for (dir_iterator dir = directions.begin(); dir != directions.end(); ++dir) { LM_TOIMPLEMENT; // nd.boundary(*dir) = true; } } } /* -------------------------------------------------------------------------- */ /* LMDESC BLOCK blocks nodes of a fem mesh, thus eliminating equations */ /* LMEXAMPLE STIMULATION boundary BLOCK INPUT femdom DIR 1 0 1 */ /* LMHERITANCE action_interface */ void StimulationBlock::declareParams() { StimulationInterface::declareParams(); /* LMKEYWORD DIR Whether to block in x direction*/ this->parseVectorKeyword("DIR", spatial_dimension, this->block, VEC_DEFAULTS(false, false, false)); } /* -------------------------------------------------------------------------- */ -DECLARE_STIMULATION_MAKE_CALL(StimulationBlock) +DECLARE_STIMULATION_MAKE_CALL(StimulationBlock, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_dislo.cc b/src/stimulation/stimulation_dislo.cc index d79c7ac..a72c20a 100644 --- a/src/stimulation/stimulation_dislo.cc +++ b/src/stimulation/stimulation_dislo.cc @@ -1,723 +1,723 @@ /** * @file stimulation_dislo.cc * * @author Guillaume Anciaux * @author Till Junge * @author Jaehyun Cho * * @date Wed Jul 09 21:59:47 2014 * * @brief This stmulation command imposes the Volterra displacement fields of * strainght screw edge or mixed dislocations * * @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_dislo.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_point_data.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ inline Real magnitude(int *vec, UInt nb) { Real mag = 0.0; for (UInt i = 0; i < nb; ++i) { mag += Real(vec[i]) * Real(vec[i]); } return sqrt(mag); } inline Real round(Real value) { Real int_value = (Real)floor(value); if (value - int_value >= 0.5) { return (Real)ceil(value); } else { return (Real)floor(value); } } /* -------------------------------------------------------------------------- */ StimulationDislo::StimulationDislo(const std::string &name) : LMObject(name) { nb_replica = 0; this->dir_replica_char = 'X'; this->burg_edge = this->rotation_matrix; this->normal = this->burg_edge + 3; this->line_dir = this->burg_edge + 6; for (UInt i = 0; i < 3; ++i) { this->pos[i] = 0.0; this->burg_edge_quant[i] = 0; this->normal_quant[i] = 0; this->line_dir_quant[i] = 0; this->dipole[i] = 0; } this->normal_quant[1] = 1; this->line_dir_quant[2] = 1; this->theta = 90; this->fix_border_tol = 0; this->linear = false; this->aniso = false; this->iso = false; for (UInt i = 0; i < 13; ++i) this->elastic_coefs[i] = 0.0; this->c11bar = 0.0; this->lambda = 0.0; this->phi = 0.0; } /* -------------------------------------------------------------------------- */ StimulationDislo::~StimulationDislo() {} /* -------------------------------------------------------------------------- */ void StimulationDislo::init() { // first step: normalize all the vectors to build a coordinate system Real mag_n = magnitude(this->normal_quant, 3); Real mag_l = magnitude(this->line_dir_quant, 3); this->burg_edge_quant[0] = this->normal_quant[1] * this->line_dir_quant[2] - this->normal_quant[2] * this->line_dir_quant[1]; this->burg_edge_quant[1] = this->normal_quant[2] * this->line_dir_quant[0] - this->normal_quant[0] * this->line_dir_quant[2]; this->burg_edge_quant[2] = this->normal_quant[0] * this->line_dir_quant[1] - this->normal_quant[1] * this->line_dir_quant[0]; for (UInt i = 0; i < 3; ++i) { this->normal[i] = Real(this->normal_quant[i]) / mag_n; this->line_dir[i] = Real(this->line_dir_quant[i]) / mag_l; } this->burg_edge[0] = this->normal[1] * this->line_dir[2] - this->normal[2] * this->line_dir[1]; this->burg_edge[1] = this->normal[2] * this->line_dir[0] - this->normal[0] * this->line_dir[2]; this->burg_edge[2] = this->normal[0] * this->line_dir[1] - this->normal[1] * this->line_dir[0]; if (this->dir_replica_char == 'X' || this->dir_replica_char == 'x') { dir_replica = 0; } else if (this->dir_replica_char == 'Y' || this->dir_replica_char == 'y') { dir_replica = 1; } else if (this->dir_replica_char == 'Z' || this->dir_replica_char == 'z') { dir_replica = 2; } else { LM_FATAL("replica direction can only be X, Y or Z"); } // security check regarding anisotropy if ((!this->aniso) && (!this->iso)) { for (UInt i = 0; i < 13; ++i) { if (this->elastic_coefs[i] != 0.0) { LM_FATAL("In case of (an)isotropy, elastic coeffcients shouldn't be " "defined."); } } } else if ((this->aniso) && (this->iso)) { LM_FATAL("Both anisotropy and isotropy cannot be defined"); } else { Real sum_elastic_coefs = 0.0; for (UInt i = 0; i < 13; ++i) { sum_elastic_coefs += this->elastic_coefs[i]; } if (sum_elastic_coefs == 0.0) { LM_FATAL( "In case of (an)isotropy, elastic coeffcients should be defined."); } if (this->linear) LM_FATAL("In case of (an)isotropy, linear field cannot be applied"); Real c11 = this->elastic_coefs[0]; Real c12 = this->elastic_coefs[1]; Real c22 = this->elastic_coefs[4]; Real c66 = this->elastic_coefs[12]; this->c11bar = sqrt(c11 * c22); if (this->aniso) { Real lambda2 = sqrt(c11 / c22); this->lambda = sqrt(lambda2); this->phi = 0.5 * acos((c12 * c12 + 2.0 * c12 * c66 - c11bar * c11bar) / (2.0 * c11bar * c66)); } else { this->lambda = 1.0; this->phi = M_PI / 2.0; } } } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::eval_q_coef(Vector<_Input::Dim> &X) { Real inside = X[0] * X[0] + 2.0 * X[0] * X[1] * this->lambda * cos(this->phi) + X[1] * X[1] * this->lambda * this->lambda; return sqrt(inside); } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::eval_t_coef(Vector<_Input::Dim> &X) { Real inside = X[0] * X[0] - 2.0 * X[0] * X[1] * this->lambda * cos(this->phi) + X[1] * X[1] * this->lambda * this->lambda; return sqrt(inside); } /* -------------------------------------------------------------------------- */ inline Real StimulationDislo::evalArctan(Real X, Real Y) { if ((X >= 0.0) && (Y > 0.0)) { return atan(Y / X); } else if ((X < 0.0) && (Y >= 0.0)) { return atan(-Y / -X) + M_PI; } else if ((X < 0.0) && (Y < 0.0)) { return atan(-Y / -X) - M_PI; } else if ((X >= 0.0) && (Y < 0.0)) { return atan(Y / X); } else { DUMP("X:" << X << ", Y: " << Y, DBG_MESSAGE); DUMP("c11bar:" << this->c11bar << ", lambda: " << lambda << ", phi: " << phi, DBG_MESSAGE); LM_FATAL("unexpected case"); } } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeBeDisplacement(Real b_edge, Vector<_Input::Dim> &X) { Real x_2 = X[0] * X[0]; Real y_2 = X[1] * X[1]; Real epsilon = 1e-15; Real res = b_edge / (2.0 * M_PI) * (atan2(X[1], X[0]) + X[0] * X[1] / (2.0 * (1 - nu) * (x_2 + y_2 + epsilon))); return res; } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeBeDisplacementAniso(Real b_edge, Vector<_Input::Dim> &X) { Real first = 0.0; first += evalArctan((X[0] + X[1] * this->lambda * cos(this->phi)), (X[1] * this->lambda * sin(this->phi))); first += evalArctan((X[0] - X[1] * this->lambda * cos(this->phi)), (X[1] * this->lambda * sin(this->phi))); Real q = eval_q_coef(X); Real t = eval_t_coef(X); Real c12 = this->elastic_coefs[1]; Real c66 = this->elastic_coefs[12]; Real logvalue = log(q / t); Real second = (this->c11bar * this->c11bar - c12 * c12) / (2.0 * this->c11bar * c66 * sin(2.0 * this->phi)) * logvalue; // Real second = 0.0; Real res = -1.0 * b_edge / (4.0 * M_PI) * (-1.0 * first + second); return res; } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeBeDisplacementIso(Real b_edge, Vector<_Input::Dim> &X) { Real first = 0.0; Real second = 0.0; Real res = 0.0; Real c12 = this->elastic_coefs[1]; Real c66 = this->elastic_coefs[12]; first += evalArctan(X[0], X[1]); first += evalArctan(X[0], X[1]); second = (this->c11bar * this->c11bar - c12 * c12) / (4.0 * this->c11bar * c66) * ((2.0 * X[0] * X[1]) / (X[0] * X[0] + X[1] * X[1])); res = b_edge / (4.0 * M_PI) * (first - second); return res; } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeNDisplacement(Real b_edge, Vector<_Input::Dim> &X) { Real x_2 = X[0] * X[0]; Real y_2 = X[1] * X[1]; // Real b_edge2 = b_edge*b_edge; Real epsilon = 1e-15; // return -b_edge/(2.0*M_PI)*((1-2.0*nu)/(4.0*(1-nu))*(log((x_2 + // y_2+epsilon)/b_edge2)) // + (x_2 - y_2)/ // (4.0*(1-nu)*(x_2 + y_2+epsilon)) // ); return -b_edge / (2.0 * M_PI) * ((1 - 2.0 * nu) / (4.0 * (1 - nu)) * log(x_2 + y_2 + epsilon) + (x_2 - y_2) / (4.0 * (1 - nu) * (x_2 + y_2 + epsilon))); } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeNDisplacementAniso(Real b_edge, Vector<_Input::Dim> &X) { Real q = eval_q_coef(X); Real t = eval_t_coef(X); Real c12 = this->elastic_coefs[1]; Real first = (this->c11bar - c12) * cos(this->phi) * log(q * t) / sin(2.0 * this->phi); Real second = (this->c11bar + c12) * sin(this->phi) * atan2((X[1] * X[1] * this->lambda * this->lambda * sin(2.0 * this->phi)), (X[0] * X[0] - this->lambda * this->lambda * X[1] * X[1] * cos(2.0 * this->phi))) / sin(2.0 * this->phi); Real res = this->lambda * b_edge / 4.0 / M_PI / this->c11bar * (first + second); return res; } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeNDisplacementIso(Real b_edge, Vector<_Input::Dim> &X) { Real first = 0.0; Real second = 0.0; Real res = 0.0; Real c12 = this->elastic_coefs[1]; first = (this->c11bar - c12) / 2. * log(X[0] * X[0] + X[1] * X[1]); second = (this->c11bar + c12) * ((X[1] * X[1]) / (X[0] * X[0] + X[1] * X[1])); res = b_edge / 4.0 / M_PI * (first + second); return res; } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeLDisplacement(Real b_screw, Vector<_Input::Dim> &X) { Real res = b_screw / (2.0 * M_PI) * (atan2(X[1], X[0])); return res; } /* -------------------------------------------------------------------------- */ template inline Real StimulationDislo::computeLDisplacementAniso(Real b_screw, Vector<_Input::Dim> &X) { Real c44 = this->elastic_coefs[9]; Real c45 = this->elastic_coefs[10]; Real c55 = this->elastic_coefs[11]; Real res = b_screw / 2.0 / M_PI * atan2(sqrt(c44 * c55 - c45 * c45) * X[1], c44 * X[0] - c45 * X[1]); return res; } /* -------------------------------------------------------------------------- */ inline Real StimulationDislo::computeLinearDisplacement(Real burg, Real *X, Real xlo, Real xhi) { Real res = burg / 2. / (xhi - xlo) * (X[0] - xhi); if (X[1] > 0) return -1.0 * res; else return 1.0 * res; } /* -------------------------------------------------------------------------- */ template void StimulationDislo::fix_boundary(typename _Input::Ref &atom, const Vector<_Input::Dim> &xmin, const Vector<_Input::Dim> &xmax) { Real upper_diff = atom.position()[this->dir_replica] - xmax[this->dir_replica]; Real lower_diff = atom.position()[this->dir_replica] - xmin[this->dir_replica]; if (fabs(upper_diff) < this->fix_border_tol) { Vector<_Input::Dim> pos = atom.position(); pos[this->dir_replica] = xmax[this->dir_replica]; atom.position() = pos; } else if (fabs(lower_diff) < this->fix_border_tol) { Vector<_Input::Dim> pos = atom.position(); pos[this->dir_replica] = xmin[this->dir_replica]; atom.position() = pos; } } /* -------------------------------------------------------------------------- */ template void StimulationDislo::stimulate(Cont &cont) { this->stimulate_to_clean(cont, current_stage); LM_TOIMPLEMENT; // if (this->fix_border_tol != 0.) { // for (auto &&at : cont) { // this->fix_boundary(at, xmin, xmax); // } // } } /* -------------------------------------------------------------------------- */ template void StimulationDislo::stimulate_to_clean(Cont &cont, UInt) { constexpr UInt Dim = Cont::Dim; Cube &cube = cont.getBoundingBox(); __attribute__((unused)) auto &&xmax = cube.getXmax(); __attribute__((unused)) auto &&xmin = cube.getXmin(); auto domainSize = cube.getSize(); Real b_edge; Real b_screw; b_edge = burgers * sin(theta * M_PI / 180.0); b_screw = -burgers * cos(theta * M_PI / 180.0); // identify the direction of propagation (Be) Real BeDir0 = normal_quant[1] * line_dir_quant[2] - normal_quant[2] * line_dir_quant[1]; Real BeDir1 = -(normal_quant[0] * line_dir_quant[2] - normal_quant[2] * line_dir_quant[0]); Real BeDir2 = normal_quant[0] * line_dir_quant[1] - normal_quant[1] * line_dir_quant[0]; Real BeMag = sqrt(BeDir0 * BeDir0 + BeDir1 * BeDir1 + BeDir2 * BeDir2); BeDir0 /= BeMag; BeDir1 /= BeMag; BeDir2 /= BeMag; // UInt bedir = 4; // if (fabs(abs(BeDir0) - 1.0) < 1e-4) // bedir = 0; // else if (fabs(abs(BeDir1) - 1.0) < 1e-4) // bedir = 1; // else if (fabs(abs(BeDir2) - 1.0) < 1e-4) // bedir = 2; // else // LM_FATAL("cannot identify propagation direction"); throw; // this was temporarily commented // Real xloBe = xmin[bedir]; // Real xhiBe = xmax[bedir]; LM_ASSERT(Dim == 3, "this stimulator works only in 3D"); Cube c = cont.getBoundingBox(); double dipolelength = 0; for (UInt i = 0; i < 3; i++) { if (this->dipole[i] == 1) { dipolelength = c.getSize(i); } } for (auto &&at : cont) { auto X0 = at.position0(); // position in physical space auto X = at.position(); // position in physical space Vector x = Vector::Zero(); // position in the coordinates of the dislo Vector Disp = Vector::Zero(); Real disp[3]; for (int rep = -1 * nb_replica; rep < nb_replica + 1; ++rep) { X = Vector(X0); for (UInt i = 0; i < 3; i++) { if (this->dipole[i] == 1) { if (X[i] < 0) { X[i] = -1.0 * (X[i] + (dipolelength / 2.0) * 0.5); } else { X[i] = 1.0 * (X[i] - (dipolelength / 2.0) * 0.5); } } } for (UInt i = 0; i < 3; ++i) { X[i] -= pos[i]; disp[i] = 0.0; } // calculate coordinates of images if (this->linear == false) X[dir_replica] += domainSize[dir_replica] * Real(rep); // rotate point for (UInt i = 0; i < 3; ++i) { x[i] = 0.0; for (UInt j = 0; j < 3; ++j) { x[i] += this->rotation_matrix[3 * i + j] * X[j]; } } if (this->linear == false) { if (this->aniso) { // compute anisotropy volterra displacement LM_TOIMPLEMENT; // disp[0] = computeBeDisplacementAniso(b_edge, x); if (rep == 0) { LM_TOIMPLEMENT; // disp[1] = computeNDisplacementAniso(b_edge, x); } LM_TOIMPLEMENT; // disp[2] = computeLDisplacementAniso(b_screw, x); } else if (this->iso) { // compute isotropy volterra displacement LM_TOIMPLEMENT; // disp[0] = computeBeDisplacementIso(b_edge, x); if (rep == 0) { LM_TOIMPLEMENT; // disp[1] = computeNDisplacementIso(b_edge, x); } LM_TOIMPLEMENT; // disp[2] = computeLDisplacementAniso(b_screw, x); } else { // compute isotropy volterra displacement LM_TOIMPLEMENT; // disp[0] = computeBeDisplacement(b_edge, x); if (rep == 0) { LM_TOIMPLEMENT; // disp[1] = computeNDisplacement(b_edge, x); } LM_TOIMPLEMENT; // disp[2] = computeLDisplacement(b_screw, x); } } else { // compute linear displacements throw; // disp[0] = computeLinearDisplacement(b_edge, x, xloBe, xhiBe); // disp[2] = computeLinearDisplacement(b_screw, x, xloBe, xhiBe); } // compensate accumulated images in x-axis if ((rep == 0) and (this->linear == false)) { if (x[1] > 0) { disp[0] -= Real(nb_replica) * b_edge / 2; disp[2] -= Real(nb_replica) * b_screw / 2; } else { disp[0] += Real(nb_replica) * b_edge / 2; disp[2] += Real(nb_replica) * b_screw / 2; } } // rotate back for (UInt i = 0; i < 3; ++i) { for (UInt j = 0; j < 3; ++j) { Disp[i] += this->rotation_matrix[3 * j + i] * disp[j]; } } if (this->linear) break; // exit replica for loops } at.position() += Disp; if ((fabs(at.position0()[0] + 26.3245) < 1e-1) && (fabs(at.position0()[1] - 102.843) < 1e-1) && (fabs(at.position0()[2] + 22.8263) < 1e-1)) { DUMP("HERE", DBG_MESSAGE); } } } /* -------------------------------------------------------------------------- */ /* LMDESC DISLO This stmulation command imposes the Volterra displacement\\ fields of straight screw, edge, or combination of two dislocations\\\\ Displacement field of Screw dislocation: \begin{equation} u_{z}(x, y) = \frac{b}{2\pi} \cdot atan2(y,x) \end{equation} Displacement fields of Edge dislocation: \begin{equation} u_{x}(x, y) = \frac{b}{2\pi} \cdot \left( atan2(y,x) + \frac{xy}{2 (1-\nu) (x^2 + y^2)} \right) \end{equation} \begin{equation} u_{y}(x, y) = -\frac{b}{2\pi} \cdot \left( \frac{1-2\nu}{4(1 - \nu)}\ln(x^2 + y^2) + \frac{x^2 - y^2}{4(1-\nu)(x^2 + y^2)} \right) \end{equation} are introduced. */ /* LMEXAMPLE STIMULATION myDislo DISLO INPUT md STAGE PRE_DUMP BURGERS burgers * POS 0 0 0 REPLICA 10 0 POISSON 0.3468 THETA 180 ONESHOT 0 */ /* LMHERITANCE action_interface */ void StimulationDislo::declareParams() { StimulationInterface::declareParams(); /* LMKEYWORD BURGERS Specify the magnitude of Burger's vector */ this->parseKeyword("BURGERS", this->burgers); /* LMKEYWORD POS Specify the coordinates of a point on the dislocation line */ this->parseVectorKeyword("POS", 3, this->pos); /* LMKEYWORD SLIP_PLANE Specifies the normal vector to the slip plane */ this->parseVectorKeyword("SLIP_PLANE", 3, this->normal_quant); /* LMKEYWORD LINE_DIR Specifies the normal vector to the slip plane */ this->parseVectorKeyword("LINE_DIR", 3, this->line_dir_quant); /* LMKEYWORD POISSON Specify Poisson's ratio of the material */ this->parseKeyword("POISSON", this->nu); // Resolution and width are temporarily taken out, because they don't do // anything, // plus, it's not clear how to spread a screw dislo, since it doesn't have a // slip plane ///* LM KEYW ORD RESOLUTION // Specify the number of increments in the spread of dislocation density // in // Z direction //*/ // parseKey word("RESOLUTION",resolution); // ///* LM KEY esWORD WIDTH // Specify the width of the core dislocation in Z direction //*/ // parseKe yword("WIDTH",width); /* LMKEYWORD NB_REPLICA Specify the number of replicas in burgers edge components */ this->parseKeyword("NB_REPLICA", this->nb_replica); /* LMKEYWORD DIR_REPLICA Specify the direction of replication (X, Y or Z) */ this->parseKeyword("DIR_REPLICA", this->dir_replica_char, 'X'); /* LMKEYWORD THETA Specify the angle of burgers vector with dislocation line in unit of degree\\ For example, 90 and 270 impose pure edge burgers vector,\\ 0 and 180 impose pure screw burgers vector, and\\ other angles impose combination of two kinds of burgers vector */ this->parseKeyword("THETA", this->theta); /* LMKEYWORD FIX_BORDER_TOL specify a tolerance whithin which border points are forces to be exacly on the border: ------------- ------------- | | | | | | | | / | => | | \ | => | | | | | | | | | | ------------- ------------- */ this->parseKeyword("FIX_BORDER_TOL", this->fix_border_tol, 0.); /* LMKEYWORD DIPOLE Specify the direction where the dislocation dipole creation in X, Y or Z. As results, two dislocation are inserted with full periodicity in all directions. %% For example, DIPOLE 0 1 0 creates: %% ------------- ------------- %% | | | | %% | | | T | %% | _|_ | => | | %% | | | _|_ | %% | | | | %% ------------- ------------- %% where T and _|_ are same burger vectors with opposite signs. */ this->parseVectorKeyword("DIPOLE", 3, this->dipole, VEC_DEFAULTS(0, 0, 0)); /* LMKEYWORD LINEAR Specify the keyword to employ the linear displacement fields instead of volterra solution */ this->parseTag("LINEAR", this->linear, false); /* LMKEYWORD ANISO Boolean to set anisotropicity: [1] Elastic Strain Fields and Dislocation Mobility, Volume 31, J.Lothe Dislocations in anisotropic media [2] Plane-strain discrete dislocation plasticity incorporating anisotropic elasticity IJSS 48(2), January 2011 */ this->parseTag("ANISO", this->aniso, false); /* L2MKEYWORD ELASTIC_COEFS Specify the thirteen elastic coefficients of the constitutive matrix of anisotropic media in the following order: C11, C12, C13, C21, C22, C23, C31, C32, C33, C44, C45, C55, C66 */ // this->parseVectorKeyword("ELASTIC_COEFS", 13, this->elastic_coefs); /* LMKEYWORD ISO Boolean to set isotropicity : [1] Elastic Strain Fields and Dislocation Mobility, Volume 31, J.Lothe Dislocations in anisotropic media [2] Plane-strain discrete dislocation plasticity incorporating anisotropic elasticity IJSS 48(2), January 2011 */ this->parseTag("ISO", this->iso, false); } /* -------------------------------------------------------------------------- */ -DECLARE_STIMULATION_MAKE_CALL(StimulationDislo) +DECLARE_STIMULATION_MAKE_CALL(StimulationDislo, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_impulse.cc b/src/stimulation/stimulation_impulse.cc index 839a61a..63b95ee 100644 --- a/src/stimulation/stimulation_impulse.cc +++ b/src/stimulation/stimulation_impulse.cc @@ -1,97 +1,98 @@ /** * @file stimulation_impulse.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This stimulator sets an initial impulse displacement field generating * a wave * * @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_impulse.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_point_data.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ StimulationImpulse::StimulationImpulse(const std::string &name) : LMObject(name), impulse_compute("computeImpulse:" + name) {} /* -------------------------------------------------------------------------- */ StimulationImpulse::~StimulationImpulse() {} /* -------------------------------------------------------------------------- */ void StimulationImpulse::init() { impulse_compute.init(); } /* -------------------------------------------------------------------------- */ template void StimulationImpulse::stimulate(Cont &cont) { constexpr UInt Dim = Cont::Dim; impulse_compute.compute(cont); auto &imposed_field = impulse_compute.evalArrayOutput(); if (imposed_field.cols() == Dim) { for (auto &&[at, field] : zip(cont, imposed_field.rowwise())) { at.displacement() = Vector(field); } } else { const int direction = impulse_compute.direction; LM_ASSERT(cont.size() == imposed_field.size(), "The impossible happened " << cont.size() << " != " << imposed_field.size()); for (auto &&[at, field] : zip(cont, imposed_field)) { at.displacement()[direction] = field; } } } /* -------------------------------------------------------------------------- */ /* LMDESC IMPULSE This stimulator sets an initial impulse displacement field that will generate - a wave. It uses the IMPULSE compute to generate a displacement field. + a wave. It uses the IMPULSE compute to generate a displacement field. + See keywords/filter/impulse for more information on the generated wave. */ /* LMHERITANCE compute_impulse */ /* LMEXAMPLE STIMULATION impulse IMPULSE INPUT central_md LWAVE WaveLength*r0 DIRECTION 0 INTENSITY 1.1e-3 ONESHOT 0 */ void StimulationImpulse::declareParams() { StimulationInterface::declareParams(); this->addSubParsableObject(impulse_compute); } /* -------------------------------------------------------------------------- */ DECLARE_STIMULATION_MAKE_CALL(StimulationImpulse) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_langevin.cc b/src/stimulation/stimulation_langevin.cc index 67b05e6..63e12b8 100644 --- a/src/stimulation/stimulation_langevin.cc +++ b/src/stimulation/stimulation_langevin.cc @@ -1,395 +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: .. math:: dW^n = \int_{x^n}^{x^{n+1}} f \cdot dx which for the case of the velocity verlet scheme is equivalent to .. 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" and 5 entries is registered. These entries are: - 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) +DECLARE_STIMULATION_MAKE_CALL(StimulationLangevin, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_rigidify.cc b/src/stimulation/stimulation_rigidify.cc index a228592..6f8a1c3 100644 --- a/src/stimulation/stimulation_rigidify.cc +++ b/src/stimulation/stimulation_rigidify.cc @@ -1,130 +1,130 @@ /** * @file stimulation_rigidify.cc * * @author Guillaume Anciaux * * @date Wed Jul 09 21:59:47 2014 * * @brief This stimulator allows to reduce the force to let a group act like a * rigid body * * @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_rigidify.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" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ StimulationRigidify::StimulationRigidify(const std::string &name) : LMObject(name) { for (UInt i = 0; i < 3; ++i) { total_force[i] = 0.; } factor = 1.; } /* -------------------------------------------------------------------------- */ StimulationRigidify::~StimulationRigidify() {} /* -------------------------------------------------------------------------- */ template void StimulationRigidify::stimulate(Cont &cont) { constexpr UInt Dim = Cont::Dim; UInt cpt = cont.size(); Vector f = Vector::Zero(); for (auto &&at : cont) { f += at.force(); } // then reduce for multiple processors CommGroup &group = cont.getCommGroup(); group.allReduce(f.data(), Dim, "reduce force", OP_SUM); group.allReduce(&cpt, 1, "reduce number atoms", OP_SUM); f += total_force.cast().block(0, 0); f /= cpt; for (auto &&at : cont) { at.force() = f / factor; } } /* -------------------------------------------------------------------------- */ /* LMDESC RIGIDIFY This stimulator allows to reduce the force acting over a set of degrees of freedom to a single scalar and assign this last to all forces of that set. Such a set will then act like a rigid body. More precisely the algorithm is:\\ \begin{equation} f_{global} = \sum_i f_i + f^{ext} \end{equation} \begin{equation} f_i := \frac{f_{global}}{\#dof \cdot factor} \end{equation} where $\#dof$ is the number of dof in the subset $f^{ext}$ (GLOBAL\_FORCE) is an external force applied onto the rigid body and $factor$ (FACTOR) is a flexibility factor to accelerate linear motion. */ /* LMHERITATE action_interface */ /* LMEXAMPLE STIMULATION rigid RIGIDIFY INPUT md GLOBAL_FORCE 1 0 0 MASS_FACTOR * 1e-3 */ //! most generic stimulate function void StimulationRigidify::declareParams() { StimulationInterface::declareParams(); this->changeDefault("STAGE", PRE_STEP3); /* LMKEYWORD TOTAL_FORCE Specifiy the external force to be applied to the rigid body. */ this->parseVectorKeyword("TOTAL_FORCE", spatial_dimension, total_force); /* LMKEYWORD FACTOR Correction factor to augment/reduce the forces applied to the rigid body. */ this->parseKeyword("FACTOR", factor, 1.); } /* -------------------------------------------------------------------------- */ -DECLARE_STIMULATION_MAKE_CALL(StimulationRigidify) +DECLARE_STIMULATION_MAKE_CALL(StimulationRigidify, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_temperature.cc b/src/stimulation/stimulation_temperature.cc index 59d01f9..0c0fe24 100644 --- a/src/stimulation/stimulation_temperature.cc +++ b/src/stimulation/stimulation_temperature.cc @@ -1,142 +1,142 @@ /** * @file stimulation_temperature.cc * * @author Guillaume Anciaux * @author Srinivasa Babu Ramisetti * * @date Wed Jul 09 21:59:47 2014 * * @brief Set a given temperature to a set of DOFs * * @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_temperature.hh" #include "compute_ekin.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ /* LMDESC TEMPERATURE This stimulator sets initial velocity to input a temperature. */ /* LMEXAMPLE STIMULATION TEMPERATURE TEMP 100 SEED 32 */ /* LMEXAMPLE STIMULATION TEMPERATURE FLUX_FLAG FLUX 1 SEED 32 */ /* LMHERITANCE action_interface */ inline void StimulationTemperature::declareParams() { StimulationInterface::declareParams(); /* 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 HEAT_RATE Set the heat rate per node (only for finite elements) */ this->parseKeyword("HEAT_RATE", heat_rate_per_node, 0.); /* LMKEYWORD FLUX_FLAG Set the flux per node (only for finite elements) */ this->parseTag("FLUX_FLAG", flux_flag, false); } /* -------------------------------------------------------------------------- */ StimulationTemperature::StimulationTemperature(const std::string &name) : LMObject(name), flux_flag(false) {} /* -------------------------------------------------------------------------- */ StimulationTemperature::~StimulationTemperature() {} /* -------------------------------------------------------------------------- */ template enable_if_md StimulationTemperature::stimulate(Cont &cont) { constexpr UInt Dim = Cont::Dim; for (auto &&at : cont) { at.velocity() = Vector::Random(); } // remove rigid body motion Vector average_v = Vector::Zero(); for (auto &&at : cont) { average_v += at.velocity(); } average_v /= cont.size(); for (auto &&at : cont) { at.velocity() -= average_v; } ComputeEKin compute_ekin("ComputeEkin:" + this->getID()); compute_ekin.compute("domain"_input = cont); Real T_ = compute_ekin.evalArrayOutput("Temperature").get(0); Real alpha = sqrt(temperature / T_); for (auto &&at : cont) { at.velocity() *= alpha; } compute_ekin.build(cont); T_ = compute_ekin.evalArrayOutput("Temperature").get(0); DUMP("temperature measured after: " << T_, DBG_MESSAGE); } /* -------------------------------------------------------------------------- */ template enable_if_mesh StimulationTemperature::stimulate(Cont & // cont ) { LM_TOIMPLEMENT; // if (!flux_flag) { // for (auto &&n : cont) { // n.temperature() = this->temperature; // } // } else { // for (auto &&n : cont) { // n.externalHeatRate() = this->heat_rate_per_node; // } // } } /* -------------------------------------------------------------------------- */ -DECLARE_STIMULATION_MAKE_CALL(StimulationTemperature) +DECLARE_STIMULATION_MAKE_CALL(StimulationTemperature, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_veluni.cc b/src/stimulation/stimulation_veluni.cc index 6025c3f..7102804 100644 --- a/src/stimulation/stimulation_veluni.cc +++ b/src/stimulation/stimulation_veluni.cc @@ -1,141 +1,141 @@ /** * @file stimulation_veluni.cc * * @author Guillaume Anciaux * @author Till Junge * * @date Wed Jul 09 21:59:47 2014 * * @brief This stimulator allows to enforce a constant velocity motion * * @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_veluni.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "ref_point_data.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template void apply_vel(Cont &cont, Real vel) { for (auto &&at : cont) { DUMP("proceed stimulation of atom " << at, DBG_DETAIL); at.velocity()[Dir] = vel; at.force()[Dir] = 0.0; } } /* -------------------------------------------------------------------------- */ template void apply_vel(Cont &cont, Vector vec) { for (auto &&at : cont) { DUMP("proceed stimulation of atom " << at, DBG_DETAIL); for (UInt i = 0; i < Dim; ++i) { at.velocity()[i] = vec[i]; } at.force().setZero(); } } /* -------------------------------------------------------------------------- */ StimulationVelUni::StimulationVelUni(const std::string &name) : LMObject(name) { for (UInt i = 0; i < 3; ++i) vec[i] = 0.; } /* -------------------------------------------------------------------------- */ -StimulationVelUni::~StimulationVelUni(){} +StimulationVelUni::~StimulationVelUni() {} /* -------------------------------------------------------------------------- */ template void StimulationVelUni::stimulate(Cont &cont) { constexpr UInt Dim = Cont::Dim; if (this->directions.size() == 0) apply_vel(cont, this->vec.block(0, 0).template cast()); else { for (UInt i = 0; i < directions.size(); ++i) { int dir = directions[i]; if (dir >= 0 && dir < int(Dim)) { switch (dir) { case 0: apply_vel(cont, this->vec[0]); break; case 1: apply_vel(cont, this->vec[1]); break; case 2: apply_vel(cont, this->vec[2]); break; default: LM_FATAL("unknown dimension"); } } else LM_FATAL("unknown dimension"); } } } /* -------------------------------------------------------------------------- */ /* LMDESC VELUNI This stimulator allows to enforce a constant velocity motion. The way chosen to enforce this is by cancelling the force and imposing a constant velocity field as provided by keyword VEC. \ \ The algorithm is:\ \ \begin{align*} \forall i \in dir, \quad f_i = 0 \quad v_i = VEC \end{align*} */ /* LMEXAMPLE STIMULATION uni VELUNI INPUT md VEC 1e-5 0 0 */ void StimulationVelUni::declareParams() { StimulationInterface::declareParams(); /* LMKEYWORD VEC Specifies the vector to be imposed to the velocity field. */ this->parseVectorKeyword("VEC", spatial_dimension, vec); /* LMKEYWORD DIR Specifies the direction of the constraint 0 is X, 1 is Y, 2 is Z and -1(default) is all directions */ this->parseKeyword("DIR", directions, std::vector{}); } /* -------------------------------------------------------------------------- */ -DECLARE_STIMULATION_MAKE_CALL(StimulationVelUni) +DECLARE_STIMULATION_MAKE_CALL(StimulationVelUni, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_zero.cc b/src/stimulation/stimulation_zero.cc index 042ca4b..e62979a 100644 --- a/src/stimulation/stimulation_zero.cc +++ b/src/stimulation/stimulation_zero.cc @@ -1,150 +1,206 @@ + + /** * @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) { +template +enable_if_md 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) { + 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]; + } + } else if (f == _torque) { + for (iter i = directions.begin(); i != directions.end(); ++i) { + at.torque()[*i] += (COEF - 1) * at.torque()[*i]; + } + } else if (f == _angular_velocity) { + for (iter i = directions.begin(); i != directions.end(); ++i) { + at.angular_velocity()[*i] += (COEF - 1) * at.angular_velocity()[*i]; + } + } + } + } +} + +template +enable_if_not_md 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) { 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: - displacement : selects the displacement field - velocity : selects the velocity field - force : selects the force field + - torque : selects the torque field (valid for DEM) + - angular_velocity : selects the angular velocity field (valid for DEM) - 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) +DECLARE_STIMULATION_MAKE_CALL(StimulationZero, ) __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_zero.hh b/src/stimulation/stimulation_zero.hh index e0de7a4..95b141e 100644 --- a/src/stimulation/stimulation_zero.hh +++ b/src/stimulation/stimulation_zero.hh @@ -1,60 +1,62 @@ /** * @file stimulation_zero.hh * * @author Guillaume Anciaux * @author Till Junge * * @date Mon Nov 25 15:05:56 2013 * * @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 . * */ #ifndef STIMULATION_ZERO_H #define STIMULATION_ZERO_H /* -------------------------------------------------------------------------- */ #include "stimulation_interface.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ class StimulationZero : public StimulationInterface { public: DECLARE_STIMULATION(StimulationZero, _or); //! most generic stimulate function - template void stimulate(_Input &cont); + template enable_if_md stimulate(Cont &cont); + template enable_if_not_md stimulate(Cont &cont); + private: std::vector field; Real COEF; bool average_flag; //* direction in which to apply the stimulation int direction; }; /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif