diff --git a/src/model/common/dof_manager/dof_manager.cc b/src/model/common/dof_manager/dof_manager.cc index 783f2a658..ac135f4bb 100644 --- a/src/model/common/dof_manager/dof_manager.cc +++ b/src/model/common/dof_manager/dof_manager.cc @@ -1,1018 +1,1015 @@ /** * @file dof_manager.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the common parts of the DOFManagers * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager.hh" #include "communicator.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "node_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver.hh" #include "periodic_node_synchronizer.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id) : id(id), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(Communicator::getStaticCommunicator()) {} /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(Mesh & mesh, const ID & id) - : id(id), mesh(&mesh), - dofs_flag(0, 1, std::string(id + ":dofs_type")), + : id(id), mesh(&mesh), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(mesh.getCommunicator()) { this->mesh->registerEventHandler(*this, _ehp_dof_manager); } /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() = default; /* -------------------------------------------------------------------------- */ std::vector DOFManager::getDOFIDs() const { std::vector keys; for (const auto & dof_data : this->dofs) { keys.push_back(dof_data.first); } return keys; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_element; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_vect.size() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); Array::const_matrix_iterator elem_it = elementary_vect.begin(nb_degree_of_freedom, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++elem_it) { UInt element = el; if (filter_it != nullptr) { // conn_it = conn_begin + *filter_it; element = *filter_it; } // const Vector & conn = *conn_it; const Matrix & elemental_val = *elem_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_node = connectivity(element, n) * nb_degree_of_freedom; Vector assemble(array_assembeled.storage() + offset_node, nb_degree_of_freedom); Vector elem_val = elemental_val(n); assemble.aXplusY(elem_val, scale_factor); } if (filter_it != nullptr) { ++filter_it; } // else // ++conn_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToResidual( - const ID & dof_id, const Array & elementary_vect, - ElementType type, GhostType ghost_type, Real scale_factor, + const ID & dof_id, const Array & elementary_vect, ElementType type, + GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToResidual(dof_id, array_localy_assembeled, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, - const ID & lumped_mtx, ElementType type, - GhostType ghost_type, Real scale_factor, - const Array & filter_elements) { + const ID & lumped_mtx, ElementType type, GhostType ghost_type, + Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToLumpedMatrix(dof_id, array_localy_assembeled, lumped_mtx, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = *pair.second; this->assembleMatMulVectToResidual(dof_id, A_id, *dof_data.dof, scale_factor); } } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { for (auto && data : this->dofs) { auto & dof_data = *data.second; dof_data.solution.resize(dof_data.dof->size() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(data.first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ void DOFManager::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getSolution(), solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); this->assembleToGlobalArray(dof_id, array_to_assemble, this->getResidual(), scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); auto & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ DOFManager::DOFData::DOFData(const ID & dof_id) : support_type(_dst_generic), group_support("__mesh__"), solution(0, 1, dof_id + ":solution"), local_equation_number(0, 1, dof_id + ":local_equation_number"), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData::~DOFData() = default; /* -------------------------------------------------------------------------- */ template auto DOFManager::countDOFsForNodes(const DOFData & dof_data, UInt nb_nodes, Func && getNode) { auto nb_local_dofs = nb_nodes; decltype(nb_local_dofs) nb_pure_local = 0; for (auto n : arange(nb_nodes)) { UInt node = getNode(n); // http://www.open-std.org/jtc1/sc22/open/n2356/conv.html // bool are by convention casted to 0 and 1 when promoted to int nb_pure_local += this->mesh->isLocalOrMasterNode(node); nb_local_dofs -= this->mesh->isPeriodicSlave(node); } const auto & dofs_array = *dof_data.dof; nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); return std::make_pair(nb_local_dofs, nb_pure_local); } /* -------------------------------------------------------------------------- */ auto DOFManager::getNewDOFDataInternal(const ID & dof_id) -> DOFData & { auto it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } std::unique_ptr dof_data_ptr = this->getNewDOFData(dof_id); DOFData & dof_data = *dof_data_ptr; this->dofs[dof_id] = std::move(dof_data_ptr); return dof_data; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = support_type; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = _dst_nodal; dofs_storage.group_support = support_group; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ std::tuple DOFManager::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { DOFData & dof_data = this->getDOFData(dof_id); dof_data.dof = &dofs_array; UInt nb_local_dofs = 0; UInt nb_pure_local = 0; const auto & support_type = dof_data.support_type; switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; std::function getNode; if (group == "__mesh__") { AKANTU_DEBUG_ASSERT( dofs_array.size() == this->mesh->getNbNodes(), "The array of dof is too short to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes( dof_data, this->mesh->getNbNodes(), [](auto && n) { return n; }); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); AKANTU_DEBUG_ASSERT( dofs_array.size() == node_group.size(), "The array of dof is too shot to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes(dof_data, node_group.size(), [&node_group](auto && n) { return node_group(n); }); } break; } case _dst_generic: { nb_local_dofs = nb_pure_local = dofs_array.size() * dofs_array.getNbComponent(); break; } - default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } + default: { + AKANTU_EXCEPTION("This type of dofs is not handled yet."); + } } dof_data.local_nb_dofs = nb_local_dofs; dof_data.pure_local_nb_dofs = nb_pure_local; dof_data.ghosts_nb_dofs = nb_local_dofs - nb_pure_local; this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; auto nb_total_pure_local = nb_pure_local; communicator.allReduce(nb_total_pure_local, SynchronizerOperation::_sum); this->system_size += nb_total_pure_local; // updating the dofs data after counting is finished switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; if (group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_local_dofs, nb_pure_local, support_nodes.size(), [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local, mesh->getNbNodes(), [](UInt node) -> UInt { return node; }); } break; } case _dst_generic: { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local); break; } } return std::make_tuple(nb_local_dofs, nb_pure_local, nb_total_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsPrevious(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.previous != nullptr) { AKANTU_EXCEPTION("The previous dofs array for " << dof_id << " has already been registered"); } dof.previous = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsIncrement(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.increment != nullptr) { AKANTU_EXCEPTION("The dofs increment array for " << dof_id << " has already been registered"); } dof.increment = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative) { DOFData & dof = this->getDOFData(dof_id); std::vector *> & derivatives = dof.dof_derivatives; if (derivatives.size() < order) { derivatives.resize(order, nullptr); } else { if (derivatives[order - 1] != nullptr) { AKANTU_EXCEPTION("The dof derivatives of order " << order << " already been registered for this dof (" << dof_id << ")"); } } derivatives[order - 1] = &dofs_derivative; } /* -------------------------------------------------------------------------- */ void DOFManager::registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs) { DOFData & dof = this->getDOFData(dof_id); if (dof.blocked_dofs != nullptr) { AKANTU_EXCEPTION("The blocked dofs array for " << dof_id << " has already been registered"); } dof.blocked_dofs = &blocked_dofs; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->matrices[matrix_id] = std::move(matrix); return ret; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix SolverVector & DOFManager::registerLumpedMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->lumped_matrices.find(matrix_id); if (it != this->lumped_matrices.end()) { AKANTU_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->lumped_matrices[matrix_id] = std::move(matrix); ret.resize(); return ret; } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::registerNonLinearSolver( const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver) { NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it != this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " already exists in " << this->id); } NonLinearSolver & ret = *non_linear_solver; this->non_linear_solvers[non_linear_solver_id] = std::move(non_linear_solver); return ret; } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::registerTimeStepSolver( const ID & time_step_solver_id, std::unique_ptr & time_step_solver) { TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it != this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " already exists in " << this->id); } TimeStepSolver & ret = *time_step_solver; this->time_step_solvers[time_step_solver_id] = std::move(time_step_solver); return ret; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::getMatrix(const ID & id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it == this->matrices.end()) { AKANTU_SILENT_EXCEPTION("The matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasMatrix(const ID & id) const { ID mtx_id = this->id + ":mtx:" + id; auto it = this->matrices.find(mtx_id); return it != this->matrices.end(); } /* -------------------------------------------------------------------------- */ SolverVector & DOFManager::getLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const SolverVector & DOFManager::getLumpedMatrix(const ID & id) const { ID matrix_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasLumpedMatrix(const ID & id) const { ID mtx_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(mtx_id); return it != this->lumped_matrices.end(); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::getNonLinearSolver(const ID & id) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it == this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasNonLinearSolver(const ID & id) const { ID solver_id = this->id + ":nls:" + id; auto it = this->non_linear_solvers.find(solver_id); return it != this->non_linear_solvers.end(); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::getTimeStepSolver(const ID & id) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it == this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasTimeStepSolver(const ID & solver_id) const { ID time_step_solver_id = this->id + ":tss:" + solver_id; auto it = this->time_step_solvers.find(time_step_solver_id); return it != this->time_step_solvers.end(); } /* -------------------------------------------------------------------------- */ void DOFManager::savePreviousDOFs(const ID & dofs_id) { this->getPreviousDOFs(dofs_id).copy(this->getDOFs(dofs_id)); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroResidual() { this->residual->zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroMatrix(const ID & mtx) { this->getMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ /* Mesh Events */ /* -------------------------------------------------------------------------- */ std::pair DOFManager::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { auto & dof_data = this->getDOFData(dof_id); UInt nb_new_local_dofs; UInt nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = countDOFsForNodes(dof_data, nodes_list.size(), [&nodes_list](auto && n) { return nodes_list(n); }); this->pure_local_system_size += nb_new_pure_local; this->local_system_size += nb_new_local_dofs; UInt nb_new_global = nb_new_pure_local; communicator.allReduce(nb_new_global, SynchronizerOperation::_sum); this->system_size += nb_new_global; dof_data.solution.resize(local_system_size); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, nodes_list.size(), [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::resizeGlobalArrays() { // resize all relevant arrays this->residual->resize(); this->solution->resize(); this->data_cache->resize(); for (auto & lumped_matrix : lumped_matrices) { lumped_matrix.second->resize(); } for (auto & matrix : matrices) { matrix.second->clearProfile(); } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesAdded(const Array & nodes_list, const NewNodesEvent & /*unused*/) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) { continue; } const auto & group = dof_data.group_support; if (group == "__mesh__") { this->updateNodalDOFs(dof_id, nodes_list); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup(); Array new_nodes_list; for (const auto & node : nodes_list) { if (node_group.find(node) != UInt(-1)) { new_nodes_list.push_back(node); } } this->updateNodalDOFs(dof_id, new_nodes_list); } } this->resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor(DOFManager::DOFData & dof_data, DOFManager & dof_manager) : dof_data(dof_data), dof_manager(dof_manager) { for (auto && pair : zip(dof_data.local_equation_number, dof_data.associated_nodes)) { UInt node; Int dof; std::tie(dof, node) = pair; dofs_per_node[node].push_back(dof); } } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { return nodes.size() * dof_data.dof->getNbComponent() * sizeof(Int); } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node.at(node); for (const auto & dof : dofs) { buffer << dof_manager.global_equation_number(dof); } } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node[node]; for (const auto & dof : dofs) { Int global_dof; buffer >> global_dof; AKANTU_DEBUG_ASSERT( (dof_manager.global_equation_number(dof) == -1 or dof_manager.global_equation_number(dof) == global_dof), "This dof already had a global_dof_id which is different from " "the received one. " << dof_manager.global_equation_number(dof) << " != " << global_dof); dof_manager.global_equation_number(dof) = global_dof; dof_manager.global_to_local_mapping[global_dof] = dof; } } } } protected: std::unordered_map> dofs_per_node; DOFManager::DOFData & dof_data; DOFManager & dof_manager; }; /* -------------------------------------------------------------------------- */ auto DOFManager::computeFirstDOFIDs(UInt nb_new_local_dofs, UInt nb_new_pure_local) { // determine the first local/global dof id to use UInt offset = 0; this->communicator.exclusiveScan(nb_new_pure_local, offset); auto first_global_dof_id = this->first_global_dof_id + offset; auto first_local_dof_id = this->local_system_size - nb_new_local_dofs; offset = nb_new_pure_local; this->communicator.allReduce(offset); this->first_global_dof_id += offset; return std::make_pair(first_local_dof_id, first_global_dof_id); } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, UInt nb_node, const std::function & getNode) { auto nb_local_dofs_added = nb_node * dof_data.dof->getNbComponent(); auto first_dof_pos = dof_data.local_equation_number.size(); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_local_dofs_added); dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_local_dofs_added); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); std::unordered_map, UInt> masters_dofs; // update per dof info UInt local_eq_num; UInt first_global_dof_id; std::tie(local_eq_num, first_global_dof_id) = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); for (auto d : arange(nb_local_dofs_added)) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof_flag = this->mesh->getNodeFlag(node); dof_data.associated_nodes.push_back(node); auto is_local_dof = this->mesh->isLocalOrMasterNode(node); auto is_periodic_slave = this->mesh->isPeriodicSlave(node); auto is_periodic_master = this->mesh->isPeriodicMaster(node); if (is_periodic_slave) { dof_data.local_equation_number.push_back(-1); continue; } // update equation numbers this->dofs_flag(local_eq_num) = dof_flag; dof_data.local_equation_number.push_back(local_eq_num); if (is_local_dof) { this->global_equation_number(local_eq_num) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = local_eq_num; ++first_global_dof_id; } else { this->global_equation_number(local_eq_num) = -1; } if (is_periodic_master) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof = d % dof_data.dof->getNbComponent(); masters_dofs.insert( std::make_pair(std::make_pair(node, dof), local_eq_num)); } ++local_eq_num; } // correct periodic slave equation numbers if (this->mesh->isPeriodic()) { auto assoc_begin = dof_data.associated_nodes.begin(); for (auto d : arange(nb_local_dofs_added)) { auto node = dof_data.associated_nodes(first_dof_pos + d); if (not this->mesh->isPeriodicSlave(node)) { continue; } auto master_node = this->mesh->getPeriodicMaster(node); auto dof = d % dof_data.dof->getNbComponent(); dof_data.local_equation_number(first_dof_pos + d) = masters_dofs[std::make_pair(master_node, dof)]; } } // synchronize the global numbering for slaves nodes if (this->mesh->isDistributed()) { GlobalDOFInfoDataAccessor data_accessor(dof_data, *this); if (this->mesh->isPeriodic()) { mesh->getPeriodicNodeSynchronizer().synchronizeOnce( data_accessor, SynchronizationTag::_giu_global_conn); } auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, SynchronizationTag::_ask_nodes); } } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local) { dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); UInt first_local_dof_id; UInt first_global_dof_id; std::tie(first_local_dof_id, first_global_dof_id) = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); // update per dof info for (auto _ [[gnu::unused]] : arange(nb_new_local_dofs)) { // update equation numbers this->dofs_flag(first_local_dof_id) = NodeFlag::_normal; dof_data.local_equation_number.push_back(first_local_dof_id); this->global_equation_number(first_local_dof_id) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = first_local_dof_id; ++first_global_dof_id; ++first_local_dof_id; } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array & /*unused*/, const Array & /*unused*/, const RemovedNodesEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array & /*unused*/, const NewElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved(const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const RemovedElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::updateGlobalBlockedDofs() { this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.reserve(this->local_system_size, 0); this->previous_global_blocked_dofs_release = this->global_blocked_dofs_release; for (auto & pair : dofs) { - if (!this->hasBlockedDOFs(pair.first)) { + if (not this->hasBlockedDOFs(pair.first)) { continue; } DOFData & dof_data = *pair.second; for (auto && data : zip(dof_data.getLocalEquationsNumbers(), make_view(*dof_data.blocked_dofs))) { const auto & dof = std::get<0>(data); const auto & is_blocked = std::get<1>(data); if (is_blocked) { this->global_blocked_dofs.push_back(dof); } } } std::sort(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); auto last = std::unique(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); this->global_blocked_dofs.resize(last - this->global_blocked_dofs.begin()); auto are_equal = global_blocked_dofs.size() == previous_global_blocked_dofs.size() and std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (not are_equal) { ++this->global_blocked_dofs_release; } } /* -------------------------------------------------------------------------- */ void DOFManager::applyBoundary(const ID & matrix_id) { auto & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getRelease()) { - auto are_equal = this->global_blocked_dofs_release == - this->previous_global_blocked_dofs_release; - // std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), - // previous_global_blocked_dofs.begin()); - - if (not are_equal) { + if (this->hasBlockedDOFsChanged()) { J.applyBoundary(); } previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getRelease(); + this->previous_global_blocked_dofs_release = + this->global_blocked_dofs_release; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToGlobalArray(const ID & dof_id, const ID & A_id, const Array & x, SolverVector & array, Real scale_factor) { auto & A = this->getMatrix(A_id); data_cache->resize(); data_cache->zero(); this->assembleToGlobalArray(dof_id, x, *data_cache, 1.); A.matVecMul(*data_cache, array, scale_factor, 1.); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { assembleMatMulVectToGlobalArray(dof_id, A_id, x, *residual, scale_factor); } } // namespace akantu diff --git a/src/model/common/dof_manager/dof_manager.hh b/src/model/common/dof_manager/dof_manager.hh index 916370287..2014fb11b 100644 --- a/src/model/common/dof_manager/dof_manager.hh +++ b/src/model/common/dof_manager/dof_manager.hh @@ -1,715 +1,721 @@ /** * @file dof_manager.hh * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Class handling the different types of dofs * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_factory.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_DOF_MANAGER_HH_ #define AKANTU_DOF_MANAGER_HH_ namespace akantu { class TermsToAssemble; class NonLinearSolver; class TimeStepSolver; class SparseMatrix; class SolverVector; class SolverCallback; } // namespace akantu namespace akantu { class DOFManager : protected MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ protected: struct DOFData; public: DOFManager(const ID & id = "dof_manager"); DOFManager(Mesh & mesh, const ID & id = "dof_manager"); ~DOFManager() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register an array of degree of freedom virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type); /// the dof as an implied type of _dst_nodal and is defined only on a subset /// of nodes virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group); /// register an array of previous values of the degree of freedom virtual void registerDOFsPrevious(const ID & dof_id, Array & dofs_array); /// register an array of increment of degree of freedom virtual void registerDOFsIncrement(const ID & dof_id, Array & dofs_array); /// register an array of derivatives for a particular dof array virtual void registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative); /// register array representing the blocked degree of freedoms virtual void registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor = 1.); /// Assemble an array to the global lumped matrix array virtual void assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor = 1.); /** * Assemble elementary values to a local array of the size nb_nodes * * nb_dof_per_node. The dof number is implicitly considered as * conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to a global array corresponding to a lumped * matrix */ virtual void assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. With 0 < * n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type = _not_ghost, const MatrixType & elemental_matrix_type = _symmetric, const Array & filter_elements = empty_filter) = 0; /// multiply a vector by a matrix and assemble the result to the residual virtual void assembleMatMulVectToArray(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor = 1) = 0; /// multiply a vector by a lumped matrix and assemble the result to the /// residual virtual void assembleLumpedMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) = 0; /// assemble coupling terms between to dofs virtual void assemblePreassembledMatrix(const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) = 0; /// multiply a vector by a matrix and assemble the result to the residual virtual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1); /// multiply the dofs by a matrix and assemble the result to the residual virtual void assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor = 1); /// updates the global blocked_dofs array virtual void updateGlobalBlockedDofs(); /// sets the residual to 0 virtual void zeroResidual(); /// sets the matrix to 0 virtual void zeroMatrix(const ID & mtx); /// sets the lumped matrix to 0 virtual void zeroLumpedMatrix(const ID & mtx); virtual void applyBoundary(const ID & matrix_id = "J"); // virtual void applyBoundaryLumped(const ID & matrix_id = "J"); /// extract a lumped matrix part corresponding to a given dof virtual void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped); /// splits the solution storage from a global view to the per dof storages void splitSolutionPerDOFs(); private: /// dispatch the creation of the dof data and register it DOFData & getNewDOFDataInternal(const ID & dof_id); protected: /// common function to help registering dofs the return values are the add new /// numbers of local dofs, pure local dofs, and system size virtual std::tuple registerDOFsInternal(const ID & dof_id, Array & dofs_array); /// minimum functionality to implement per derived version of the DOFManager /// to allow the splitSolutionPerDOFs function to work virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array); /// fill a Vector with the equation numbers corresponding to the given /// connectivity static inline void extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & element_equation_number); /// Assemble a array to a global one void assembleMatMulVectToGlobalArray(const ID & dof_id, const ID & A_id, const Array & x, SolverVector & array, Real scale_factor = 1.); /// common function that can be called by derived class with proper matrice /// types template void assemblePreassembledMatrix_(Mat & A, const ID & dof_id_m, const ID & dof_id_n, const TermsToAssemble & terms); template void assembleElementalMatricesToMatrix_( Mat & A, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements); template void assembleMatMulVectToArray_(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the location type of a given dof inline bool isLocalOrMasterDOF(UInt local_dof_num); /// Answer to the question is a dof a slave dof ? inline bool isSlaveDOF(UInt local_dof_num); /// Answer to the question is a dof a slave dof ? inline bool isPureGhostDOF(UInt local_dof_num); /// tells if the dof manager knows about a global dof bool hasGlobalEquationNumber(Int global) const; /// return the local index of the global equation number inline Int globalToLocalEquationNumber(Int global) const; /// converts local equation numbers to global equation numbers; inline Int localToGlobalEquationNumber(Int local) const; /// get the array of dof types (use only if you know what you do...) inline NodeFlag getDOFFlag(Int local_id) const; + /// defines if the boundary changed + bool hasBlockedDOFsChanged() { + return this->global_blocked_dofs_release != + this->previous_global_blocked_dofs_release; + } + /// Global number of dofs AKANTU_GET_MACRO(SystemSize, this->system_size, UInt); /// Local number of dofs AKANTU_GET_MACRO(LocalSystemSize, this->local_system_size, UInt); /// Pure local number of dofs AKANTU_GET_MACRO(PureLocalSystemSize, this->pure_local_system_size, UInt); /// Retrieve all the registered DOFs std::vector getDOFIDs() const; /* ------------------------------------------------------------------------ */ /* DOFs and derivatives accessors */ /* ------------------------------------------------------------------------ */ /// Get a reference to the registered dof array for a given id inline Array & getDOFs(const ID & dofs_id); /// Get the support type of a given dof inline DOFSupportType getSupportType(const ID & dofs_id) const; /// are the dofs registered inline bool hasDOFs(const ID & dof_id) const; /// Get a reference to the registered dof derivatives array for a given id inline Array & getDOFsDerivatives(const ID & dofs_id, UInt order); /// Does the dof has derivatives inline bool hasDOFsDerivatives(const ID & dofs_id, UInt order) const; /// Get a reference to the blocked dofs array registered for the given id inline const Array & getBlockedDOFs(const ID & dofs_id) const; /// Does the dof has a blocked array inline bool hasBlockedDOFs(const ID & dofs_id) const; /// Get a reference to the registered dof increment array for a given id inline Array & getDOFsIncrement(const ID & dofs_id); /// Does the dof has a increment array inline bool hasDOFsIncrement(const ID & dofs_id) const; /// Does the dof has a previous array inline Array & getPreviousDOFs(const ID & dofs_id); /// Get a reference to the registered dof array for previous step values a /// given id inline bool hasPreviousDOFs(const ID & dofs_id) const; /// saves the values from dofs to previous dofs virtual void savePreviousDOFs(const ID & dofs_id); /// Get a reference to the solution array registered for the given id inline const Array & getSolution(const ID & dofs_id) const; /// Get a reference to the solution array registered for the given id inline Array & getSolution(const ID & dofs_id); /// Get the blocked dofs array AKANTU_GET_MACRO(GlobalBlockedDOFs, global_blocked_dofs, const Array &); /// Get the blocked dofs array AKANTU_GET_MACRO(PreviousGlobalBlockedDOFs, previous_global_blocked_dofs, const Array &); /* ------------------------------------------------------------------------ */ /* Matrices accessors */ /* ------------------------------------------------------------------------ */ /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) = 0; /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) = 0; /// Get the equation numbers corresponding to a dof_id. This might be used to /// access the matrix. inline const Array & getLocalEquationsNumbers(const ID & dof_id) const; protected: /// get the array of dof types (use only if you know what you do...) inline const Array & getDOFsAssociatedNodes(const ID & dof_id) const; protected: /* ------------------------------------------------------------------------ */ /// register a matrix SparseMatrix & registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix); /// register a lumped matrix (aka a Vector) SolverVector & registerLumpedMatrix(const ID & matrix_id, std::unique_ptr & matrix); /// register a non linear solver instantiated by a derived class NonLinearSolver & registerNonLinearSolver(const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver); /// register a time step solver instantiated by a derived class TimeStepSolver & registerTimeStepSolver(const ID & time_step_solver_id, std::unique_ptr & time_step_solver); template NonLinearSolver & registerNonLinearSolver(DMType & dm, const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; std::unique_ptr nls = std::make_unique( dm, type, non_linear_solver_id); return this->registerNonLinearSolver(non_linear_solver_id, nls); } template TimeStepSolver & registerTimeStepSolver(DMType & dm, const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique(dm, type, non_linear_solver, solver_callback, time_step_solver_id); return this->registerTimeStepSolver(time_step_solver_id, tss); } template SparseMatrix & registerSparseMatrix(DMType & dm, const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; std::unique_ptr sm = std::make_unique(dm, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } template SparseMatrix & registerSparseMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; auto & sm_to_copy = aka::as_type(this->getMatrix(matrix_to_copy_id)); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } template SolverVector & registerLumpedMatrix(DMType & dm, const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; std::unique_ptr sm = std::make_unique(dm, matrix_id); return this->registerLumpedMatrix(matrix_id, sm); } protected: virtual void makeConsistentForPeriodicity(const ID & dof_id, SolverVector & array) = 0; virtual void assembleToGlobalArray(const ID & dof_id, const Array & array_to_assemble, SolverVector & global_array, Real scale_factor) = 0; public: /// extract degrees of freedom (identified by ID) from a global solver array virtual void getArrayPerDOFs(const ID & dof_id, const SolverVector & global, Array & local) = 0; /// Get the reference of an existing matrix SparseMatrix & getMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasMatrix(const ID & matrix_id) const; /// Get an instance of a new lumped matrix virtual SolverVector & getNewLumpedMatrix(const ID & matrix_id) = 0; /// Get the lumped version of a given matrix const SolverVector & getLumpedMatrix(const ID & matrix_id) const; /// Get the lumped version of a given matrix SolverVector & getLumpedMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasLumpedMatrix(const ID & matrix_id) const; /* ------------------------------------------------------------------------ */ /* Non linear system solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver( const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type) = 0; /// get instance of a non linear solver virtual NonLinearSolver & getNonLinearSolver(const ID & nls_solver_id); /// check if the given solver exists bool hasNonLinearSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver virtual TimeStepSolver & getNewTimeStepSolver(const ID & time_step_solver_id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback) = 0; /// get instance of a time step solver virtual TimeStepSolver & getTimeStepSolver(const ID & time_step_solver_id); /// check if the given solver exists bool hasTimeStepSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ const Mesh & getMesh() { if (mesh != nullptr) { return *mesh; } AKANTU_EXCEPTION("No mesh registered in this dof manager"); } /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Communicator, communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, communicator, auto &); /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Solution, *(solution.get()), const auto &); AKANTU_GET_MACRO_NOT_CONST(Solution, *(solution.get()), auto &); AKANTU_GET_MACRO(Residual, *(residual.get()), const auto &); AKANTU_GET_MACRO_NOT_CONST(Residual, *(residual.get()), auto &); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ protected: friend class GlobalDOFInfoDataAccessor; /// helper function for the DOFManager::onNodesAdded method virtual std::pair updateNodalDOFs(const ID & dof_id, const Array & nodes_list); template auto countDOFsForNodes(const DOFData & dof_data, UInt nb_nodes, Func && getNode); void updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, UInt nb_nodes, const std::function & getNode); void updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local); auto computeFirstDOFIDs(UInt nb_new_local_dofs, UInt nb_new_pure_local); /// resize all the global information and takes the needed measure like /// cleaning matrices profiles virtual void resizeGlobalArrays(); public: /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /// function to implement to react on akantu::RemovedNodesEvent void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event) override; /// function to implement to react on akantu::NewElementsEvent void onElementsAdded(const Array & elements_list, const NewElementsEvent & event) override; /// function to implement to react on akantu::RemovedElementsEvent void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; /// function to implement to react on akantu::ChangedElementsEvent void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event) override; protected: inline DOFData & getDOFData(const ID & dof_id); inline const DOFData & getDOFData(const ID & dof_id) const; template inline DOFData_ & getDOFDataTyped(const ID & dof_id); template inline const DOFData_ & getDOFDataTyped(const ID & dof_id) const; virtual std::unique_ptr getNewDOFData(const ID & dof_id) = 0; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// dof representations in the dof manager struct DOFData { DOFData() = delete; explicit DOFData(const ID & dof_id); virtual ~DOFData(); /// DOF support type (nodal, general) this is needed to determine how the /// dof are shared among processors DOFSupportType support_type; ID group_support; /// Degree of freedom array Array * dof{nullptr}; /// Blocked degree of freedoms array Array * blocked_dofs{nullptr}; /// Degree of freedoms increment Array * increment{nullptr}; /// Degree of freedoms at previous step Array * previous{nullptr}; /// Solution associated to the dof Array solution; /* ---------------------------------------------------------------------- */ /* data for dynamic simulations */ /* ---------------------------------------------------------------------- */ /// Degree of freedom derivatives arrays std::vector *> dof_derivatives; /* ---------------------------------------------------------------------- */ /// number of dofs to consider locally for this dof id UInt local_nb_dofs{0}; /// Number of purely local dofs UInt pure_local_nb_dofs{0}; /// number of ghost dofs UInt ghosts_nb_dofs{0}; /// local numbering equation numbers Array local_equation_number; /// associated node for _dst_nodal dofs only Array associated_nodes; virtual Array & getLocalEquationsNumbers() { return local_equation_number; } }; /// type to store dofs information using DOFStorage = std::map>; /// type to store all the matrices using SparseMatricesMap = std::map>; /// type to store all the lumped matrices using LumpedMatricesMap = std::map>; /// type to store all the non linear solver using NonLinearSolversMap = std::map>; /// type to store all the time step solver using TimeStepSolversMap = std::map>; ID id; /// store a reference to the dof arrays DOFStorage dofs; /// list of sparse matrices that where created SparseMatricesMap matrices; /// list of lumped matrices LumpedMatricesMap lumped_matrices; /// non linear solvers storage NonLinearSolversMap non_linear_solvers; /// time step solvers storage TimeStepSolversMap time_step_solvers; /// reference to the underlying mesh Mesh * mesh{nullptr}; /// Total number of degrees of freedom (size with the ghosts) UInt local_system_size{0}; /// Number of purely local dofs (size without the ghosts) UInt pure_local_system_size{0}; /// Total number of degrees of freedom UInt system_size{0}; /// rhs to the system of equation corresponding to the residual linked to the /// different dofs std::unique_ptr residual; /// solution of the system of equation corresponding to the different dofs std::unique_ptr solution; /// a vector that helps internally to perform some tasks std::unique_ptr data_cache; /// define the dofs type, local, shared, ghost Array dofs_flag; /// equation number in global numbering Array global_equation_number; using equation_numbers_map = std::unordered_map; /// dual information of global_equation_number equation_numbers_map global_to_local_mapping; /// Communicator used for this manager, should be the same as in the mesh if a /// mesh is registered Communicator & communicator; /// accumulator to know what would be the next global id to use UInt first_global_dof_id{0}; /// Release at last apply boundary on jacobian UInt jacobian_release{0}; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array global_blocked_dofs; UInt global_blocked_dofs_release{0}; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array previous_global_blocked_dofs; UInt previous_global_blocked_dofs_release{0}; private: /// This is for unit testing friend class DOFManagerTester; }; using DefaultDOFManagerFactory = Factory; using DOFManagerFactory = Factory; } // namespace akantu #include "dof_manager_inline_impl.hh" #endif /* AKANTU_DOF_MANAGER_HH_ */ diff --git a/src/model/common/integration_scheme/generalized_trapezoidal.cc b/src/model/common/integration_scheme/generalized_trapezoidal.cc index 89f8d24b1..906f9322e 100644 --- a/src/model/common/integration_scheme/generalized_trapezoidal.cc +++ b/src/model/common/integration_scheme/generalized_trapezoidal.cc @@ -1,193 +1,195 @@ /** * @file generalized_trapezoidal.cc * * @author Nicolas Richart * * @date creation: Fri Oct 23 2015 * @date last modification: Wed Jan 31 2018 * * @brief implementation of inline functions * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "generalized_trapezoidal.hh" #include "aka_array.hh" #include "dof_manager.hh" #include "mesh.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ GeneralizedTrapezoidal::GeneralizedTrapezoidal(DOFManager & dof_manager, const ID & dof_id, Real alpha) : IntegrationScheme1stOrder(dof_manager, dof_id), alpha(alpha) { this->registerParam("alpha", this->alpha, alpha, _pat_parsmod, "The alpha parameter"); } /* -------------------------------------------------------------------------- */ void GeneralizedTrapezoidal::predictor(Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.size(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt d = 0; d < nb_degree_of_freedom; d++) { if (!(*blocked_dofs_val)) { *u_val += (1. - alpha) * delta_t * *u_dot_val; } u_val++; u_dot_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GeneralizedTrapezoidal::corrector(const SolutionType & type, Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); switch (type) { case _temperature: this->allCorrector<_temperature>(delta_t, u, u_dot, blocked_dofs, delta); break; case _temperature_rate: this->allCorrector<_temperature_rate>(delta_t, u, u_dot, blocked_dofs, delta); break; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real GeneralizedTrapezoidal::getTemperatureCoefficient( const SolutionType & type, Real delta_t) const { switch (type) { case _temperature: return 1.; case _temperature_rate: return alpha * delta_t; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ Real GeneralizedTrapezoidal::getTemperatureRateCoefficient( const SolutionType & type, Real delta_t) const { switch (type) { case _temperature: return 1. / (alpha * delta_t); case _temperature_rate: return 1.; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ template void GeneralizedTrapezoidal::allCorrector(Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.size(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real e = getTemperatureCoefficient(type, delta_t); Real d = getTemperatureRateCoefficient(type, delta_t); Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); Real * delta_val = delta.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { if (!(*blocked_dofs_val)) { *u_val += e * *delta_val; *u_dot_val += d * *delta_val; } u_val++; u_dot_val++; delta_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GeneralizedTrapezoidal::assembleJacobian(const SolutionType & type, Real delta_t) { AKANTU_DEBUG_IN(); SparseMatrix & J = this->dof_manager.getMatrix("J"); const SparseMatrix & M = this->dof_manager.getMatrix("M"); const SparseMatrix & K = this->dof_manager.getMatrix("K"); bool does_j_need_update = false; does_j_need_update |= M.getRelease() != m_release; does_j_need_update |= K.getRelease() != k_release; - if (!does_j_need_update) { + does_j_need_update |= this->dof_manager.hasBlockedDOFsChanged(); + + if (not does_j_need_update) { AKANTU_DEBUG_OUT(); return; } J.copyProfile(K); // J.zero(); Real c = this->getTemperatureRateCoefficient(type, delta_t); Real e = this->getTemperatureCoefficient(type, delta_t); J.add(M, e); J.add(K, c); m_release = M.getRelease(); k_release = K.getRelease(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/common/integration_scheme/newmark-beta.cc b/src/model/common/integration_scheme/newmark-beta.cc index 375e0a7dd..9340dc3fc 100644 --- a/src/model/common/integration_scheme/newmark-beta.cc +++ b/src/model/common/integration_scheme/newmark-beta.cc @@ -1,260 +1,262 @@ /** * @file newmark-beta.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Oct 23 2015 * @date last modification: Wed Jan 31 2018 * * @brief implementation of the newmark-@f$\beta@f$ integration scheme. This * implementation is taken from Méthodes numériques en mécanique des solides by * Alain Curnier \note{ISBN: 2-88074-247-1} * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "newmark-beta.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NewmarkBeta::NewmarkBeta(DOFManager & dof_manager, const ID & dof_id, Real alpha, Real beta) : IntegrationScheme2ndOrder(dof_manager, dof_id), beta(beta), alpha(alpha), k(0.), h(0.), m_release(0), k_release(0), c_release(0) { this->registerParam("alpha", this->alpha, alpha, _pat_parsmod, "The alpha parameter"); this->registerParam("beta", this->beta, beta, _pat_parsmod, "The beta parameter"); } /* -------------------------------------------------------------------------- */ /* * @f$ \tilde{u_{n+1}} = u_{n} + \Delta t \dot{u}_n + \frac{\Delta t^2}{2} * \ddot{u}_n @f$ * @f$ \tilde{\dot{u}_{n+1}} = \dot{u}_{n} + \Delta t \ddot{u}_{n} @f$ * @f$ \tilde{\ddot{u}_{n}} = \ddot{u}_{n} @f$ */ void NewmarkBeta::predictor(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.size(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); Real * u_dot_dot_val = u_dot_dot.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt d = 0; d < nb_degree_of_freedom; d++) { if (!(*blocked_dofs_val)) { Real dt_a_n = delta_t * *u_dot_dot_val; *u_val += (1 - k * alpha) * delta_t * *u_dot_val + (.5 - h * alpha * beta) * delta_t * dt_a_n; *u_dot_val = (1 - k) * *u_dot_val + (1 - h * beta) * dt_a_n; *u_dot_dot_val = (1 - h) * *u_dot_dot_val; } u_val++; u_dot_val++; u_dot_dot_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NewmarkBeta::corrector(const SolutionType & type, Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); switch (type) { case _acceleration: { this->allCorrector<_acceleration>(delta_t, u, u_dot, u_dot_dot, blocked_dofs, delta); break; } case _velocity: { this->allCorrector<_velocity>(delta_t, u, u_dot, u_dot_dot, blocked_dofs, delta); break; } case _displacement: { this->allCorrector<_displacement>(delta_t, u, u_dot, u_dot_dot, blocked_dofs, delta); break; } default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real NewmarkBeta::getAccelerationCoefficient(const SolutionType & type, Real delta_t) const { switch (type) { case _acceleration: return 1.; case _velocity: return 1. / (beta * delta_t); case _displacement: return 1. / (alpha * beta * delta_t * delta_t); default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ Real NewmarkBeta::getVelocityCoefficient(const SolutionType & type, Real delta_t) const { switch (type) { case _acceleration: return beta * delta_t; case _velocity: return 1.; case _displacement: return 1. / (alpha * delta_t); default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ Real NewmarkBeta::getDisplacementCoefficient(const SolutionType & type, Real delta_t) const { switch (type) { case _acceleration: return alpha * beta * delta_t * delta_t; case _velocity: return alpha * delta_t; case _displacement: return 1.; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ template void NewmarkBeta::allCorrector(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.size(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real c = getAccelerationCoefficient(type, delta_t); Real d = getVelocityCoefficient(type, delta_t); Real e = getDisplacementCoefficient(type, delta_t); Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); Real * u_dot_dot_val = u_dot_dot.storage(); Real * delta_val = delta.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { if (!(*blocked_dofs_val)) { *u_val += e * *delta_val; *u_dot_val += d * *delta_val; *u_dot_dot_val += c * *delta_val; } u_val++; u_dot_val++; u_dot_dot_val++; delta_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NewmarkBeta::assembleJacobian(const SolutionType & type, Real delta_t) { AKANTU_DEBUG_IN(); SparseMatrix & J = this->dof_manager.getMatrix("J"); const SparseMatrix & M = this->dof_manager.getMatrix("M"); const SparseMatrix & K = this->dof_manager.getMatrix("K"); bool does_j_need_update = false; does_j_need_update |= M.getRelease() != m_release; does_j_need_update |= K.getRelease() != k_release; if (this->dof_manager.hasMatrix("C")) { const SparseMatrix & C = this->dof_manager.getMatrix("C"); does_j_need_update |= C.getRelease() != c_release; } + does_j_need_update |= this->dof_manager.hasBlockedDOFsChanged(); + if (!does_j_need_update) { AKANTU_DEBUG_OUT(); return; } J.copyProfile(K); // J.zero(); Real c = this->getAccelerationCoefficient(type, delta_t); Real e = this->getDisplacementCoefficient(type, delta_t); if (!(e == 0.)) { // in explicit this coefficient is exactly 0. J.add(K, e); } J.add(M, c); m_release = M.getRelease(); k_release = K.getRelease(); if (this->dof_manager.hasMatrix("C")) { Real d = this->getVelocityCoefficient(type, delta_t); const SparseMatrix & C = this->dof_manager.getMatrix("C"); J.add(C, d); c_release = C.getRelease(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/common/integration_scheme/pseudo_time.cc b/src/model/common/integration_scheme/pseudo_time.cc index 6986cba06..4f5f78950 100644 --- a/src/model/common/integration_scheme/pseudo_time.cc +++ b/src/model/common/integration_scheme/pseudo_time.cc @@ -1,84 +1,89 @@ /** * @file pseudo_time.cc * * @author Nicolas Richart * * @date creation: Fri Feb 19 2016 * @date last modification: Wed Jan 31 2018 * * @brief Implementation of a really simple integration scheme * * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "pseudo_time.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PseudoTime::PseudoTime(DOFManager & dof_manager, const ID & dof_id) : IntegrationScheme(dof_manager, dof_id, 0), k_release(0) {} /* -------------------------------------------------------------------------- */ std::vector PseudoTime::getNeededMatrixList() { return {"K"}; } /* -------------------------------------------------------------------------- */ void PseudoTime::predictor(Real /*delta_t*/) {} /* -------------------------------------------------------------------------- */ void PseudoTime::corrector(const SolutionType & /*type*/, Real /*delta_t*/) { auto & us = this->dof_manager.getDOFs(this->dof_id); const auto & deltas = this->dof_manager.getSolution(this->dof_id); const auto & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); for (auto && tuple : zip(make_view(us), deltas, make_view(blocked_dofs))) { auto & u = std::get<0>(tuple); const auto & delta = std::get<1>(tuple); const auto & bld = std::get<2>(tuple); if (not bld) { u += delta; } } } /* -------------------------------------------------------------------------- */ void PseudoTime::assembleJacobian(const SolutionType & /*type*/, Real /*delta_t*/) { SparseMatrix & J = this->dof_manager.getMatrix("J"); const SparseMatrix & K = this->dof_manager.getMatrix("K"); - if (K.getRelease() == k_release) { + bool does_j_need_update = false; + does_j_need_update |= K.getRelease() != k_release; + does_j_need_update |= this->dof_manager.hasBlockedDOFsChanged(); + + if (not does_j_need_update) { + AKANTU_DEBUG_OUT(); return; } J.copyProfile(K); // J.zero(); J.add(K); k_release = K.getRelease(); } /* -------------------------------------------------------------------------- */ void PseudoTime::assembleResidual(bool /*is_lumped*/) {} /* -------------------------------------------------------------------------- */ } // namespace akantu