diff --git a/src/continuum/akantu/akantu_dof_wrapper.hh b/src/continuum/akantu/akantu_dof_wrapper.hh index f5cd19c..b95ad3e 100644 --- a/src/continuum/akantu/akantu_dof_wrapper.hh +++ b/src/continuum/akantu/akantu_dof_wrapper.hh @@ -1,139 +1,139 @@ /** * @file akantu_dof_wrapper.hh * * @author Guillaume Anciaux * @author Till Junge * @author Srinivasa Babu Ramisetti * * @date Mon Jun 16 17:13:26 2014 * * @brief This is a helper class to manage Akantu models * * @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 __LIBMULTISCALE_AKANTU_DOF_WRAPPER_HH__ #define __LIBMULTISCALE_AKANTU_DOF_WRAPPER_HH__ /* -------------------------------------------------------------------------- */ #include "lm_common.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { #ifdef AKANTU_HEAT_TRANSFER class HeatTransferModel; #endif } // namespace akantu /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ template class AkantuDOFWrapper { public: using MyFEM = typename akantu::SolidMechanicsModel::MyFEEngineType; AkantuDOFWrapper(); ~AkantuDOFWrapper(); //! get solidmechanics model object akantu::SolidMechanicsModel &getSolidMechanicsModel(); const akantu::SolidMechanicsModel &getSolidMechanicsModel() const; //! set solidmechanics model object virtual void setSolidMechanicsModel(akantu::SolidMechanicsModel &model, akantu::Array &force); //! set force vector // virtual void setForce(akantu::Array &force); //! get force vector akantu::Array &getForce(); #ifdef AKANTU_HEAT_TRANSFER //! get heattransfer model object akantu::HeatTransferModel &getHeatTransferModel(); //! set heattransfer model object virtual void setHeatTransferModel(akantu::HeatTransferModel &model); #endif //! return the initialization state bool isInitialized() const; //! return the element type const inline akantu::ElementType &getType() { return this->type; } inline UInt getNbNodesPerElem() { return this->nb_nodes_per_element; } inline UInt *getConnectivity() { return this->connectivity; } //! return the id of the material of the element index inline const UInt &getMaterialID(UInt index) { return (*this->material_per_element)(index); } //! return the id in its material of element index inline const UInt &getIDInMaterial(UInt index) { return (*this->material_local_numbering)(index); } inline const UInt &getNbQuadraturePoints() { return this->nb_quadrature_points; } const akantu::Array &getGradU(UInt material_index); const akantu::Array &getStress(UInt material_index); /// function to print the contain of the class virtual void printself(std::ostream &stream) const; private: akantu::SolidMechanicsModel *solid_mechanics_model; #ifdef AKANTU_HEAT_TRANSFER akantu::HeatTransferModel *heat_transfer_model; #endif //! element type in the mesh akantu::ElementType type; //! material id per element - akantu::Array *material_per_element; + const akantu::Array *material_per_element; //! gives the index used in material for each element - akantu::Array *material_local_numbering; + const akantu::Array *material_local_numbering; //! pointers of strain vectors ordered by material std::vector *> grad_u; //! pointers of stress vectors ordered by material std::vector *> stresses; //! pointers of stress vectors ordered by material akantu::Array *force; //! connectivity array UInt *connectivity; //! nb nodes per element UInt nb_nodes_per_element; UInt nb_quadrature_points; }; /* -------------------------------------------------------------------------- */ template inline std::ostream &operator<<(std::ostream &os, AkantuDOFWrapper &q) { q.printself(os); return os; } __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_AKANTU_DOF_WRAPPER_HH__ */ diff --git a/src/continuum/akantu/domain_akantu.cc b/src/continuum/akantu/domain_akantu.cc index ea6248b..467f584 100644 --- a/src/continuum/akantu/domain_akantu.cc +++ b/src/continuum/akantu/domain_akantu.cc @@ -1,842 +1,842 @@ /** * @file domain_akantu.cc * * @author Guillaume Anciaux * @author Till Junge * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * * @date Mon Jul 21 10:32:33 2014 * * @brief This is the model wrapping Akantu * * @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 "domain_akantu.hh" #include "lib_geometry.hh" #include "lm_common.hh" #include "math_tools.hh" #include "units.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template DomainAkantu::~DomainAkantu() {} /* -------------------------------------------------------------------------- */ template DomainAkantu::DomainAkantu(LMID ID, CommGroup &GID) : LMObject(ID), DomainContinuum, ContainerNodesAkantu, Dim>(GID), surface_tag(""), boundary_start(0), boundary_is_traction(false), boundary_has_been_applied(false), total_force(0, Dim) { for (UInt i = 0; i < 3; ++i) { pbc[i] = 0; scale[i] = 0.; shift[i] = 0.; for (UInt j = 0; j < Dim; ++j) { this->boundary_stress[Dim * i + j] = std::numeric_limits::max(); } } } /* -------------------------------------------------------------------------- */ template void DomainAkantu::scaleMesh() { auto &nodes = mesh->getNodes(); for (auto &n : akantu::make_view(nodes, Dim)) { for (UInt i = 0; i < Dim; ++i) n[i] *= scale[i]; } } /* -------------------------------------------------------------------------- */ template void DomainAkantu::shiftMesh() { auto &nodes = mesh->getNodes(); for (auto &n : akantu::make_view(nodes, Dim)) { for (UInt i = 0; i < Dim; ++i) n[i] += shift[i]; } } /* -------------------------------------------------------------------------- */ template void DomainAkantu::initMesh() { - mesh = std::make_unique(Dim, this->getID(), 0); + mesh = std::make_unique(Dim, this->getID()); UInt prank = this->getCommGroup().getMyRank(); if (prank == 0) { if (mesh_filename != "") { mesh->read(mesh_filename); #ifndef LM_OPTIMIZED UInt ntypes = 0; for (auto &type [[gnu::unused]] : mesh->elementTypes(Dim)) { ++ntypes; } LM_ASSERT(ntypes <= 1, "At current time Akantu plugin cannot support to load" << " more than one element type for volume elements"); #endif // LM_OPTIMIZED scaleMesh(); shiftMesh(); // UInt nb_fragments = mesh->createClusters(Dim); // DUMP("the number of fragments is " << nb_fragments,DBG_MESSAGE); } else { LM_ASSERT(Dim == 1, "automatic mesh generation only possible in 1D" << "a mesh file should have been provided using the" << " MESH_FILENAME keyword"); const akantu::ElementType type = akantu::_segment_2; mesh->addConnectivityType(type); Ball &b = dynamic_cast( *GeometryManager::getManager().getGeometry(this->geometries[0])); // number of element on each side of the zero UInt nx = (UInt)(nearbyint((b.rMax() - b.rMin()) / elem_size)); Real center = b.getCenter(0); Real rmin = b.rMin(); Real rmax = b.rMax(); Real range = rmax - rmin; Real e_size = range / nx; DUMP("element size " << e_size, DBG_MESSAGE); UInt nb_nodes; nb_nodes = 2 * nx + 1; if (rmin > 0) ++nb_nodes; UInt nb_elems; nb_elems = 2 * nx; akantu::Array &nodes = akantu::MeshAccessor(*mesh).getNodes(); nodes.resize(nb_nodes); // create nodes from -1 to 0 for (UInt i = 0; i < nx + 1; ++i) { nodes(i) = static_cast(i) / static_cast(nx) - 1.0; nodes(i) = nodes(i) * range - rmin + center; DUMP("creating node at position " << nodes(i), DBG_DETAIL) } // create nodes from 0 to 1 if (rmin == 0) { for (UInt i = 0; i < nx + 1; ++i) { nodes(i + nx) = static_cast(i) / static_cast(nx); nodes(i + nx) = nodes(i + nx) * range + rmin + center; DUMP("creating node at position " << nodes(i + nx), DBG_DETAIL); } } else for (UInt i = 0; i < nx + 1; ++i) { nodes(i + nx + 1) = static_cast(i) / static_cast(nx); nodes(i + nx + 1) = nodes(i + nx + 1) * range + rmin + center; DUMP("creating node at position " << nodes(i + nx + 1), DBG_DETAIL) } akantu::Array &connectivity = akantu::MeshAccessor(*mesh).getConnectivity(type); connectivity.resize(nb_elems); for (UInt i = 0; i < nx; ++i) { connectivity(i, 0) = i; connectivity(i, 1) = i + 1; } if (rmin == 0) for (UInt i = 0; i < nx; ++i) { connectivity(i + nx, 0) = i + nx; connectivity(i + nx, 1) = i + nx + 1; } else for (UInt i = 0; i < nx; ++i) { connectivity(i + nx, 0) = i + nx + 1; connectivity(i + nx, 1) = i + nx + 2; } akantu::MeshAccessor(*mesh).makeReady(); } } } /* ------------------------------------------------------------------------- */ template void DomainAkantu::initMaterial() { if (this->material_filename == "") { LM_ASSERT(Dim == 1, "automatic material property definition only possible in 1D" << "a material file should have been provided using the" << " MATERIAL_FILENAME keyword"); DUMP("Automatic parameter generation for LJ and lattice provided " "parameters", DBG_WARNING); akantu::Material &mat = model->registerNewMaterial("linear_LJ", "elastic", ""); UInt ordre = static_cast(rcut / r0); Real coeff = MathTools::ipow(sigma / r0, 6); Real E = 24.0 * epsilon * coeff / r0 * (-7.0 * MathTools::zeta(6, ordre) + 26.0 * coeff * MathTools::zeta(12, ordre)); mat.setParam("E", E); mat.setParam("nu", 0.0); mat.setParam("rho", Real(atomic_mass / r0)); DUMP("mass density " << (Real)(atomic_mass / r0), DBG_MESSAGE); DUMP("r0 " << r0, DBG_MESSAGE); DUMP("atomic_mass " << atomic_mass, DBG_MESSAGE); this->model->initMaterials(); } } /* ------------------------------------------------------------------------- */ template void DomainAkantu::init() { auto &group = this->getCommGroup(); if (!group.amIinGroup()) return; if (this->material_filename != "") akantu::getStaticParser().parse(this->material_filename); akantu::debug::setDebugLevel(akantu::dbl0); // akantu::debug::setDebugLevel(akantu::dblTrace); initMesh(); DUMP("created " << this->mesh->getNodes().size() << " nodes ", DBG_INFO_STARTUP); aka_comm = std::make_unique(); aka::as_type(aka_comm->getCommunicatorData()) .setMPICommunicator(group.getMPIComm()); this->mesh->distribute(akantu::_communicator = *aka_comm); bool pbc_activated = false; for (UInt i = 0; i < Dim; ++i) { pbc_activated |= pbc[i]; } if (pbc[0]) mesh->makePeriodic(akantu::_x); if (pbc[1]) mesh->makePeriodic(akantu::_y); if (pbc[2]) mesh->makePeriodic(akantu::_z); if (pbc_activated) { // auto &akantu_pairs = this->model->getPBCPairs(); // for (auto &&pair : akantu_pairs) { // this->pbc_pairs.push_back(pair); // } // dans mesh.cc // std::unordered_map periodic_slave_master; // std::unordered_multimap periodic_master_slave; } this->model = std::make_unique(*mesh, Dim, this->getID()); this->model->initFull( akantu::SolidMechanicsModelOptions(akantu::_explicit_lumped_mass)); this->initMaterial(); this->model->setF_M2A(code_unit_system.ft_m2v); this->total_force.resize(this->model->getInternalForce().size()); this->mesh_container.getContainerElems().setSolidMechanicsModel( *this->model, this->total_force); this->mesh_container.getContainerNodes().setSolidMechanicsModel( *this->model, this->total_force); Real time_step_suggested = this->model->getStableTimeStep() / sqrt(code_unit_system.e_m2dd_tt); DUMP("Akantu's critical timestep evaluation: " << time_step_suggested, DBG_MESSAGE); if (this->timeStep > time_step_suggested) LM_FATAL("Required timestep larger than critical timestep"); this->model->setTimeStep(this->timeStep); this->checkBoundaryInputSanity(); this->enableNeumanBndyConditions(); } /* ------------------------------------------------------------------------ */ // template void DomainAkantu::performStep1() { // STARTTIMER("AkantuStepPosition"); // // if (this->prestressed_periodicity_flag && // // !this->slave_displacements_computed) { // // this->buildSlaveDisplacements(); // // this->slave_displacements_computed = true; // //} // if ((this->haveToApplyNeumann()) && !this->boundary_has_been_applied && // (current_step >= this->boundary_start)) { // this->enableNeumanBndyConditions(); // } // UInt nb_nodes = this->model->getFEEngine().getMesh().getNbNodes(); // for (UInt n = 0; n < nb_nodes; ++n) { // for (UInt i = 0; i < Dim; ++i) { // this->getV()[Dim * n + i] += // code_unit_system.ft_m2v * .5 * this->getInternalForce()[Dim * n + // i] / this->getMass()[Dim * n + i] * this->timeStep; // LM_ASSERT(!isnan(this->getV()[Dim * n + i]) && // !isinf(this->getV()[Dim * n + i]), // "problem with velocity on node " // << n << " Dim " << i << " " // << this->getInternalForce()[Dim * n + i] << " " // << this->getMass()[Dim * n + i]); // this->getU()[Dim * n + i] += this->timeStep * this->getV()[Dim * n + // i]; // } // } // STOPTIMER("AkantuStepPosition"); // #ifdef AKANTU_HEAT_TRANSFER // if (heat_transfer_flag) // heat_transfer_model->explicitPred(); // #endif // LM_TOIMPLEMENT; // // this->mesh_container.incRelease(); // } /* -------------------------------------------------------------------------- */ // template void DomainAkantu::performStep3() { // STARTTIMER("AkantuStepVelocity"); // UInt nb_nodes = this->model->getMesh().getNbNodes(); // for (UInt n = 0; n < nb_nodes; ++n) { // for (UInt i = 0; i < Dim; ++i) { // this->getV()[Dim * n + i] += // code_unit_system.ft_m2v * .5 * this->getInternalForce()[Dim * n + // i] / this->getMass()[Dim * n + i] * this->timeStep; // LM_ASSERT(!isnan(this->getV()[Dim * n + i]) && // !isinf(this->getV()[Dim * n + i]), // "problem with velocity on node " // << n << " Dim " << i << " " << this->getV()[Dim * n + // i] // << " " << this->getInternalForce()[Dim * n + i] << " " // << this->getMass()[Dim * n + i] << " " << // this->timeStep); // } // } // STOPTIMER("AkantuStepVelocity"); // #ifdef AKANTU_HEAT_TRANSFER // if (heat_transfer_flag) // heat_transfer_model->explicitCorr(); // #endif // LM_TOIMPLEMENT; // // this->mesh_container.incRelease(); // } /* ------------------------------------------------------------------------ */ template UInt DomainAkantu::getCurrentStep() { LM_TOIMPLEMENT; } /* ------------------------------------------------------------------------ */ // template void DomainAkantu::buildSlaveDisplacements() { // typedef std::vector>::iterator pair_it; // pair_it pbc_pair = this->pbc_pairs.begin(); // pair_it end = this->pbc_pairs.end(); // // upper bounds are masters, lower bounds slaves // for (; pbc_pair != end; ++pbc_pair) { // for (UInt direction = 0; direction < Dim; ++direction) { // // this->slave_displacements.push_back( // // this->getU()[Dim * pbc_pair->first + direction] - // // this->getU()[Dim * pbc_pair->second + direction]); // } // } // } /* ------------------------------------------------------------------------ */ template void DomainAkantu::checkBoundaryInputSanity() { Real realmax = std::numeric_limits::max(); this->boundary_is_traction = ((this->boundary_stress[0] != realmax) && (this->boundary_stress[Dim] == realmax)); this->boundary_is_stress = (this->boundary_stress[Dim] != realmax); if (this->boundary_is_traction && this->boundary_is_stress) { LM_FATAL("you cannot specify both a traction (TRACTION_BNDY_VECTOR) and " << "a stress (TRACTION_BNDY_TENSOR) for traction " << "boundary conditions"); } // } else if (!this->boundary_is_traction && !boundary_is_stress) { // LM_FATAL("you specified a TRACTION_BNDY_TAG, but neither a stress " // << "(TRACTION_BNDY_TENSOR) nor a traction " // << "(TRACTION_BNDY_VECTOR)"); // } DUMP("boundary_is_traction: " << this->boundary_is_traction, DBG_MESSAGE); DUMP("boundary_is_stress: " << this->boundary_is_stress, DBG_MESSAGE); } /* ------------------------------------------------------------------------ */ template void DomainAkantu::enableNeumanBndyConditions() { try { if (this->boundary_is_traction) { akantu::Vector surface_traction(Dim); for (UInt i = 0; i < Dim; ++i) { surface_traction(i) = this->boundary_stress[i]; } model->applyBC(akantu::BC::Neumann::FromSameDim(surface_traction), surface_tag); } else if (this->boundary_is_stress) { akantu::Matrix surface_stress(Dim, Dim); for (UInt i = 0; i < Dim; ++i) { for (UInt j = 0; j < Dim; ++j) { surface_stress(i, j) = this->boundary_stress[Dim * i + j]; } } model->applyBC(akantu::BC::Neumann::FromHigherDim(surface_stress), surface_tag); } } catch (akantu::debug::Exception &e) { LM_FATAL(e.what()); } this->boundary_has_been_applied = true; } /* -------------------------------------------------------------------------- */ template void DomainAkantu::setTimeStep(Real ts) { model->setTimeStep(ts); this->timeStep = ts; } /* -------------------------------------------------------------------------- */ template void DomainAkantu::setCurrentStep(UInt step) { current_step = step; } /* -------------------------------------------------------------------------- */ template Real DomainAkantu::getEpot() { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ template Cube DomainAkantu::getSubDomainBoundingBox() const { akantu::Vector _xmin; akantu::Vector _xmax; _xmin = mesh->getLocalLowerBounds(); _xmax = mesh->getLocalUpperBounds(); Cube c; for (UInt k = 0; k < Dim; ++k) { c.getXmin()[k] = _xmin[k]; c.getXmax()[k] = _xmax[k]; } return c; } /* -------------------------------------------------------------------------- */ template Cube DomainAkantu::getDomainBoundingBox() const { auto _xmin = akantu::Vector::zeros(Dim); auto _xmax = akantu::Vector::zeros(Dim); if (mesh) { _xmin = mesh->getLowerBounds(); _xmax = mesh->getUpperBounds(); } Cube c; for (UInt k = 0; k < Dim; ++k) { c.getXmin()[k] = _xmin[k]; c.getXmax()[k] = _xmax[k]; } return c; } /* -------------------------------------------------------------------------- */ template void DomainAkantu::initHeatTransfer() { #ifndef AKANTU_HEAT_TRANSFER LM_FATAL("Akantu was not compiled using the heat transfer model"); #else heat_transfer_model = new akantu::HeatTransferModel(*this->mesh, Dim, this->getID() + "heat"); // initialize PBC bool pbc_activated = false; for (UInt i = 0; i < Dim; ++i) { pbc_activated |= this->pbc[i]; } if (pbc_activated) { heat_transfer_model->setPBC(this->pbc[0], this->pbc[1], this->pbc[2]); std::map &akantu_pairs = heat_transfer_model->getPBCPairs(); std::map::iterator it = akantu_pairs.begin(); std::map::iterator end = akantu_pairs.end(); while (it != end) { this->pbc_pairs.push_back((*it)); ++it; } } // initialize everything if (heat_transfer_material_filename != "") { this->akantu_parser.parse(this->heat_transfer_material_filename); heat_transfer_model->setParser(this->akantu_parser); heat_transfer_model->initFull(); } else LM_FATAL("Heat transfer material file not provided."); this->elems.setHeatTransferModel(*heat_transfer_model); this->nodes.setHeatTransferModel(*heat_transfer_model); this->elems.setHeatTransferFlag(heat_transfer_flag); this->nodes.setHeatTransferFlag(heat_transfer_flag); // assemble the lumped capacity heat_transfer_model->assembleCapacityLumped(); // get stable time step akantu::Real suggested_time_step_heat_transfer = heat_transfer_model->getStableTimeStep() * 0.8; DUMP("suggested timestep for heat transfer model is " << suggested_time_step_heat_transfer, DBG_MESSAGE); if (this->timeStep > suggested_time_step_heat_transfer) LM_FATAL("provided timestep is larger than the critical timestep for the " "heat transfer model: please adapt it"); heat_transfer_model->setTimeStep(this->timeStep); /* -------------------- */ // initialize the vectors /* -------------------- */ temperature = new VecAkantu(heat_transfer_model->getTemperature()); heat_flux = new VecAkantu(heat_transfer_model->getExternalHeatRate()); temperature_rate = new VecAkantu(heat_transfer_model->getTemperatureRate()); heat_transfer_model->updateResidual(); akantu::Array &boundary = heat_transfer_model->getBlockedDOFs(); UInt nb_nodes = boundary.getSize(); bool *bound = boundary.storage(); if (heat_boundary_geom != invalidGeom) { Geometry *geom_dirichlet = GeometryManager::getManager().getGeometry(heat_boundary_geom); for (UInt n = 0; n < nb_nodes; ++n, ++bound) { Real pos[Dim]; this->nodes.get(n).getPositions0(pos); if (geom_dirichlet->contains(pos)) { *bound = true; } } } #endif } /* -------------------------------------------------------------------------- */ template void DomainAkantu::updateGradient() { if (!this->getCommGroup().amIinGroup()) return; STARTTIMER("AkantuStepForce"); auto &model = *this->model; auto &dof_manager = model.getDOFManager(); auto &residual = dof_manager.getResidual(); residual.zero(); static_cast(model).assembleResidual(); dof_manager.getArrayPerDOFs("displacement", residual, total_force); // this->model->assembleInternalForces(); STOPTIMER("AkantuStepForce"); // #ifdef AKANTU_HEAT_TRANSFER // if (heat_transfer_flag) // heat_transfer_model->updateResidual(); // #endif this->changeRelease(); } /* -------------------------------------------------------------------------- */ template void DomainAkantu::updateAcceleration() { if (!this->getCommGroup().amIinGroup()) return; auto &acceleration = model->getAcceleration(); auto &lumped_mass = model->getMass(); for (UInt i = 0; i < acceleration.size(); ++i) { for (UInt k = 0; k < Dim; ++k) { Real f = total_force(i, k); acceleration(i, k) = code_unit_system.ft_m2v * f / lumped_mass(i, k); } } } /* -------------------------------------------------------------------------- */ template Real DomainAkantu::potentialEnergy() { return model->getEnergy("potential"); } /* -------------------------------------------------------------------------- */ template void DomainAkantu::enforceCompatibility() { auto &blocked_dofs = model->getBlockedDOFs(); auto &acceleration = model->getAcceleration(); auto &velocity = model->getVelocity(); for (UInt i = 0; i < acceleration.size(); ++i) { for (UInt k = 0; k < Dim; ++k) { if (blocked_dofs(i, k)) { acceleration(i, k) = 0.; velocity(i, k) = 0.; } } } } /* -------------------------------------------------------------------------- */ template void DomainAkantu::setDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ template void DomainAkantu::setSubDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; } /* -------------------------------------------------------------------------- */ template ArrayView DomainAkantu::getField(FieldType field_type) { switch (field_type) { case _position0: return make_array_from_akantu(mesh->getNodes()); case _velocity: return make_array_from_akantu(model->getVelocity()); case _displacement: return make_array_from_akantu(model->getDisplacement()); default: LM_FATAL("Not accessible field: " << field_type); } } /* -------------------------------------------------------------------------- */ // template void DomainAkantu::performStep2() { // STARTTIMER("AkantuStepForce"); // LM_TOIMPLEMENT; // // this->model->initializeUpdateResidualData(); // // if (this->prestressed_periodicity_flag) { // // typedef std::vector>::iterator pair_it; // // pair_it pbc_pair = this->pbc_pairs.begin(); // // pair_it end = this->pbc_pairs.end(); // // UInt node_counter = 0; // // for (; pbc_pair != end; ++pbc_pair, ++node_counter) { // // for (UInt direction = 0; direction < Dim; ++direction) { // // this->getU()[Dim * pbc_pair->first + direction] += // // this->slave_displacements[Dim * node_counter + direction]; // // } // // } // // } // LM_TOIMPLEMENT; // // this->model->updateResidual(false); // // this->model->updateAcceleration(); // STOPTIMER("AkantuStepForce"); // #ifdef AKANTU_HEAT_TRANSFER // if (heat_transfer_flag) // heat_transfer_model->updateResidual(); // #endif // LM_TOIMPLEMENT; // // this->mesh_container.incRelease(); // } /* -------------------------------------------------------------------------- */ /* LMDESC AKANTU_BASE This domain implements the plugin with Akantu Finite Element library. */ /* LMHERITANCE domain_continuum */ template void DomainAkantu::declareParams() { Real realmax = std::numeric_limits::max(); DomainContinuum, ContainerNodesAkantu, Dim>::declareParams(); /* LMKEYWORD MESH_FILENAME Specify the mesh file to be loaded using the MeshIO utility of Akantu */ this->parseKeyword("MESH_FILENAME", mesh_filename, ""); /* LMKEYWORD MATERIAL_FILENAME Specify the material file to be loaded for the Akantu SolidMechanicsModel object */ this->parseKeyword("MATERIAL_FILENAME", material_filename, ""); /* LMKEYWORD PBC Specify the directions in which Periodic Boundary Conditions should be activated */ this->parseVectorKeyword("PBC", Dim, pbc, VEC_DEFAULTS(0u, 0u, 0u)); /* LMKEYWORD SCALE Specify factors for each directions to be used in rescaling the mesh nodal positions */ this->parseVectorKeyword("SCALE", Dim, scale, VEC_DEFAULTS(1., 1., 1.)); /* LMKEYWORD SHIFT Specify factors for each directions to be used in shifting the mesh nodal positions */ this->parseVectorKeyword("SHIFT", Dim, shift, VEC_DEFAULTS(0., 0., 0.)); /* LMKEYWORD R0 In the case of 1D mesh allows to use a linearized LJ material definition fitting monoatomic chain: This provides the interatomic distance. */ this->parseKeyword("R0", r0, 0.); /* LMKEYWORD RCUT In the case of 1D mesh allows to use a linearized LJ material definition fitting monoatomic chain: This provides the cutoff radius for the potential evaluation. */ this->parseKeyword("RCUT", rcut, 0.); /* LMKEYWORD EPSILON In the case of 1D mesh allows to use a linearized LJ material definition fitting monoatomic chain: This provides the :math:`\epsilon` parameter in the LJ potential definition. */ this->parseKeyword("EPSILON", epsilon, 0.); /* LMKEYWORD SIGMA In the case of 1D mesh allows to use a linearized LJ material definition fitting monoatomic chain: This provides the :math:`\sigma` parameter in the LJ potential definition. */ this->parseKeyword("SIGMA", sigma, 0.); /* LMKEYWORD MASS In the case of 1D mesh allows to use a linearized LJ material definition fitting monoatomic chain: This provides the atomic mass. */ this->parseKeyword("MASS", atomic_mass, 0.); /* LMKEYWORD ELEM_SIZE In the case of 1D mesh generates automatically a homogeneous 1D first order mesh. */ this->parseKeyword("ELEM_SIZE", elem_size, 0.); /* LMKEYWORD TRACTION_BNDY_TENSOR Specify a full stress tensor to be applied to a boundary */ this->parseVectorKeyword( "TRACTION_BNDY_TENSOR", Dim * Dim, this->boundary_stress, VEC_DEFAULTS(realmax, realmax, realmax, realmax, realmax, realmax, realmax, realmax, realmax)); // /* L2MKEYWORD TRACTION_BNDY_VECTOR // Specify a traction vector to be applied to a boundary */ // this->parseVectorKeyword("TRACTION_BNDY_VECTOR", Dim, // this->boundary_stress, // {0.,0.,0.}); /* LMKEYWORD TRACTION_BNDY_TAG Specify the tag of the boundary on which the traction should be applied. This refers to :math:`tag_1` in the msh format */ this->parseKeyword("TRACTION_BNDY_TAG", this->surface_tag, ""); /* LMKEYWORD TRACTION_BNDY_START Specify the time step at which to start applying the traction boundary conditions */ this->parseKeyword("TRACTION_BNDY_START", this->boundary_start, 0u); /* LMKEYWORD HEAT_TRANSFER To switch on the heat transfer model object */ this->parseTag("HEAT_TRANSFER", this->heat_transfer_flag, false); /* LMKEYWORD HEAT_TRANSFER_MATERIAL_FILENAME Specify the material file to be loaded for the Akantu HeatTransferModel object */ this->parseKeyword("HEAT_TRANSFER_MATERIAL_FILENAME", heat_transfer_material_filename, ""); /* LMKEYWORD HEAT_BOUNDARY Specify the geometry for which a Dirichlet boundary condition should be used for the HeatTransfer model */ this->parseKeyword("HEAT_BOUNDARY", heat_boundary_geom, invalidGeom); } /* -------------------------------------------------------------------------- */ template class DomainAkantu<3>; template class DomainAkantu<2>; template class DomainAkantu<1>; __END_LIBMULTISCALE__ diff --git a/src/stimulation/stimulation_zero.cc b/src/stimulation/stimulation_zero.cc index 771e5e8..042ca4b 100644 --- a/src/stimulation/stimulation_zero.cc +++ b/src/stimulation/stimulation_zero.cc @@ -1,151 +1,150 @@ /** * @file stimulation_raz.cc * * @author Guillaume Anciaux * @author Till Junge * * @date Wed Jul 09 21:59:47 2014 * * @brief Remise a zero (reset to zero) stimulator * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * LibMultiScale is free software: you can redistribute it and/or modify it * under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * LibMultiScale is distributed in the hope that it will be useful, but * WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with LibMultiScale. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "stimulation_zero.hh" #include "lib_continuum.hh" #include "lib_dd.hh" #include "lib_md.hh" #include "lm_common.hh" #include "lm_communicator.hh" #include "ref_point_data.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /* -------------------------------------------------------------------------- */ StimulationZero::StimulationZero(const std::string &name) : LMObject(name), COEF(0), average_flag(false), direction(-1) {} /* -------------------------------------------------------------------------- */ StimulationZero::~StimulationZero() {} /* -------------------------------------------------------------------------- */ template void StimulationZero::stimulate(Cont &cont) { constexpr UInt Dim = Cont::Dim; std::vector directions; if (this->direction < 0) { for (UInt i = 0; i < Dim; ++i) { directions.push_back(i); } } else { directions.push_back(this->direction); } typedef std::vector::iterator iter; DUMP("Stimulation ZERO fields " << " COEF : " << COEF << " pointer this " << this, DBG_INFO); for (UInt i = 0; i < field.size(); ++i) { FieldType f = field[i]; for (auto &&at : cont) { if (f == _velocity) { for (iter i = directions.begin(); i != directions.end(); ++i) { at.velocity()[*i] += (COEF - 1) * at.velocity()[*i]; } } else if (f == _displacement) { for (iter i = directions.begin(); i != directions.end(); ++i) { - throw; - // at.displacement()[*i] += (COEF-1) * at.displacement()[*i]; + 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 - position0 : selects the initial position field. Actually it is relevant for atomic systems assigning the initial positions to the current one. This is useful when a restart file needs to be generated with current positions and a zero displacement field. This command can be called several times for each field that need to be altered. */ this->parseKeyword("FIELD", field); /* LMKEYWORD COEF This is effective only for displacement, velocity or force field. Instead of reset to zero it multiplies by the provided coefficient. It enables, for instance to perform velocity damping. */ this->parseKeyword("COEF", COEF, 0.); /* LMKEYWORD COMPONENT Just apply stimulation in direction COMPONENT. By default all directions are stimulated. */ this->parseKeyword("COMPONENT", direction, -1); } /* -------------------------------------------------------------------------- */ DECLARE_STIMULATION_MAKE_CALL(StimulationZero) __END_LIBMULTISCALE__