diff --git a/src/model/dof_manager.cc b/src/model/dof_manager.cc index 9f282345d..2af5c8269 100644 --- a/src/model/dof_manager.cc +++ b/src/model/dof_manager.cc @@ -1,605 +1,607 @@ /** * @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 * * @section LICENSE * * 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 "element_group.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), communicator(Communicator::getStaticCommunicator()) {} /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(Mesh & mesh, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), mesh(&mesh), local_system_size(0), pure_local_system_size(0), system_size(0), communicator(mesh.getCommunicator()) { this->mesh->registerEventHandler(*this, _ehp_dof_manager); } /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() = default; /* -------------------------------------------------------------------------- */ // void DOFManager::getEquationsNumbers(const ID &, Array &) { // AKANTU_TO_IMPLEMENT(); // } /* -------------------------------------------------------------------------- */ 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, const ElementType & type, const 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_vector_iterator conn_begin = // connectivity.begin(nb_nodes_per_element); // Array::const_vector_iterator conn_it = conn_begin; 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, const ElementType & type, const 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.clear(); 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, const ElementType & type, const 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.clear(); 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); } } /* -------------------------------------------------------------------------- */ DOFManager::DOFData::DOFData(const ID & dof_id) : support_type(_dst_generic), group_support("__mesh__"), dof(nullptr), blocked_dofs(nullptr), increment(nullptr), previous(nullptr), solution(0, 1, dof_id + ":solution"), local_equation_number(0, 1, dof_id + ":local_equation_number") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData::~DOFData() = default; /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManager::getNewDOFData(const ID & dof_id) { auto it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } std::unique_ptr dofs_storage = std::make_unique(dof_id); this->dofs[dof_id] = std::move(dofs_storage); return *dofs_storage; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { DOFData & dofs_storage = this->getDOFData(dof_id); dofs_storage.dof = &dofs_array; UInt nb_local_dofs = 0; UInt nb_pure_local = 0; const DOFSupportType & support_type = dofs_storage.support_type; switch (support_type) { case _dst_nodal: { UInt nb_nodes = 0; const ID & group = dofs_storage.group_support; NodeGroup * node_group = nullptr; if (group == "__mesh__") { nb_nodes = this->mesh->getNbNodes(); } else { node_group = &this->mesh->getElementGroup(group).getNodeGroup(); nb_nodes = node_group->size(); } nb_local_dofs = nb_nodes; AKANTU_DEBUG_ASSERT( dofs_array.size() == nb_local_dofs, "The array of dof is too shot to be associated to nodes."); for (UInt n = 0; n < nb_nodes; ++n) { UInt node = n; if (node_group) node = node_group->getNodes()(n); nb_pure_local += this->mesh->isLocalOrMasterNode(node) ? 1 : 0; } nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); 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."); } } this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; communicator.allReduce(nb_pure_local, SynchronizerOperation::_sum); this->system_size += nb_pure_local; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { DOFData & dofs_storage = this->getNewDOFData(dof_id); dofs_storage.support_type = support_type; this->registerDOFsInternal(dof_id, dofs_array); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group) { DOFData & dofs_storage = this->getNewDOFData(dof_id); dofs_storage.support_type = _dst_nodal; dofs_storage.group_support = support_group; this->registerDOFsInternal(dof_id, dofs_array); } /* -------------------------------------------------------------------------- */ 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; } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { auto it = this->dofs.begin(); auto end = this->dofs.end(); for (; it != end; ++it) { DOFData & dof_data = *it->second; dof_data.solution.resize(dof_data.dof->size() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(it->first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix) { SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } SparseMatrix & ret = *matrix; this->matrices[matrix_id] = std::move(matrix); return ret; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix Array & DOFManager::getNewLumpedMatrix(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_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } auto mtx = std::make_unique>(this->local_system_size, 1, matrix_id); this->lumped_matrices[matrix_id] = std::move(mtx); return *this->lumped_matrices[matrix_id]; } /* -------------------------------------------------------------------------- */ 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(); } /* -------------------------------------------------------------------------- */ Array & 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 Array & 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)); } /* -------------------------------------------------------------------------- */ /* 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 = 0; UInt nb_new_pure_local = 0; nb_new_local_dofs = nodes_list.size(); for (const auto & node : nodes_list) { nb_new_pure_local += this->mesh->isLocalOrMasterNode(node) ? 1 : 0; } const auto & dof_array = *dof_data.dof; nb_new_pure_local *= dof_array.getNbComponent(); nb_new_local_dofs *= dof_array.getNbComponent(); 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(dof_data.solution.size() + nb_new_local_dofs); + return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesAdded(const Array & nodes_list, const NewNodesEvent &) { 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); } } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array &, const NewElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved(const Array &, const ElementTypeMapArray &, const RemovedElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) {} /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/dof_manager.hh b/src/model/dof_manager.hh index 4539874d8..7bf9d46b8 100644 --- a/src/model/dof_manager.hh +++ b/src/model/dof_manager.hh @@ -1,476 +1,479 @@ /** * @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 * * @section LICENSE * * 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 "aka_memory.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; } // namespace akantu namespace akantu { class DOFManager : protected Memory, protected MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFManager(const ID & id = "dof_manager", const MemoryID & memory_id = 0); DOFManager(Mesh & mesh, const ID & id = "dof_manager", const MemoryID & memory_id = 0); ~DOFManager() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ private: /// common function to help registering dofs void registerDOFsInternal(const ID & dof_id, Array & dofs_array); 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 & group_support); /// 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, const Array & array_to_assemble, Real scale_factor = 1.) = 0; /// Assemble an array to the global lumped matrix array virtual void assembleToLumpedMatrix(const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor = 1.) = 0; /** * 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, const ElementType & type, const 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, const ElementType & type, const 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, const ElementType & type, const 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, const ElementType & type, const 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 assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) = 0; /// multiply the dofs by a matrix and assemble the result to the residual virtual void assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor = 1); /// 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; /// sets the residual to 0 virtual void clearResidual() = 0; /// sets the matrix to 0 virtual void clearMatrix(const ID & mtx) = 0; /// sets the lumped matrix to 0 virtual void clearLumpedMatrix(const ID & mtx) = 0; /// splits the solution storage from a global view to the per dof storages void splitSolutionPerDOFs(); /// extract a lumped matrix part corresponding to a given dof virtual void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) = 0; protected: /// 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) = 0; protected: /* ------------------------------------------------------------------------ */ /// register a matrix SparseMatrix & registerSparseMatrix(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); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the equation numbers corresponding to a dof_id. This might be used to /// access the matrix. inline const Array & getEquationsNumbers(const ID & dof_id) const; /// 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); /// 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 & dofs_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); + /* ------------------------------------------------------------------------ */ /* 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 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 Array & getNewLumpedMatrix(const ID & matrix_id); /// Get the lumped version of a given matrix const Array & getLumpedMatrix(const ID & matrix_id) const; /// Get the lumped version of a given matrix Array & 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) = 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) { return *mesh; } else { AKANTU_EXCEPTION("No mesh registered in this dof manager"); } } /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Communicator, communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, communicator, auto &); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ protected: /// helper function for the DOFManager::onNodesAdded method virtual std::pair updateNodalDOFs(const ID & dof_id, const Array & nodes_list); 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: struct DOFData; 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 DOFData & getNewDOFData(const ID & dof_id); /* ------------------------------------------------------------------------ */ /* 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; /// Blocked degree of freedoms array Array * blocked_dofs; /// Degree of freedoms increment Array * increment; /// Degree of freedoms at previous step Array * previous; /// Solution associated to the dof Array solution; /// local numbering equation numbers Array local_equation_number; /* ---------------------------------------------------------------------- */ /* data for dynamic simulations */ /* ---------------------------------------------------------------------- */ /// Degree of freedom derivatives arrays std::vector *> dof_derivatives; }; /// 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>; /// 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}; /// Communicator used for this manager, should be the same as in the mesh if a /// mesh is registered Communicator & communicator; }; using DefaultDOFManagerFactory = Factory; using DOFManagerFactory = Factory; } // namespace akantu #include "dof_manager_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_HH__ */ diff --git a/src/model/dof_manager_inline_impl.cc b/src/model/dof_manager_inline_impl.cc index 71ebf49fd..e82dd6155 100644 --- a/src/model/dof_manager_inline_impl.cc +++ b/src/model/dof_manager_inline_impl.cc @@ -1,154 +1,159 @@ /** * @file dof_manager_inline_impl.cc * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Jan 31 2018 * * @brief inline functions of the dof manager * * @section LICENSE * * Copyright (©) 2014-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" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_INLINE_IMPL_CC__ #define __AKANTU_DOF_MANAGER_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasDOFs(const ID & dof_id) const { auto it = this->dofs.find(dof_id); return it != this->dofs.end(); } /* -------------------------------------------------------------------------- */ inline DOFManager::DOFData & DOFManager::getDOFData(const ID & dof_id) { auto it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof " << dof_id << " does not exists in " << this->id); } return *it->second; } /* -------------------------------------------------------------------------- */ const DOFManager::DOFData & DOFManager::getDOFData(const ID & dof_id) const { auto it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof " << dof_id << " does not exists in " << this->id); } return *it->second; } /* -------------------------------------------------------------------------- */ const Array & DOFManager::getEquationsNumbers(const ID & dof_id) const { return getDOFData(dof_id).local_equation_number; } /* -------------------------------------------------------------------------- */ template inline _DOFData & DOFManager::getDOFDataTyped(const ID & dof_id) { return dynamic_cast<_DOFData &>(this->getDOFData(dof_id)); } /* -------------------------------------------------------------------------- */ template inline const _DOFData & DOFManager::getDOFDataTyped(const ID & dof_id) const { return dynamic_cast(this->getDOFData(dof_id)); } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFs(const ID & dofs_id) { return *(this->getDOFData(dofs_id).dof); } /* -------------------------------------------------------------------------- */ inline DOFSupportType DOFManager::getSupportType(const ID & dofs_id) const { return this->getDOFData(dofs_id).support_type; } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getPreviousDOFs(const ID & dofs_id) { return *(this->getDOFData(dofs_id).previous); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasPreviousDOFs(const ID & dofs_id) const { return (this->getDOFData(dofs_id).previous != nullptr); } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFsIncrement(const ID & dofs_id) { return *(this->getDOFData(dofs_id).increment); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasDOFsIncrement(const ID & dofs_id) const { return (this->getDOFData(dofs_id).increment != nullptr); } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFsDerivatives(const ID & dofs_id, UInt order) { std::vector *> & derivatives = this->getDOFData(dofs_id).dof_derivatives; if ((order > derivatives.size()) || (derivatives[order - 1] == nullptr)) AKANTU_EXCEPTION("No derivatives of order " << order << " present in " << this->id << " for dof " << dofs_id); return *derivatives[order - 1]; } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasDOFsDerivatives(const ID & dofs_id, UInt order) const { const std::vector *> & derivatives = this->getDOFData(dofs_id).dof_derivatives; return ((order < derivatives.size()) && (derivatives[order - 1] != nullptr)); } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getSolution(const ID & dofs_id) const { return this->getDOFData(dofs_id).solution; } +/* -------------------------------------------------------------------------- */ +inline Array & DOFManager::getSolution(const ID & dofs_id) { + return this->getDOFData(dofs_id).solution; +} + /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getBlockedDOFs(const ID & dofs_id) const { return *(this->getDOFData(dofs_id).blocked_dofs); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasBlockedDOFs(const ID & dofs_id) const { return (this->getDOFData(dofs_id).blocked_dofs != nullptr); } /* -------------------------------------------------------------------------- */ } // akantu #endif /* __AKANTU_DOF_MANAGER_INLINE_IMPL_CC__ */