diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 1012c39e6..2837ce216 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,878 +1,881 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date Tue Aug 11 16:21:01 2015 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2010-2011 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_default.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver_default.hh" #include "sparse_matrix_aij.hh" #include "static_communicator.hh" #include "terms_to_assemble.hh" #include "time_step_solver_default.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = i; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { if (c_jcn >= c_irn) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addElementalMatrixToUnsymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), jacobian_release(0), global_equation_number(0, 1, "global_equation_number"), first_global_dof_id(0), synchronizer(nullptr) {} /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), jacobian_release(0), global_equation_number(0, 1, "global_equation_number"), first_global_dof_id(0), synchronizer(nullptr) { if (this->mesh->isDistributed()) this->synchronizer = new DOFSynchronizer(*this, this->id + ":dof_synchronizer", this->memory_id, this->mesh->getCommunicator()); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() { delete this->synchronizer; delete this->global_residual; } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = array_to_assemble.getSize() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.getSize() == nb_degree_of_freedoms, "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); typename Array::const_scalar_iterator arr_it = array_to_assemble.begin_reinterpret(nb_degree_of_freedoms); Array::const_scalar_iterator equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++arr_it, ++equ_it) { global_array(*equ_it) += scale_factor * (*arr_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id) : DOFData(dof_id), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManagerDefault::getNewDOFData(const ID & dof_id) { DOFDataDefault * dofs_storage = new DOFDataDefault(dof_id); this->dofs[dof_id] = dofs_storage; return *dofs_storage; } /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: typedef typename std::unordered_map>::size_type size_type; GlobalDOFInfoDataAccessor() = default; void addDOFToNode(UInt node, UInt dof) { dofs_per_node[node].push_back(dof); } UInt getNthDOFForNode(UInt nth_dof, UInt node) const { return dofs_per_node.find(node)->second[nth_dof]; } virtual UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const { if (tag == _gst_size) { return nodes.getSize() * sizeof(size_type); } if (tag == _gst_update) { UInt total_size = 0; for (auto node : nodes) { auto it = dofs_per_node.find(node); if (it != dofs_per_node.end()) total_size += CommunicationBuffer::sizeInBuffer(it->second); } return total_size; } return 0; } virtual void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const { for (auto node : nodes) { auto it = dofs_per_node.find(node); if (tag == _gst_size) { if (it != dofs_per_node.end()) { buffer << it->second.size(); } else { buffer << 0; } } else if (tag == _gst_update) { if (it != dofs_per_node.end()) buffer << it->second; } } } virtual void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) { for (auto node : nodes) { auto it = dofs_per_node.find(node); if (tag == _gst_size) { size_type size; buffer >> size; if (size != 0) dofs_per_node[node].resize(size); } else if (tag == _gst_update) { if (it != dofs_per_node.end()) buffer >> it->second; } } } protected: std::unordered_map> dofs_per_node; }; /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs) { // Count the number of pure local dofs per proc const StaticCommunicator * comm = nullptr; UInt prank = 0; UInt psize = 1; if (mesh) { comm = &mesh->getCommunicator(); prank = comm->whoAmI(); psize = comm->getNbProc(); } else { comm = &(StaticCommunicator::getStaticCommunicator()); } // access the relevant data to update DOFDataDefault & dof_data = this->getDOFDataTyped(dof_id); const DOFSupportType & support_type = dof_data.support_type; const ID & group = dof_data.group_support; GlobalDOFInfoDataAccessor data_accessor; // resize all relevant arrays this->residual.resize(this->local_system_size); this->dofs_type.resize(local_system_size); this->global_solution.resize(this->local_system_size); this->global_blocked_dofs.resize(this->local_system_size); this->previous_global_blocked_dofs.resize(this->local_system_size); this->global_equation_number.resize(this->local_system_size); dof_data.local_equation_number.resize(nb_dofs); // determine the first local/global dof id to use Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = nb_pure_local_dofs; comm->allGather(nb_dofs_per_proc); this->first_global_dof_id += std::accumulate( nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0); UInt first_dof_id = this->local_system_size - nb_dofs; const Array * support_nodes = nullptr; if (support_type == _dst_nodal) { if (group != "mesh") { support_nodes = &this->mesh->getElementGroup(group).getNodeGroup().getNodes(); } dof_data.associated_nodes.resize(nb_dofs); } // update per dof info for (UInt d = 0; d < nb_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number(d) = local_eq_num; bool is_local_dof = true; // determine the dof type switch (support_type) { case _dst_nodal: { UInt node = d / dof_data.dof->getNbComponent(); if (support_nodes) node = (*support_nodes)(node); this->dofs_type(local_eq_num) = this->mesh->getNodeType(node); dof_data.associated_nodes(d) = node; is_local_dof = this->mesh->isLocalOrMasterNode(node); if (is_local_dof) { data_accessor.addDOFToNode(node, first_global_dof_id); } + break; } case _dst_generic: { this->dofs_type(local_eq_num) = _nt_normal; break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } // update global id for local dofs if (is_local_dof) { this->global_equation_number(local_eq_num) = this->first_global_dof_id; this->global_to_local_mapping[this->first_global_dof_id] = local_eq_num; ++this->first_global_dof_id; } else { this->global_equation_number(local_eq_num) = 0; } } if (support_type == _dst_nodal) { auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, _gst_size); node_synchronizer.synchronizeOnce(data_accessor, _gst_update); std::vector counters(nb_dofs); for (UInt d = 0; d < nb_dofs; ++d) { UInt local_eq_num = first_dof_id + d; if (this->isSlaveDOF(local_eq_num)) { UInt node = d / dof_data.dof->getNbComponent(); UInt dof_count = counters[node]++; UInt global_dof_id = data_accessor.getNthDOFForNode(dof_count, node); this->global_equation_number(local_eq_num) = global_dof_id; this->global_to_local_mapping[global_dof_id] = local_eq_num; } } } // update the synchronizer if needed if (this->synchronizer) this->synchronizer->registerDOFs(dof_id); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, support_type); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, group_support); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrix * sm = new SparseMatrixAIJ(*this, matrix_type, matrix_id); this->registerSparseMatrix(matrix_id, *sm); return *sm; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrixAIJ & sm_to_copy = this->getMatrix(matrix_to_copy_id); SparseMatrix * sm = new SparseMatrixAIJ(sm_to_copy, matrix_id); this->registerSparseMatrix(matrix_id, *sm); return *sm; } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) { SparseMatrix & matrix = DOFManager::getMatrix(id); return dynamic_cast(matrix); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerDefault::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolver * nls = NULL; switch (type) { #if defined(AKANTU_IMPLICIT) case _nls_newton_raphson: case _nls_newton_raphson_modified: { nls = new NonLinearSolverNewtonRaphson(*this, type, non_linear_solver_id, this->memory_id); break; } case _nls_linear: { nls = new NonLinearSolverLinear(*this, type, non_linear_solver_id, this->memory_id); break; } #endif case _nls_lumped: { nls = new NonLinearSolverLumped(*this, type, non_linear_solver_id, this->memory_id); break; } default: AKANTU_EXCEPTION("The asked type of non linear solver is not supported by " "this dof manager"); } this->registerNonLinearSolver(non_linear_solver_id, *nls); return *nls; } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolver * tss = new TimeStepSolverDefault( *this, type, non_linear_solver, time_step_solver_id, this->memory_id); this->registerTimeStepSolver(time_step_solver_id, *tss); return *tss; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearResidual() { this->residual.resize(this->local_system_size); this->residual.clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearMatrix(const ID & mtx) { this->getMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateGlobalBlockedDofs() { DOFStorage::iterator it = this->dofs.begin(); DOFStorage::iterator end = this->dofs.end(); this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.resize(this->local_system_size); this->global_blocked_dofs.clear(); for (; it != end; ++it) { - if (!this->hasBlockedDOFs(it->first)) continue; + if (!this->hasBlockedDOFs(it->first)) + continue; DOFData & dof_data = *it->second; this->assembleToGlobalArray(it->first, *dof_data.blocked_dofs, this->global_blocked_dofs, true); } } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = equation_number.getSize(); local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent()); auto loc_it = local_array.begin_reinterpret(nb_degree_of_freedoms); auto equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++loc_it, ++equ_it) { (*loc_it) = global_array(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getEquationsNumbers(const ID & dof_id, Array & equation_numbers) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->global_equation_number, equation_numbers); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->global_solution, solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::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 DOFManagerDefault::assembleToResidual( const ID & dof_id, const Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); this->assembleToGlobalArray(dof_id, array_to_assemble, this->residual, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToLumpedMatrix( const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); Array & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { SparseMatrixAIJ & A = this->getMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, 1.); A.matVecMul(data_cache, this->residual, scale_factor, 1.); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleLumpedMatMulVectToResidual( const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { const Array & A = this->getLumpedMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, scale_factor); Array::const_scalar_iterator A_it = A.begin(); Array::const_scalar_iterator A_end = A.end(); Array::const_scalar_iterator x_it = data_cache.begin(); Array::scalar_iterator r_it = this->residual.begin(); for (; A_it != A_end; ++A_it, ++x_it, ++r_it) { *r_it += *A_it * *x_it; } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); this->addToProfile(matrix_id, dof_id, type, ghost_type); DOFData & dof_data = this->getDOFData(dof_id); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt nb_element; if (ghost_type == _not_ghost) { nb_element = this->mesh->getNbElement(type); } else { AKANTU_DEBUG_TO_IMPLEMENT(); } UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filter_it = filter_elements.storage(); } else { if (dof_data.group_support != "mesh") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); nb_element = group_elements.getSize(); filter_it = group_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } } AKANTU_DEBUG_ASSERT(elementary_mat.getSize() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = dof_data.dof->getNbComponent(); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); auto conn_begin = connectivity.begin(nb_nodes_per_element); auto conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (filter_it) ++filter_it; else ++conn_it; if (A.getMatrixType() == _symmetric) if (elemental_matrix_type == _symmetric) this->addSymmetricElementalMatrixToSymmetric( A, *el_mat_it, element_eq_nb, A.getSize()); else this->addUnsymmetricElementalMatrixToSymmetric( A, *el_mat_it, element_eq_nb, A.getSize()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb, A.getSize()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { const Array & equation_number_m = this->getLocalEquationNumbers(dof_id_m); const Array & equation_number_n = this->getLocalEquationNumbers(dof_id_n); SparseMatrixAIJ & A = this->getMatrix(matrix_id); for (const auto & term : terms) { - A.addToMatrix(equation_number_m(term.i()), equation_number_n(term.j()), - term); + UInt gi = this->localToGlobalEquationNumber(equation_number_m(term.i())); + UInt gj = this->localToGlobalEquationNumber(equation_number_n(term.j())); + A.addToMatrix(gi, gj, term); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const DOFData & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) return; auto mat_dof = std::make_pair(matrix_id, dof_id); auto type_pair = std::make_pair(type, ghost_type); auto prof_it = this->matrix_profiled_dofs.find(mat_dof); if (prof_it != this->matrix_profiled_dofs.end() && std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) != prof_it->second.end()) return; UInt nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt size = A.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto & connectivity = this->mesh->getConnectivity(type, ghost_type); auto cbegin = connectivity.begin(nb_nodes_per_element); auto cit = cbegin; UInt nb_elements = connectivity.getSize(); UInt * ge_it = nullptr; if (dof_data.group_support != "mesh") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); ge_it = group_elements.storage(); nb_elements = group_elements.getSize(); } UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; Vector element_eq_nb(size_mat); for (UInt e = 0; e < nb_elements; ++e) { if (ge_it) cit = cbegin + *ge_it; this->extractElementEquationNumber( equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (ge_it) ++ge_it; else ++cit; for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = element_eq_nb(i); if (c_irn < size) { for (UInt j = 0; j < size_mat; ++j) { UInt c_jcn = element_eq_nb(j); if (c_jcn < size) { A.addToProfile(c_irn, c_jcn); } } } } } this->matrix_profiled_dofs[mat_dof].push_back(type_pair); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::applyBoundary(const ID & matrix_id) { this->updateGlobalBlockedDofs(); SparseMatrixAIJ & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getValueRelease()) { Array::const_scalar_iterator it = global_blocked_dofs.begin(); Array::const_scalar_iterator end = global_blocked_dofs.end(); Array::const_scalar_iterator pit = previous_global_blocked_dofs.begin(); for (; it != end && *it == *pit; ++it, ++pit) ; if (it != end) J.applyBoundary(); } else { J.applyBoundary(); } this->jacobian_release = J.getValueRelease(); } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getGlobalResidual() { if (this->synchronizer) { if (not this->global_residual) { this->global_residual = new Array(0, 1, "global_residual"); } if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->global_residual->resize(this->system_size); this->synchronizer->gather(this->residual, *this->global_residual); } else { this->synchronizer->gather(this->residual); } return *this->global_residual; } else { return this->residual; } } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getResidual() const { return this->residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::setGlobalSolution(const Array & solution) { if (this->synchronizer) { if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->synchronizer->scatter(this->global_solution, solution); } else { this->synchronizer->scatter(this->global_solution); } } else { AKANTU_DEBUG_ASSERT(solution.getSize() == this->global_solution.getSize(), "Sequential call to this function needs the solution " "to be the same size as the global_solution"); this->global_solution.copy(solution); } } /* -------------------------------------------------------------------------- */ __END_AKANTU__ // LocalWords: dof dofs diff --git a/src/model/dof_manager_default.hh b/src/model/dof_manager_default.hh index dc5ac4eba..4384b5301 100644 --- a/src/model/dof_manager_default.hh +++ b/src/model/dof_manager_default.hh @@ -1,333 +1,336 @@ /** * @file dof_manager_default.hh * * @author Nicolas Richart * * @date Tue Aug 11 14:06:18 2015 * * @brief Default implementation of the dof manager * * @section LICENSE * * Copyright (©) 2010-2011 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_DEFAULT_HH__ #define __AKANTU_DOF_MANAGER_DEFAULT_HH__ namespace akantu { class SparseMatrixAIJ; class NonLinearSolverDefault; class TimeStepSolverDefault; class DOFSynchronizer; } namespace akantu { class DOFManagerDefault : public DOFManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFManagerDefault(const ID & id = "dof_manager_default", const MemoryID & memory_id = 0); DOFManagerDefault(Mesh & mesh, const ID & id = "dof_manager_default", const MemoryID & memory_id = 0); virtual ~DOFManagerDefault(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ private: void registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs); 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); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor = 1.); /// 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.); /** * Assemble elementary values to the global matrix. 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, const MatrixType & elemental_matrix_type, const Array & filter_elements); /// 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 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); /// 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); protected: /// Assemble an array to the global residual array template void assembleToGlobalArray(const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor); public: /// clear the residual virtual void clearResidual(); /// sets the matrix to 0 virtual void clearMatrix(const ID & mtx); /// sets the lumped matrix to 0 virtual void clearLumpedMatrix(const ID & mtx); /// update the global dofs vector virtual void updateGlobalBlockedDofs(); /// apply boundary conditions to jacobian matrix virtual void applyBoundary(const ID & matrix_id = "J"); virtual void getEquationsNumbers(const ID & dof_id, Array & equation_numbers); protected: /// Get local part of an array corresponding to a given dofdata template void getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const; /// Get the part of the solution corresponding to the dof_id virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array); /// fill a Vector with the equation numbers corresponding to the given /// connectivity inline void extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & local_equation_number); public: /// extract a lumped matrix part corresponding to a given dof virtual void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped); private: /// Add a symmetric matrices to a symmetric sparse matrix void addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a unsymmetric matrices to a symmetric sparse matrix (i.e. cohesive /// elements) void addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a matrices to a unsymmetric sparse matrix void addElementalMatrixToUnsymmetric(SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); void addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type); /// 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); /// Get the reference of an existing matrix SparseMatrixAIJ & getMatrix(const ID & matrix_id); /* ------------------------------------------------------------------------ */ /* Non Linear Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver(const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type); /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver TimeStepSolver & getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver); /* ------------------------------------------------------------------------ */ /// Get the solution array AKANTU_GET_MACRO_NOT_CONST(GlobalSolution, global_solution, Array &); /// Set the global solution array void setGlobalSolution(const Array & solution); /// Get the global residual array across processors (SMP call) const Array & getGlobalResidual(); /// Get the residual array const Array & getResidual() const; /// 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 &); /// Get the location type of a given dof inline bool isLocalOrMasterDOF(UInt dof_num); /// Answer to the question is a dof a slave dof ? inline bool isSlaveDOF(UInt dof_num); /// get the equation numbers (in local numbering) corresponding to a dof ID inline const Array & getLocalEquationNumbers(const ID & dof_id) const; + /// tells if the dof manager knows about a global dof + bool hasGlobalEquationNumber(UInt global) const; + /// return the local index of the global equation number inline UInt globalToLocalEquationNumber(UInt global) const; /// converts local equation numbers to global equation numbers; inline UInt localToGlobalEquationNumber(UInt local) const; /// get the array of dof types (use only if you know what you do...) inline Int getDOFType(UInt local_id) const; /// get the array of dof types (use only if you know what you do...) inline const Array & getDOFsAssociatedNodes(const ID & dof_id) const; /// access the internal dof_synchronizer AKANTU_GET_MACRO_NOT_CONST(Synchronizer, *synchronizer, DOFSynchronizer &); /// access the internal dof_synchronizer bool hasSynchronizer() const { return synchronizer != nullptr; } protected: virtual DOFData & getNewDOFData(const ID & dof_id); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: struct DOFDataDefault : public DOFData { DOFDataDefault(const ID & dof_id); /// associated node for _dst_nodal dofs only Array associated_nodes; }; typedef std::map AIJMatrixMap; typedef std::map DefaultNonLinearSolversMap; typedef std::map DefaultTimeStepSolversMap; typedef std::map, std::vector>> DOFToMatrixProfile; /// contains the the dofs that where added to the profile of a given matrix. DOFToMatrixProfile matrix_profiled_dofs; /// rhs to the system of equation corresponding to the residual linked to the /// different dofs Array residual; /// rhs used only on root proc in case of parallel computing, this is the full /// gathered rhs array Array * global_residual; /// solution of the system of equation corresponding to the different dofs Array global_solution; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array global_blocked_dofs; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array previous_global_blocked_dofs; /// define the dofs type, local, shared, ghost Array dofs_type; /// Map of the different matrices stored in the dof manager AIJMatrixMap aij_matrices; /// Map of the different time step solvers stored with there real type DefaultTimeStepSolversMap default_time_step_solver_map; /// Memory cache, this is an array to keep the temporary memory needed for /// some operations, it is meant to be resized or cleared when needed Array data_cache; /// Release at last apply boundary on jacobian UInt jacobian_release; /// equation number in global numbering Array global_equation_number; typedef unordered_map::type equation_numbers_map; /// dual information of global_equation_number equation_numbers_map global_to_local_mapping; /// accumulator to know what would be the next global id to use UInt first_global_dof_id; /// synchronizer to maintain coherency in dof fields DOFSynchronizer * synchronizer; }; } // akantu #include "dof_manager_default_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_DEFAULT_HH__ */ diff --git a/src/model/dof_manager_default_inline_impl.cc b/src/model/dof_manager_default_inline_impl.cc index 6216b6d99..b26b77f4c 100644 --- a/src/model/dof_manager_default_inline_impl.cc +++ b/src/model/dof_manager_default_inline_impl.cc @@ -1,100 +1,105 @@ /** * @file dof_manager_default_inline_impl.cc * * @author Nicolas Richart * * @date Wed Aug 12 10:57:47 2015 * * @brief Implementation of the DOFManagerDefault inline functions * * @section LICENSE * * Copyright (©) 2010-2011 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_default.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_DEFAULT_INLINE_IMPL_CC__ #define __AKANTU_DOF_MANAGER_DEFAULT_INLINE_IMPL_CC__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline bool DOFManagerDefault::isLocalOrMasterDOF(UInt dof_num) { Int dof_type = this->dofs_type(dof_num); return (dof_type == Int(_nt_normal)) || (dof_type == Int(_nt_master)); } /* -------------------------------------------------------------------------- */ inline bool DOFManagerDefault::isSlaveDOF(UInt dof_num) { Int dof_type = this->dofs_type(dof_num); return (dof_type >= 0); } /* -------------------------------------------------------------------------- */ inline const Array & DOFManagerDefault::getLocalEquationNumbers(const ID & dof_id) const { return this->getDOFData(dof_id).local_equation_number; } inline const Array & DOFManagerDefault::getDOFsAssociatedNodes(const ID & dof_id) const { - const DOFDataDefault & dof_data = this->getDOFDataTyped(dof_id); + const DOFDataDefault & dof_data = + this->getDOFDataTyped(dof_id); return dof_data.associated_nodes; } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & element_equation_number) { for (UInt i = 0, ld = 0; i < connectivity.size(); ++i) { UInt n = connectivity(i); for (UInt d = 0; d < nb_degree_of_freedom; ++d, ++ld) { element_equation_number(ld) = equation_numbers(n * nb_degree_of_freedom + d); } } } /* -------------------------------------------------------------------------- */ inline UInt DOFManagerDefault::localToGlobalEquationNumber(UInt local) const { return this->global_equation_number(local); } +/* -------------------------------------------------------------------------- */ +inline bool DOFManagerDefault::hasGlobalEquationNumber(UInt global) const { + auto it = this->global_to_local_mapping.find(global); + return (it != this->global_to_local_mapping.end()); +} + /* -------------------------------------------------------------------------- */ inline UInt DOFManagerDefault::globalToLocalEquationNumber(UInt global) const { - equation_numbers_map::const_iterator it = - this->global_to_local_mapping.find(global); + auto it = this->global_to_local_mapping.find(global); AKANTU_DEBUG_ASSERT(it != this->global_to_local_mapping.end(), "This global equation number " << global << " does not exists in " << this->id); return it->second; } /* -------------------------------------------------------------------------- */ inline Int DOFManagerDefault::getDOFType(UInt local_id) const { return this->dofs_type(local_id); } /* -------------------------------------------------------------------------- */ - __END_AKANTU__ #endif /* __AKANTU_DOF_MANAGER_DEFAULT_INLINE_IMPL_CC_ */ diff --git a/src/model/non_linear_solver_linear.hh b/src/model/non_linear_solver_linear.hh index 7d14507b9..09f8b2915 100644 --- a/src/model/non_linear_solver_linear.hh +++ b/src/model/non_linear_solver_linear.hh @@ -1,76 +1,78 @@ /** * @file non_linear_solver_linear.hh * * @author Nicolas Richart * * @date Tue Aug 25 00:48:07 2015 * * @brief Default implementation of NonLinearSolver, in case no external library * is there to do the job * * @section LICENSE * * Copyright (©) 2010-2011 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 "non_linear_solver.hh" #include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LINEAR_SOLVER_LINEAR_HH__ #define __AKANTU_NON_LINEAR_SOLVER_LINEAR_HH__ namespace akantu { class DOFManagerDefault; } __BEGIN_AKANTU__ class NonLinearSolverLinear : public NonLinearSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NonLinearSolverLinear(DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id = "non_linear_solver_linear", UInt memory_id = 0); virtual ~NonLinearSolverLinear(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Function that solve the non linear system described by the dof manager and /// the solver callback functions virtual void solve(SolverCallback & solver_callback); + AKANTU_GET_MACRO_NOT_CONST(Solver, solver, SparseSolverMumps &); + AKANTU_GET_MACRO(Solver, solver, const SparseSolverMumps &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: DOFManagerDefault & dof_manager; /// Sparse solver used for the linear solves SparseSolverMumps solver; }; __END_AKANTU__ #endif /* __AKANTU_NON_LINEAR_SOLVER_LINEAR_HH__ */ diff --git a/src/model/non_linear_solver_newton_raphson.hh b/src/model/non_linear_solver_newton_raphson.hh index d3e929d49..7faf40241 100644 --- a/src/model/non_linear_solver_newton_raphson.hh +++ b/src/model/non_linear_solver_newton_raphson.hh @@ -1,101 +1,104 @@ /** * @file non_linear_solver_newton_raphson.hh * * @author Nicolas Richart * * @date Tue Aug 25 00:48:07 2015 * * @brief Default implementation of NonLinearSolver, in case no external library * is there to do the job * * @section LICENSE * * Copyright (©) 2010-2011 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 "non_linear_solver.hh" #include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LINEAR_SOLVER_NEWTON_RAPHSON_HH__ #define __AKANTU_NON_LINEAR_SOLVER_NEWTON_RAPHSON_HH__ namespace akantu { class DOFManagerDefault; } __BEGIN_AKANTU__ class NonLinearSolverNewtonRaphson : public NonLinearSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NonLinearSolverNewtonRaphson(DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id = "non_linear_solver_newton_raphson", UInt memory_id = 0); virtual ~NonLinearSolverNewtonRaphson(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Function that solve the non linear system described by the dof manager and /// the solver callback functions virtual void solve(SolverCallback & solver_callback); + AKANTU_GET_MACRO_NOT_CONST(Solver, solver, SparseSolverMumps &); + AKANTU_GET_MACRO(Solver, solver, const SparseSolverMumps &); + protected: /// test the convergence compare norm of array to convergence_criteria bool testConvergence(const Array & array); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: DOFManagerDefault & dof_manager; /// Sparse solver used for the linear solves SparseSolverMumps solver; /// Type of convergence criteria SolveConvergenceCriteria convergence_criteria_type; /// convergence threshold Real convergence_criteria; /// Max number of iterations int max_iterations; /// Number of iterations at last solve call int n_iter; /// Convergence error at last solve call Real error; /// Did the last call to solve reached convergence bool converged; /// Force a re-computation of the jacobian matrix bool force_linear_recompute; }; __END_AKANTU__ #endif /* __AKANTU_NON_LINEAR_SOLVER_NEWTON_RAPHSON_HH__ */ diff --git a/src/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc index 760148ccb..526552327 100644 --- a/src/solver/sparse_matrix_aij.cc +++ b/src/solver/sparse_matrix_aij.cc @@ -1,228 +1,225 @@ /** * @file sparse_matrix_aij.cc * * @author Nicolas Richart * * @date Tue Aug 18 16:31:23 2015 * * @brief Implementation of the AIJ sparse matrix * * @section LICENSE * * Copyright (©) 2010-2011 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 "sparse_matrix_aij.hh" #include "dof_manager_default.hh" -#include "terms_to_assemble.hh" #include "dof_synchronizer.hh" +#include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(DOFManagerDefault & dof_manager, const MatrixType & matrix_type, const ID & id) : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager), irn(0, 1, id + ":irn"), jcn(0, 1, id + ":jcn"), a(0, 1, id + ":a"), profile_release(1), value_release(1) {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id) : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager), irn(matrix.irn, true, id + ":irn"), jcn(matrix.jcn, true, id + ":jcn"), a(matrix.a, true, id + ":a"), profile_release(1), value_release(1) {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::~SparseMatrixAIJ() {} /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::applyBoundary(Real block_val) { AKANTU_DEBUG_IN(); - const Array & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); + // clang-format off + const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); - Array::const_scalar_iterator irn_val = irn.begin(); - Array::const_scalar_iterator jcn_val = jcn.begin(); - Array::scalar_iterator a_val = a.begin(); + auto irn_it = irn.begin(); + auto jcn_it = jcn.begin(); + auto a_it = a.begin(); + auto a_end = a.end(); - for (UInt i = 0; i < nb_non_zero; ++i) { - UInt ni = this->dof_manager.globalToLocalEquationNumber(*irn_val - 1); - UInt nj = this->dof_manager.globalToLocalEquationNumber(*jcn_val - 1); + for (;a_it != a_end; ++a_it, ++irn_it, ++jcn_it) { + UInt ni = this->dof_manager.globalToLocalEquationNumber(*irn_it - 1); + UInt nj = this->dof_manager.globalToLocalEquationNumber(*jcn_it - 1); if (blocked_dofs(ni) || blocked_dofs(nj)) { - if (*irn_val != *jcn_val) { - *a_val = 0; - } else { - if (this->dof_manager.isLocalOrMasterDOF(ni)) { - *a_val = block_val; - } else { - *a_val = 0; - } - } + *a_it = *irn_it != *jcn_it ? 0. + : this->dof_manager.isLocalOrMasterDOF(ni) ? block_val + : 0.; } - ++irn_val; - ++jcn_val; - ++a_val; } this->value_release++; + // clang-format on AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); + UInt m = this->size; outfile << "%%MatrixMarket matrix coordinate pattern"; - if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; - - UInt m = this->size; outfile << m << " " << m << " " << this->nb_non_zero << std::endl; for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn.storage()[i] << " " << this->jcn.storage()[i] << " 1" << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); + // open and set the properties of the stream std::ofstream outfile; - outfile.precision(std::numeric_limits::digits10); - outfile.open(filename.c_str()); + outfile.precision(std::numeric_limits::digits10); + // write header outfile << "%%MatrixMarket matrix coordinate real"; - if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; - outfile << this->size << " " << this->size << " " << this->nb_non_zero << std::endl; + // write content for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn(i) << " " << this->jcn(i) << " " << this->a(i) << std::endl; } + // time to end outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const Array & x, Array & y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); y *= beta; - Array::const_scalar_iterator i_it = this->irn.storage(); - Array::const_scalar_iterator j_it = this->jcn.storage(); - Array::const_scalar_iterator a_it = this->a.storage(); - Array::const_scalar_iterator x_it = x.storage(); - Array::scalar_iterator y_it = y.storage(); + auto i_it = this->irn.begin(); + auto j_it = this->jcn.begin(); + auto a_it = this->a.begin(); + auto a_end = this->a.end(); + auto x_it = x.begin(); + auto y_it = y.begin(); - for (UInt k = 0; k < this->nb_non_zero; ++k, ++i_it, ++j_it, ++a_it) { + for (; a_it != a_end; ++i_it, ++j_it, ++a_it) { Int i = this->dof_manager.globalToLocalEquationNumber(*i_it - 1); Int j = this->dof_manager.globalToLocalEquationNumber(*j_it - 1); const Real & A = *a_it; y_it[i] += alpha * A * x_it[j]; if ((this->matrix_type == _symmetric) && (i != j)) y_it[j] += alpha * A * x_it[i]; } if (this->dof_manager.hasSynchronizer()) this->dof_manager.getSynchronizer().reduceSynchronize(y); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyContent(const SparseMatrix & matrix) { AKANTU_DEBUG_IN(); const SparseMatrixAIJ & mat = dynamic_cast(matrix); AKANTU_DEBUG_ASSERT(nb_non_zero == mat.getNbNonZero(), "The to matrix don't have the same profiles"); memcpy(a.storage(), mat.getA().storage(), nb_non_zero * sizeof(Real)); this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template +template void SparseMatrixAIJ::addMeToTemplated(MatrixType & B, Real alpha) const { - Array::const_scalar_iterator i_it = this->irn.begin(); - Array::const_scalar_iterator j_it = this->jcn.begin(); - Array::const_scalar_iterator a_it = this->a.begin(); - for (UInt n = 0; n < this->nb_non_zero; ++n, ++a_it, ++i_it, ++j_it) { + auto i_it = this->irn.begin(); + auto j_it = this->jcn.begin(); + auto a_it = this->a.begin(); + auto a_end = this->a.end(); + + for (; a_it != a_end; ++a_it, ++i_it, ++j_it) { const Int & i = *i_it; const Int & j = *j_it; const Real & A_ij = *a_it; + B.addToMatrix(i - 1, j - 1, alpha * A_ij); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::addMeTo(SparseMatrix & B, Real alpha) const { - if(SparseMatrixAIJ * B_aij = dynamic_cast(&B)) { + if (SparseMatrixAIJ * B_aij = dynamic_cast(&B)) { this->addMeToTemplated(*B_aij, alpha); } else { // this->addMeToTemplated(*this, alpha); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::mul(Real alpha) { this->a *= alpha; this->value_release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::clear() { a.set(0.); this->value_release++; } __END_AKANTU__ diff --git a/src/solver/terms_to_assemble.hh b/src/solver/terms_to_assemble.hh index 2c6662781..61e968029 100644 --- a/src/solver/terms_to_assemble.hh +++ b/src/solver/terms_to_assemble.hh @@ -1,90 +1,99 @@ /** * @file terms_to_assemble.hh * * @author Nicolas Richart * * @date creation Tue Dec 20 2016 * * @brief List of terms to assemble to a matrix * * @section LICENSE * * Copyright (©) 2010-2011 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_common.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TERMS_TO_ASSEMBLE_HH__ #define __AKANTU_TERMS_TO_ASSEMBLE_HH__ namespace akantu { class TermsToAssemble { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: TermsToAssemble() {} virtual ~TermsToAssemble() {} class TermToAssemble { public: - TermToAssemble(UInt i, UInt j) : _i(i), _j(j) {} - inline void operator+=(Real val) { this->val = val; } - inline operator Real() const { return val; } + TermToAssemble(UInt i, UInt j) : _i(i), _j(j), val(0.) {} + inline TermToAssemble & operator=(Real val) { + this->val = val; + return *this; + } + inline TermToAssemble operator+=(Real val) { + this->val += val; + return *this; + } + inline operator Real() const { return val; } inline UInt i() const { return _i; } inline UInt j() const { return _j; } + private: UInt _i, _j; Real val; }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline TermToAssemble & operator()(UInt i, UInt j) { terms.emplace_back(i, j); return terms.back(); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ private: typedef std::vector TermsContainer; + public: typedef TermsContainer::const_iterator const_terms_iterator; const_terms_iterator begin() const { return terms.begin(); } const_terms_iterator end() const { return terms.end(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: TermsContainer terms; }; } // akantu #endif /* __AKANTU_TERMS_TO_ASSEMBLE_HH__ */ diff --git a/test/test_model/test_model_solver/test_model_solver.cc b/test/test_model/test_model_solver/test_model_solver.cc index 6e2c5c590..325e996a4 100644 --- a/test/test_model/test_model_solver/test_model_solver.cc +++ b/test/test_model/test_model_solver/test_model_solver.cc @@ -1,141 +1,164 @@ /** * @file test_dof_manager_default.cc * * @author Nicolas Richart * * @date Wed Feb 24 12:28:44 2016 * * @brief Test default dof manager * * @section LICENSE * * Copyright (©) 2010-2011 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_random_generator.hh" #include "dof_manager.hh" #include "dof_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "model_solver.hh" +#include "non_linear_solver_newton_raphson.hh" #include "sparse_matrix.hh" +#include "sparse_solver_mumps.hh" + /* -------------------------------------------------------------------------- */ #include "test_model_solver_my_model.hh" /* -------------------------------------------------------------------------- */ -#include #include /* -------------------------------------------------------------------------- */ using namespace akantu; static void genMesh(Mesh & mesh, UInt nb_nodes); static void printResults(MyModel & model, UInt nb_nodes); Real F = -10; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); UInt prank = StaticCommunicator::getStaticCommunicator().whoAmI(); std::cout << std::setprecision(7); - UInt global_nb_nodes = 11; + UInt global_nb_nodes = 100; Mesh mesh(1); - if (prank == 0) + RandGenerator::seed(1); + + if (prank == 0) { genMesh(mesh, global_nb_nodes); + } + + std::cout << prank << RandGenerator::seed() << std::endl; mesh.distribute(); MyModel model(F, mesh, false); + UInt node = 0; + for (auto pos : mesh.getNodes()) { + std::cout << prank << " " << node << " pos: " << pos << " [" + << mesh.getNodeGlobalId(node) << "] " << mesh.getNodeType(node) + << " -- " << model.forces(node) << " " << model.blocked(node) + << std::endl; + ++node; + } + model.getNewSolver("static", _tsst_static, _nls_newton_raphson); model.setIntegrationScheme("static", "disp", _ist_pseudo_time); + NonLinearSolver & solver = + model.getDOFManager().getNonLinearSolver("static"); + solver.set("max_iterations", 2); + + model.solveStep(); + dynamic_cast( + model.getDOFManager().getNonLinearSolver("static")) + .getSolver() + .analysis(); + printResults(model, global_nb_nodes); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); } conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } } /* -------------------------------------------------------------------------- */ void printResults(MyModel & model, UInt nb_nodes) { UInt prank = StaticCommunicator::getStaticCommunicator().whoAmI(); auto & sync = dynamic_cast(model.getDOFManager()) .getSynchronizer(); if (prank == 0) { Array global_displacement(nb_nodes); Array global_forces(nb_nodes); Array global_blocked(nb_nodes); sync.gather(model.forces, global_forces); - - auto force_it = global_forces.begin(); - bool a = std::abs(global_forces(nb_nodes - 1) - F) > 1e-10; - while(a) {} - sync.gather(model.displacement, global_displacement); sync.gather(model.blocked, global_blocked); + auto force_it = global_forces.begin(); auto disp_it = global_displacement.begin(); auto blocked_it = global_blocked.begin(); std::cout << "node" << ", " << std::setw(8) << "disp" << ", " << std::setw(8) << "force" << ", " << std::setw(8) << "blocked" << std::endl; UInt node = 0; for (; disp_it != global_displacement.end(); ++disp_it, ++force_it, ++blocked_it, ++node) { std::cout << node << ", " << std::setw(8) << *disp_it << ", " << std::setw(8) << *force_it << ", " << std::setw(8) << *blocked_it << std::endl; std::cout << std::flush; } } else { sync.gather(model.forces); sync.gather(model.displacement); sync.gather(model.blocked); } } diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic.cc b/test/test_model/test_model_solver/test_model_solver_dynamic.cc index aa78e5390..2c046d2d7 100644 --- a/test/test_model/test_model_solver/test_model_solver_dynamic.cc +++ b/test/test_model/test_model_solver/test_model_solver_dynamic.cc @@ -1,171 +1,171 @@ /** * @file test_dof_manager_default.cc * * @author Nicolas Richart * * @date Wed Feb 24 12:28:44 2016 * * @brief Test default dof manager * * @section LICENSE * * Copyright (©) 2010-2011 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 "data_accessor.hh" #include "dof_manager.hh" #include "dof_manager_default.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "model_solver.hh" #include "non_linear_solver.hh" #include "sparse_matrix.hh" #include "static_communicator.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ #include "test_model_solver_my_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef EXPLICIT #define EXPLICIT true #endif using namespace akantu; static void genMesh(Mesh & mesh, UInt nb_nodes); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); UInt prank = StaticCommunicator::getStaticCommunicator().whoAmI(); - UInt global_nb_nodes = 21; + UInt global_nb_nodes = 201; UInt max_steps = 200; Real time_step = 0.001; Mesh mesh(1); Real F = -9.81; bool _explicit = EXPLICIT; if (prank == 0) genMesh(mesh, global_nb_nodes); mesh.distribute(); MyModel model(F, mesh, _explicit); if (!_explicit) { model.getNewSolver("dynamic", _tsst_dynamic, _nls_newton_raphson); model.setIntegrationScheme("dynamic", "disp", _ist_trapezoidal_rule_2, IntegrationScheme::_displacement); } else { model.getNewSolver("dynamic", _tsst_dynamic_lumped, _nls_lumped); model.setIntegrationScheme("dynamic", "disp", _ist_central_difference, IntegrationScheme::_acceleration); } model.setTimeStep(time_step); // #if EXPLICIT == true // std::ofstream output("output_dynamic_explicit.csv"); // #else // std::ofstream output("output_dynamic_implicit.csv"); // #endif if (prank == 0) { std::cout << std::scientific; std::cout << std::setw(14) << "time" << "," << std::setw(14) << "wext" << "," << std::setw(14) << "epot" << "," << std::setw(14) << "ekin" << "," << std::setw(14) << "total" << std::endl; } Real wext = 0.; model.getDOFManager().clearResidual(); model.assembleResidual(); Real epot = model.getPotentialEnergy(); Real ekin = model.getKineticEnergy(); Real einit = ekin + epot; Real etot = ekin + epot - wext - einit; if (prank == 0) { std::cout << std::setw(14) << 0. << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << std::endl; } #if EXPLICIT == false NonLinearSolver & solver = model.getDOFManager().getNonLinearSolver("dynamic"); - solver.set("max_iterations", 2); + solver.set("max_iterations", 20); #endif for (UInt i = 1; i < max_steps + 1; ++i) { model.solveStep(); #if EXPLICIT == false - int nb_iter = solver.get("nb_iterations"); - Real error = solver.get("error"); - bool converged = solver.get("converged"); - if (prank == 0) - std::cerr << error << " " << nb_iter << " -> " << converged << std::endl; + //int nb_iter = solver.get("nb_iterations"); + //Real error = solver.get("error"); + //bool converged = solver.get("converged"); + //if (prank == 0) + // std::cerr << error << " " << nb_iter << " -> " << converged << std::endl; #endif epot = model.getPotentialEnergy(); ekin = model.getKineticEnergy(); wext += model.getExternalWorkIncrement(); etot = ekin + epot - wext - einit; if (prank == 0) { std::cout << std::setw(14) << time_step * i << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << std::endl; } - if (std::abs(etot) > 1e1) { + if (std::abs(etot) > 1e12) { AKANTU_DEBUG_ERROR("The total energy of the system is not conserved!"); } } // output.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); } conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } } diff --git a/test/test_model/test_model_solver/test_model_solver_my_model.hh b/test/test_model/test_model_solver/test_model_solver_my_model.hh index 75ac8385f..2c1e95045 100644 --- a/test/test_model/test_model_solver/test_model_solver_my_model.hh +++ b/test/test_model/test_model_solver/test_model_solver_my_model.hh @@ -1,443 +1,443 @@ /** * @file test_dof_manager_default.cc * * @author Nicolas Richart * * @date Wed Feb 24 12:28:44 2016 * * @brief Test default dof manager * * @section LICENSE * * Copyright (©) 2010-2011 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 "data_accessor.hh" #include "dof_manager_default.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "model_solver.hh" #include "sparse_matrix.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { #ifndef __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ #define __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ /** * =\o-----o-----o-> F * | | * |---- L ----| */ class MyModel : public ModelSolver, public DataAccessor { public: MyModel(Real F, Mesh & mesh, bool lumped) - : ModelSolver(mesh, "model_solver", 0), mesh(mesh), + : ModelSolver(mesh, "model_solver", 0), nb_dofs(mesh.getNbNodes()), nb_elements(mesh.getNbElement()), E(1.), - A(1.), rho(1.), lumped(lumped), displacement(nb_dofs, 1, "disp"), + A(1.), rho(1.), lumped(lumped), mesh(mesh), displacement(nb_dofs, 1, "disp"), velocity(nb_dofs, 1, "velo"), acceleration(nb_dofs, 1, "accel"), blocked(nb_dofs, 1, "blocked"), forces(nb_dofs, 1, "force_ext"), internal_forces(nb_dofs, 1, "force_int"), stresses(nb_elements, 1, "stress"), strains(nb_elements, 1, "strain"), initial_lengths(nb_elements, 1, "L0") { this->initDOFManager(); this->getDOFManager().registerDOFs("disp", displacement, _dst_nodal); this->getDOFManager().registerDOFsDerivative("disp", 1, velocity); this->getDOFManager().registerDOFsDerivative("disp", 2, acceleration); this->getDOFManager().registerBlockedDOFs("disp", blocked); this->getDOFManager().getNewMatrix("K", _symmetric); this->getDOFManager().getNewMatrix("M", "K"); this->getDOFManager().getNewMatrix("J", "K"); this->getDOFManager().getNewLumpedMatrix("M"); if (lumped) { this->assembleLumpedMass(); } else { this->assembleMass(); this->assembleStiffness(); } this->assembleJacobian(); displacement.set(0.); velocity.set(0.); acceleration.set(0.); forces.set(0.); blocked.set(false); UInt global_nb_nodes = mesh.getNbGlobalNodes(); for (UInt n = 0; n < nb_dofs; ++n) { - if (mesh.getGlobalNodesIds()(n) == (global_nb_nodes - 1)) + if (mesh.getNodeGlobalId(n) == (global_nb_nodes - 1)) forces(n, _x) = F; if (mesh.getGlobalNodesIds()(n) == 0) blocked(n, _x) = true; } auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++L_it) { const Vector & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); *L_it = std::abs(p2 - p1); } this->registerDataAccessor(*this); this->registerSynchronizer( const_cast(this->mesh.getElementSynchronizer()), _gst_user_1); } void assembleLumpedMass() { Array & M = this->getDOFManager().getLumpedMatrix("M"); M.clear(); this->assembleLumpedMass(_not_ghost); if (this->mesh.getNbElement(_segment_2, _ghost) > 0) this->assembleLumpedMass(_ghost); } void assembleLumpedMass(const GhostType & ghost_type) { Array & M = this->getDOFManager().getLumpedMatrix("M"); Array m_all_el(this->mesh.getNbElement(_segment_2, ghost_type), 2); Array::vector_iterator m_it = m_all_el.begin(2); auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2); auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2); for (; cit != cend; ++cit, ++m_it) { const Vector & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); Real M_n = rho * A * L / 2; (*m_it)(0) = (*m_it)(1) = M_n; } this->getDOFManager().assembleElementalArrayLocalArray( m_all_el, M, _segment_2, ghost_type); } void assembleMass() { SparseMatrix & M = this->getDOFManager().getMatrix("M"); M.clear(); Array m_all_el(this->nb_elements, 4); Array::matrix_iterator m_it = m_all_el.begin(2, 2); auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); Matrix m(2, 2); m(0, 0) = m(1, 1) = 2; m(0, 1) = m(1, 0) = 1; // under integrated // m(0, 0) = m(1, 1) = 3./2.; // m(0, 1) = m(1, 0) = 3./2.; // lumping the mass matrix // m(0, 0) += m(0, 1); // m(1, 1) += m(1, 0); // m(0, 1) = m(1, 0) = 0; for (; cit != cend; ++cit, ++m_it) { const Vector & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); Matrix & m_el = *m_it; m_el = m; m_el *= rho * A * L / 6.; } this->getDOFManager().assembleElementalMatricesToMatrix( "M", "disp", m_all_el, _segment_2); } void assembleJacobian() {} void assembleStiffness() { SparseMatrix & K = this->getDOFManager().getMatrix("K"); K.clear(); Matrix k(2, 2); k(0, 0) = k(1, 1) = 1; k(0, 1) = k(1, 0) = -1; Array k_all_el(this->nb_elements, 4); auto k_it = k_all_el.begin(2, 2); auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); for (; cit != cend; ++cit, ++k_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); auto & k_el = *k_it; k_el = k; k_el *= E * A / L; } this->getDOFManager().assembleElementalMatricesToMatrix( "K", "disp", k_all_el, _segment_2); } void assembleResidual() { this->getDOFManager().assembleToResidual("disp", forces); internal_forces.clear(); this->assembleResidual(_not_ghost); this->synchronize(_gst_user_1); this->getDOFManager().assembleToResidual("disp", internal_forces, -1.); - auto & comm = StaticCommunicator::getStaticCommunicator(); - const auto & dof_manager_default = - dynamic_cast(this->getDOFManager()); - const auto & residual = dof_manager_default.getResidual(); - int prank = comm.whoAmI(); - int psize = comm.getNbProc(); - - for (int p = 0; p < psize; ++p) { - if (prank == p) { - UInt local_dof = 0; - for (auto res : residual) { - UInt global_dof = dof_manager_default.localToGlobalEquationNumber(local_dof); - std::cout << local_dof << " [" << global_dof << " - " - << dof_manager_default.getDOFType(local_dof) << "]: " << res - << std::endl; - ++local_dof; - } - std::cout << std::flush; - } - comm.barrier(); - } - - comm.barrier(); - if(prank == 0) std::cout << "===========================" << std::endl; + // auto & comm = StaticCommunicator::getStaticCommunicator(); + // const auto & dof_manager_default = + // dynamic_cast(this->getDOFManager()); + // const auto & residual = dof_manager_default.getResidual(); + // int prank = comm.whoAmI(); + // int psize = comm.getNbProc(); + + // for (int p = 0; p < psize; ++p) { + // if (prank == p) { + // UInt local_dof = 0; + // for (auto res : residual) { + // UInt global_dof = dof_manager_default.localToGlobalEquationNumber(local_dof); + // std::cout << local_dof << " [" << global_dof << " - " + // << dof_manager_default.getDOFType(local_dof) << "]: " << res + // << std::endl; + // ++local_dof; + // } + // std::cout << std::flush; + // } + // comm.barrier(); + // } + + // comm.barrier(); + // if(prank == 0) std::cout << "===========================" << std::endl; } void assembleResidual(const GhostType & ghost_type) { Array forces_internal_el( this->mesh.getNbElement(_segment_2, ghost_type), 2); auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2); auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2); auto f_it = forces_internal_el.begin(2); auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++f_it, ++strain_it, ++stress_it, ++L_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real u1 = this->displacement(n1, _x); Real u2 = this->displacement(n2, _x); *strain_it = (u2 - u1) / *L_it; *stress_it = E * *strain_it; Real f_n = A * *stress_it; Vector & f = *f_it; f(0) = -f_n; f(1) = f_n; } this->getDOFManager().assembleElementalArrayLocalArray( forces_internal_el, internal_forces, _segment_2, ghost_type); } Real getPotentialEnergy() { Real res = 0; if (!lumped) { res = this->mulVectMatVect(this->displacement, this->getDOFManager().getMatrix("K"), this->displacement); } else { auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto strain_end = this->strains.end(); auto L_it = this->initial_lengths.begin(); for (; strain_it != strain_end; ++strain_it, ++stress_it, ++L_it) { res += *strain_it * *stress_it * A * *L_it; } StaticCommunicator::getStaticCommunicator().allReduce(res, _so_sum); } return res / 2.; } Real getKineticEnergy() { Real res = 0; if (!lumped) { res = this->mulVectMatVect( this->velocity, this->getDOFManager().getMatrix("M"), this->velocity); } else { auto & m = this->getDOFManager().getLumpedMatrix("M"); auto it = velocity.begin(); auto end = velocity.end(); auto m_it = m.begin(); for (UInt node = 0; it != end; ++it, ++m_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += *m_it * *it * *it; } StaticCommunicator::getStaticCommunicator().allReduce(res, _so_sum); } return res / 2.; } Real getExternalWorkIncrement() { Real res = 0; auto it = velocity.begin(); auto end = velocity.end(); auto if_it = internal_forces.begin(); auto ef_it = forces.begin(); auto b_it = blocked.begin(); for (UInt node = 0; it != end; ++it, ++if_it, ++ef_it, ++b_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += (*b_it ? -*if_it : *ef_it) * *it; } StaticCommunicator::getStaticCommunicator().allReduce(res, _so_sum); return res * this->getTimeStep(); } Real mulVectMatVect(const Array & x, const SparseMatrix & A, const Array & y) { Array Ay(this->nb_dofs, 1, 0.); A.matVecMul(y, Ay); Real res = 0.; Array::const_scalar_iterator it = Ay.begin(); Array::const_scalar_iterator end = Ay.end(); Array::const_scalar_iterator x_it = x.begin(); for (UInt node = 0; it != end; ++it, ++x_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += *x_it * *it; } StaticCommunicator::getStaticCommunicator().allReduce(res, _so_sum); return res; } void predictor() {} void corrector() {} /* ------------------------------------------------------------------------ */ UInt getNbData(const Array & elements, const SynchronizationTag &) const { return elements.getSize() * sizeof(Real); } void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag == _gst_user_1) { for (const auto & el : elements) { buffer << this->stresses(el.element); } } } void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag == _gst_user_1) { auto cit = this->mesh.getConnectivity(_segment_2, _ghost).begin(2); for (const auto & el : elements) { Real stress; buffer >> stress; Real f = A * stress; Vector conn = cit[el.element]; this->internal_forces(conn(0), _x) += -f; this->internal_forces(conn(1), _x) += f; } } } private: - Mesh & mesh; UInt nb_dofs; UInt nb_elements; Real E, A, rho; bool lumped; public: + Mesh & mesh; Array displacement; Array velocity; Array acceleration; Array blocked; Array forces; Array internal_forces; Array stresses; Array strains; Array initial_lengths; }; #endif /* __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ */ } // akantu diff --git a/test/test_solver/CMakeLists.txt b/test/test_solver/CMakeLists.txt index 9d465f19b..507452dc6 100644 --- a/test/test_solver/CMakeLists.txt +++ b/test/test_solver/CMakeLists.txt @@ -1,96 +1,97 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Fri Sep 03 2010 # @date last modification: Wed Dec 16 2015 # # @brief configuration for solver tests # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 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 . # # @section DESCRIPTION # #=============================================================================== add_mesh(test_solver_mesh triangle.geo 2 1) add_mesh(test_matrix_mesh square.geo 2 1) add_mesh(test_solver_petsc_mesh 1D_bar.geo 1 1) register_test(test_sparse_matrix_profile SOURCES test_sparse_matrix_profile.cc DEPENDS test_solver_mesh PACKAGE implicit ) register_test(test_sparse_matrix_assemble SOURCES test_sparse_matrix_assemble.cc DEPENDS test_solver_mesh PACKAGE implicit ) register_test(test_sparse_matrix_product SOURCES test_sparse_matrix_product.cc FILES_TO_COPY bar.msh PACKAGE implicit ) register_test(test_sparse_solver_mumps SOURCES test_sparse_solver_mumps.cc PACKAGE mumps + PARALLEL ) register_test(test_petsc_matrix_profile SOURCES test_petsc_matrix_profile.cc DEPENDS test_matrix_mesh PACKAGE petsc ) register_test(test_petsc_matrix_profile_parallel SOURCES test_petsc_matrix_profile_parallel.cc DEPENDS test_matrix_mesh PACKAGE petsc ) register_test(test_petsc_matrix_diagonal SOURCES test_petsc_matrix_diagonal.cc DEPENDS test_solver_mesh PACKAGE petsc ) register_test(test_petsc_matrix_apply_boundary SOURCES test_petsc_matrix_apply_boundary.cc DEPENDS test_solver_mesh PACKAGE petsc ) register_test(test_solver_petsc SOURCES test_solver_petsc.cc DEPENDS test_solver_petsc_mesh PACKAGE petsc ) register_test(test_solver_petsc_parallel SOURCES test_solver_petsc.cc DEPENDS test_solver_petsc_mesh PACKAGE petsc ) diff --git a/test/test_solver/test_sparse_solver_mumps.cc b/test/test_solver/test_sparse_solver_mumps.cc index c36506e9a..b6cc0cba6 100644 --- a/test/test_solver/test_sparse_solver_mumps.cc +++ b/test/test_solver/test_sparse_solver_mumps.cc @@ -1,142 +1,159 @@ /** * @file test_sparse_matrix_product.cc * * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Sun Oct 19 2014 * * @brief test the matrix vector product in parallel * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_common.hh" #include "dof_synchronizer.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_partition_scotch.hh" #include "sparse_matrix_aij.hh" #include "sparse_solver_mumps.hh" +#include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); const UInt spatial_dimension = 1; - const UInt nb_global_dof = 100; + const UInt nb_global_dof = 11; StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); // Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); Mesh mesh(spatial_dimension); - if (prank == 0) + if (prank == 0) { genMesh(mesh, nb_global_dof); + RandGenerator::seed(1496137735); + } else { + RandGenerator::seed(2992275470); + } mesh.distribute(); - + UInt node = 0; + for (auto pos : mesh.getNodes()) { + std::cout << prank << " " << node << " pos: " << pos << " [" + << mesh.getNodeGlobalId(node) << "] " << mesh.getNodeType(node) + << std::endl; + ++node; + } UInt nb_nodes = mesh.getNbNodes(); DOFManagerDefault dof_manager(mesh, "test_dof_manager"); Array x(nb_nodes); dof_manager.registerDOFs("x", x, _dst_nodal); - Array local_equation_number(nb_nodes); - dof_manager.getEquationsNumbers("x", local_equation_number); + Array global_equation_number(nb_nodes); + dof_manager.getEquationsNumbers("x", global_equation_number); auto & A = dof_manager.getNewMatrix("A", _symmetric); Array b(nb_nodes); + TermsToAssemble terms; for (UInt i = 0; i < nb_nodes; ++i) { if (dof_manager.isLocalOrMasterDOF(i)) { - UInt eqn = local_equation_number(i); - A.addToMatrix(eqn, eqn, 1. / (eqn+1)); + UInt gi = global_equation_number(i); + terms(i, i) = 1. / (1. + gi); } } + dof_manager.assemblePreassembledMatrix("x", "x", "A", terms); + std::stringstream str; str << "Matrix_" << prank << ".mtx"; A.saveMatrix(str.str()); for (UInt n = 0; n < nb_nodes; ++n) { b(n) = 1.; } SparseSolverMumps solver(dof_manager, "A"); solver.solve(x, b); auto & sync = dynamic_cast(dof_manager).getSynchronizer(); if (prank == 0) { Array x_gathered(dof_manager.getSystemSize()); sync.gather(x, x_gathered); debug::setDebugLevel(dblTest); std::cout << x_gathered << std::endl; debug::setDebugLevel(dblWarning); UInt d = 1.; - for(auto x : x_gathered) { - if ( std::abs(x - d)/d > 1e-15 ) AKANTU_EXCEPTION("Error in the solution"); + for (auto x : x_gathered) { + if (std::abs(x - d) / d > 1e-15) + AKANTU_EXCEPTION("Error in the solution: " << x << " != " << d << " [" + << (std::abs(x - d) / d) + << "]."); ++d; } } else { sync.gather(x); } finalize(); return 0; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); } conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } }