diff --git a/packages/mumps.cmake b/packages/mumps.cmake index 229d5ff40..2db2e8225 100644 --- a/packages/mumps.cmake +++ b/packages/mumps.cmake @@ -1,112 +1,112 @@ #=============================================================================== # @file mumps.cmake # # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for mumps support # # @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 . # #=============================================================================== package_declare(Mumps EXTERNAL DESCRIPTION "Add Mumps support in akantu" SYSTEM ON third-party/cmake/mumps.cmake ) package_declare_sources(Mumps - solver/solver_mumps.cc - solver/solver_mumps.hh + solver/sparse_solver_mumps.cc + solver/sparse_solver_mumps.hh ) set(_mumps_float_type ${AKANTU_FLOAT_TYPE}) if(AKANTU_FLOAT_TYPE STREQUAL "float" OR AKANTU_FLOAT_TYPE STREQUAL "double") set(_mumps_components ${AKANTU_FLOAT_TYPE}) else() if(DEFINED AKANTU_FLOAT_TYPE) message(FATAL_ERROR "MUMPS doea not support floating point type \"${AKANTU_FLOAT_TYPE}\"") endif() endif() package_get_option_name(parallel _par_option) if(${_par_option}) package_set_find_package_extra_options(Mumps ARGS COMPONENTS "parallel" ${_mumps_components}) package_add_third_party_script_variable(Mumps MUMPS_TYPE "par") package_set_package_system_dependency(Mumps deb libmumps) package_set_package_system_dependency(Mumps deb-src libmumps-dev) else() package_set_find_package_extra_options(Mumps ARGS COMPONENTS "sequential" ${_mumps_components}) package_add_third_party_script_variable(Mumps MUMPS_TYPE "seq") package_set_package_system_dependency(Mumps deb libmumps-seq) package_set_package_system_dependency(Mumps deb-src libmumps-seq-dev) endif() package_use_system(Mumps _use_system) if(NOT _use_system) enable_language(Fortran) set(AKANTU_USE_MUMPS_VERSION "4.10.0" CACHE STRING "Default Mumps version to compile") mark_as_advanced(AKANTU_USE_MUMPS_VERSION) set_property(CACHE AKANTU_USE_MUMPS_VERSION PROPERTY STRINGS "4.9.2" "4.10.0" "5.0.0") package_get_option_name(MPI _mpi_option) if(${_mpi_option}) package_add_dependencies(Mumps ScaLAPACK MPI) endif() package_add_dependencies(Mumps Scotch BLAS) endif() package_declare_documentation(Mumps "This package enables the \\href{http://mumps.enseeiht.fr/}{MUMPS} parallel direct solver for sparce matrices." "This is necessary to solve static or implicit problems." "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "" "\\begin{command}" " > sudo apt-get install libmumps-seq-dev # for sequential" " > sudo apt-get install libmumps-dev # for parallel" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{command}" " > sudo port install mumps" "\\end{command}" "" "If you activate the advanced option AKANTU\\_USE\\_THIRD\\_PARTY\\_MUMPS the make system of akantu can automatically compile MUMPS. For this you will have to download MUMPS from \\url{http://mumps.enseeiht.fr/} or \\url{http://graal.ens-lyon.fr/MUMPS} and place it in \\shellcode{/third-party}" ) package_declare_extra_files_to_package(MUMPS PROJECT third-party/MUMPS_4.10.0_make.inc.cmake third-party/MUMPS_5.0.0.patch third-party/MUMPS_4.10.0.patch third-party/MUMPS_4.9.2_make.inc.cmake third-party/cmake/mumps.cmake cmake/Modules/FindMumps.cmake ) diff --git a/src/model/dof_manager.hh b/src/model/dof_manager.hh index dec0b599a..7b89eee14 100644 --- a/src/model/dof_manager.hh +++ b/src/model/dof_manager.hh @@ -1,464 +1,467 @@ /** * @file dof_manager.hh * * @author Nicolas Richart * * @date Wed Jul 22 11:43:43 2015 * * @brief Class handling the different types of dofs * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_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 { 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); virtual ~DOFManager(); /* ------------------------------------------------------------------------ */ /* 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 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 void registerSparseMatrix(const ID & matrix_id, SparseMatrix & matrix); /// register a non linear solver instantiated by a derived class void registerNonLinearSolver(const ID & non_linear_solver_id, NonLinearSolver & non_linear_solver); /// register a time step solver instantiated by a derived class void registerTimeStepSolver(const ID & time_step_solver_id, TimeStepSolver & time_step_solver); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the equation numbers corresponding to a dof_id. This might be used to /// access the matrix. virtual void getEquationsNumbers(const ID & dof_id, Array & equation_numbers); /// 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; /* ------------------------------------------------------------------------ */ /* 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 StaticCommunicator &); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ private: /// fills the nodes_to_elements structure void fillNodesToElements(); public: /// function to implement to react on akantu::NewNodesEvent virtual void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event); /// function to implement to react on akantu::RemovedNodesEvent virtual void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event); /// function to implement to react on akantu::NewElementsEvent virtual void onElementsAdded(const Array & elements_list, const NewElementsEvent & event); /// function to implement to react on akantu::RemovedElementsEvent virtual void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /// function to implement to react on akantu::ChangedElementsEvent virtual void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event); 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; }; typedef Array *> NodesToElements; /// This info is stored to simplify the dynamic changes NodesToElements nodes_to_elements; /// type to store dofs information typedef std::map DOFStorage; /// type to store all the matrices typedef std::map SparseMatricesMap; /// type to store all the lumped matrices typedef std::map *> LumpedMatricesMap; /// type to store all the non linear solver typedef std::map NonLinearSolversMap; /// type to store all the time step solver typedef std::map TimeStepSolversMap; /// 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; /// Total number of degrees of freedom (size with the ghosts) UInt local_system_size; /// Number of purely local dofs (size without the ghosts) UInt pure_local_system_size; /// Total number of degrees of freedom UInt system_size; /// Communicator used for this manager, should be the same as in the mesh if a /// mesh is registered const StaticCommunicator & communicator; }; } // akantu #include "dof_manager_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_HH__ */ diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 703383a18..1012c39e6 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,865 +1,878 @@ /** * @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; + 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) + 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; + 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); - this->localToGlobalEquationNumber(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); } } /* -------------------------------------------------------------------------- */ 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); - this->localToGlobalEquationNumber(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() { +void DOFManagerDefault::applyBoundary(const ID & matrix_id) { this->updateGlobalBlockedDofs(); - SparseMatrixAIJ & J = this->getMatrix("J"); + 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() { 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 aaa9eb308..dc5ac4eba 100644 --- a/src/model/dof_manager_default.hh +++ b/src/model/dof_manager_default.hh @@ -1,333 +1,333 @@ /** * @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(); + 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 + /// Get the global residual array across processors (SMP call) const Array & getGlobalResidual(); /// Get the residual array - const Array & getResidual(); + 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; /// return the local index of the global equation number inline UInt globalToLocalEquationNumber(UInt global) const; /// converts local equation numbers to global equation numbers; - template inline void localToGlobalEquationNumber(S & inout); + 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 e60884127..6216b6d99 100644 --- a/src/model/dof_manager_default_inline_impl.cc +++ b/src/model/dof_manager_default_inline_impl.cc @@ -1,104 +1,100 @@ /** * @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); 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); } } } /* -------------------------------------------------------------------------- */ -template -inline void DOFManagerDefault::localToGlobalEquationNumber(S & inout) { - for (UInt i = 0; i < inout.size(); ++i) { - inout(i) = this->global_equation_number(inout(i)); - } -} - -/* -------------------------------------------------------------------------- */ -template<> -inline void DOFManagerDefault::localToGlobalEquationNumber(UInt & inout) { - inout = this->global_equation_number(inout); +inline UInt DOFManagerDefault::localToGlobalEquationNumber(UInt local) const { + return this->global_equation_number(local); } /* -------------------------------------------------------------------------- */ inline UInt DOFManagerDefault::globalToLocalEquationNumber(UInt global) const { equation_numbers_map::const_iterator 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__ */ +#endif /* __AKANTU_DOF_MANAGER_DEFAULT_INLINE_IMPL_CC_ */ diff --git a/src/model/dof_manager_inline_impl.cc b/src/model/dof_manager_inline_impl.cc index 7fa134ceb..7d963b7c2 100644 --- a/src/model/dof_manager_inline_impl.cc +++ b/src/model/dof_manager_inline_impl.cc @@ -1,143 +1,148 @@ /** * @file dof_manager_inline_impl.cc * * @author Nicolas Richart * * @date Wed Aug 12 11:07:01 2015 * * @brief * * @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_INLINE_IMPL_CC__ #define __AKANTU_DOF_MANAGER_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasDOFs(const ID & dof_id) const { DOFStorage::const_iterator it = this->dofs.find(dof_id); return it != this->dofs.end(); } /* -------------------------------------------------------------------------- */ inline DOFManager::DOFData & DOFManager::getDOFData(const ID & dof_id) { DOFStorage::iterator 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 { DOFStorage::const_iterator 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; } /* -------------------------------------------------------------------------- */ 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 != NULL); } /* -------------------------------------------------------------------------- */ 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 != NULL); } /* -------------------------------------------------------------------------- */ 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] == NULL)) 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] != NULL)); } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getSolution(const ID & dofs_id) const { 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__ */ diff --git a/src/model/integration_scheme/generalized_trapezoidal.cc b/src/model/integration_scheme/generalized_trapezoidal.cc index e75fb6df0..cef1e7214 100644 --- a/src/model/integration_scheme/generalized_trapezoidal.cc +++ b/src/model/integration_scheme/generalized_trapezoidal.cc @@ -1,179 +1,192 @@ /** * @file generalized_trapezoidal.cc * * @author Nicolas Richart * * @date creation: Mon Jul 04 2011 * @date last modification: Thu Jun 05 2014 * * @brief implementation of inline functions * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "generalized_trapezoidal.hh" #include "mesh.hh" #include "aka_array.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ GeneralizedTrapezoidal::GeneralizedTrapezoidal(DOFManager & dof_manager, const ID & dof_id, Real alpha) : IntegrationScheme1stOrder(dof_manager, dof_id), alpha(alpha) { this->registerParam("alpha", this->alpha, alpha, _pat_parsmod, "The alpha parameter"); } /* -------------------------------------------------------------------------- */ void GeneralizedTrapezoidal::predictor(Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.getSize(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt d = 0; d < nb_degree_of_freedom; d++) { if (!(*blocked_dofs_val)) { *u_val += (1. - alpha) * delta_t * *u_dot_val; } u_val++; u_dot_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GeneralizedTrapezoidal::corrector(const SolutionType & type, Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); switch (type) { case _temperature: this->allCorrector<_temperature>(delta_t, u, u_dot, blocked_dofs, delta); break; case _temperature_rate: this->allCorrector<_temperature_rate>(delta_t, u, u_dot, blocked_dofs, delta); break; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real GeneralizedTrapezoidal::getTemperatureCoefficient( const SolutionType & type, Real delta_t) const { switch (type) { case _temperature: return 1.; case _temperature_rate: return alpha * delta_t; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ Real GeneralizedTrapezoidal::getTemperatureRateCoefficient( const SolutionType & type, Real delta_t) const { switch (type) { case _temperature: return 1. / (alpha * delta_t); case _temperature_rate: return 1.; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ template void GeneralizedTrapezoidal::allCorrector(Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.getSize(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real e = getTemperatureCoefficient(type, delta_t); Real d = getTemperatureRateCoefficient(type, delta_t); Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); Real * delta_val = delta.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { if (!(*blocked_dofs_val)) { *u_val += e ** delta_val; *u_dot_val += d ** delta_val; } u_val++; u_dot_val++; delta_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GeneralizedTrapezoidal::assembleJacobian(const SolutionType & type, Real delta_t) { AKANTU_DEBUG_IN(); SparseMatrix & J = this->dof_manager.getMatrix("J"); const SparseMatrix & M = this->dof_manager.getMatrix("M"); const SparseMatrix & K = this->dof_manager.getMatrix("K"); + bool does_j_need_update = false; + does_j_need_update |= M.getRelease() != m_release; + does_j_need_update |= K.getRelease() != k_release; + if (!does_j_need_update) { + AKANTU_DEBUG_OUT(); + return; + } + + J.clear(); + Real c = this->getTemperatureRateCoefficient(type, delta_t); Real e = this->getTemperatureCoefficient(type, delta_t); J.add(M, e); J.add(K, c); + m_release = M.getRelease(); + k_release = K.getRelease(); + AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/model/integration_scheme/generalized_trapezoidal.hh b/src/model/integration_scheme/generalized_trapezoidal.hh index 873d2c7d6..7a2256941 100644 --- a/src/model/integration_scheme/generalized_trapezoidal.hh +++ b/src/model/integration_scheme/generalized_trapezoidal.hh @@ -1,160 +1,166 @@ /** * @file generalized_trapezoidal.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jul 04 2011 * @date last modification: Fri Oct 23 2015 * * @brief Generalized Trapezoidal Method. This implementation is taken from * Méthodes numériques en mécanique des solides by Alain Curnier \note{ISBN: * 2-88074-247-1} * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__ #define __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__ #include "integration_scheme_1st_order.hh" __BEGIN_AKANTU__ /** * The two differentiate equation (thermal and kinematic) are : * @f{eqnarray*}{ * C\dot{u}_{n+1} + Ku_{n+1} = q_{n+1}\\ * u_{n+1} = u_{n} + (1-\alpha) \Delta t \dot{u}_{n} + \alpha \Delta t *\dot{u}_{n+1} * @f} * * To solve it : * Predictor : * @f{eqnarray*}{ * u^0_{n+1} &=& u_{n} + (1-\alpha) \Delta t v_{n} \\ * \dot{u}^0_{n+1} &=& \dot{u}_{n} * @f} * * Solve : * @f[ (a C + b K^i_{n+1}) w = q_{n+1} - f^i_{n+1} - C \dot{u}^i_{n+1} @f] * * Corrector : * @f{eqnarray*}{ * \dot{u}^{i+1}_{n+1} &=& \dot{u}^{i}_{n+1} + a w \\ * u^{i+1}_{n+1} &=& u^{i}_{n+1} + b w * @f} * * a and b depends on the resolution method : temperature (u) or temperature *rate (@f$\dot{u}@f$) * * For temperature : @f$ w = \delta u, a = 1 / (\alpha \Delta t) , b = 1 @f$ @n * For temperature rate : @f$ w = \delta \dot{u}, a = 1, b = \alpha \Delta t @f$ */ class GeneralizedTrapezoidal : public IntegrationScheme1stOrder { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GeneralizedTrapezoidal(DOFManager & dof_manager, const ID & dof_id, Real alpha = 0); virtual ~GeneralizedTrapezoidal(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void predictor(Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs) const; virtual void corrector(const SolutionType & type, Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs, const Array & delta) const; virtual void assembleJacobian(const SolutionType & type, Real time_step); public: /// the coeffichent @f{b@f} in the description Real getTemperatureCoefficient(const SolutionType & type, Real delta_t) const; /// the coeffichent @f{a@f} in the description Real getTemperatureRateCoefficient(const SolutionType & type, Real delta_t) const; private: template void allCorrector(Real delta_t, Array & u, Array & u_dot, const Array & blocked_dofs, const Array & delta) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Alpha, alpha, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// the @f$\alpha@f$ parameter Real alpha; + + /// last release of M matrix + UInt m_release; + + /// last release of K matrix + UInt k_release; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Forward Euler (explicit) -> condition on delta_t */ class ForwardEuler : public GeneralizedTrapezoidal { public: ForwardEuler(DOFManager & dof_manager, const ID & dof_id) : GeneralizedTrapezoidal(dof_manager, dof_id, 0.){}; }; /** * Trapezoidal rule (implicit), midpoint rule or Crank-Nicolson */ class TrapezoidalRule1 : public GeneralizedTrapezoidal { public: TrapezoidalRule1(DOFManager & dof_manager, const ID & dof_id) : GeneralizedTrapezoidal(dof_manager, dof_id, .5){}; }; /** * Backward Euler (implicit) */ class BackwardEuler : public GeneralizedTrapezoidal { public: BackwardEuler(DOFManager & dof_manager, const ID & dof_id) : GeneralizedTrapezoidal(dof_manager, dof_id, 1.){}; }; /* -------------------------------------------------------------------------- */ __END_AKANTU__ #endif /* __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__ */ diff --git a/src/model/integration_scheme/newmark-beta.cc b/src/model/integration_scheme/newmark-beta.cc index a2669b46d..9527dfc6c 100644 --- a/src/model/integration_scheme/newmark-beta.cc +++ b/src/model/integration_scheme/newmark-beta.cc @@ -1,240 +1,260 @@ /** * @file newmark-beta.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Tue Oct 05 2010 * @date last modification: Thu Jun 05 2014 * * @brief implementation of the newmark-@f$\beta@f$ integration scheme. This * implementation is taken from Méthodes numériques en mécanique des solides by * Alain Curnier \note{ISBN: 2-88074-247-1} * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "newmark-beta.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ NewmarkBeta::NewmarkBeta(DOFManager & dof_manager, const ID & dof_id, Real alpha, Real beta) : IntegrationScheme2ndOrder(dof_manager, dof_id), beta(beta), alpha(alpha), - k(0.), h(0.){ + k(0.), h(0.), m_release(0), k_release(0), c_release(0) { - this->registerParam("alpha", this->alpha, alpha, _pat_parsmod, "The alpha parameter"); - this->registerParam("beta", this->beta, beta, _pat_parsmod, "The beta parameter"); + this->registerParam("alpha", this->alpha, alpha, _pat_parsmod, + "The alpha parameter"); + this->registerParam("beta", this->beta, beta, _pat_parsmod, + "The beta parameter"); } /* -------------------------------------------------------------------------- */ /* * @f$ \tilde{u_{n+1}} = u_{n} + \Delta t \dot{u}_n + \frac{\Delta t^2}{2} * \ddot{u}_n @f$ * @f$ \tilde{\dot{u}_{n+1}} = \dot{u}_{n} + \Delta t \ddot{u}_{n} @f$ * @f$ \tilde{\ddot{u}_{n}} = \ddot{u}_{n} @f$ */ void NewmarkBeta::predictor(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.getSize(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); Real * u_dot_dot_val = u_dot_dot.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt d = 0; d < nb_degree_of_freedom; d++) { if (!(*blocked_dofs_val)) { Real dt_a_n = delta_t * *u_dot_dot_val; *u_val += (1 - k * alpha) * delta_t * *u_dot_val + (.5 - h * alpha * beta) * delta_t * dt_a_n; *u_dot_val = (1 - k) * *u_dot_val + (1 - h * beta) * dt_a_n; *u_dot_dot_val = (1 - h) * *u_dot_dot_val; } u_val++; u_dot_val++; u_dot_dot_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NewmarkBeta::corrector(const SolutionType & type, Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); switch (type) { case _acceleration: { this->allCorrector<_acceleration>(delta_t, u, u_dot, u_dot_dot, blocked_dofs, delta); break; } case _velocity: { this->allCorrector<_velocity>(delta_t, u, u_dot, u_dot_dot, blocked_dofs, delta); break; } case _displacement: { this->allCorrector<_displacement>(delta_t, u, u_dot, u_dot_dot, blocked_dofs, delta); break; } default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real NewmarkBeta::getAccelerationCoefficient(const SolutionType & type, Real delta_t) const { switch (type) { case _acceleration: return 1.; case _velocity: return 1. / (beta * delta_t); case _displacement: return 1. / (alpha * beta * delta_t * delta_t); default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ Real NewmarkBeta::getVelocityCoefficient(const SolutionType & type, Real delta_t) const { switch (type) { case _acceleration: return beta * delta_t; case _velocity: return 1.; case _displacement: return 1. / (alpha * delta_t); default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ Real NewmarkBeta::getDisplacementCoefficient(const SolutionType & type, Real delta_t) const { switch (type) { case _acceleration: return alpha * beta * delta_t * delta_t; case _velocity: return alpha * delta_t; case _displacement: return 1.; default: AKANTU_EXCEPTION("The corrector type : " << type << " is not supported by this type of integration scheme"); } } /* -------------------------------------------------------------------------- */ template void NewmarkBeta::allCorrector(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const { AKANTU_DEBUG_IN(); UInt nb_nodes = u.getSize(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Real c = getAccelerationCoefficient(type, delta_t); Real d = getVelocityCoefficient(type, delta_t); Real e = getDisplacementCoefficient(type, delta_t); Real * u_val = u.storage(); Real * u_dot_val = u_dot.storage(); Real * u_dot_dot_val = u_dot_dot.storage(); Real * delta_val = delta.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); for (UInt dof = 0; dof < nb_degree_of_freedom; dof++) { if (!(*blocked_dofs_val)) { - *u_val += e ** delta_val; - *u_dot_val += d ** delta_val; - *u_dot_dot_val += c ** delta_val; + *u_val += e * *delta_val; + *u_dot_val += d * *delta_val; + *u_dot_dot_val += c * *delta_val; } u_val++; u_dot_val++; u_dot_dot_val++; delta_val++; blocked_dofs_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NewmarkBeta::assembleJacobian(const SolutionType & type, Real delta_t) { AKANTU_DEBUG_IN(); SparseMatrix & J = this->dof_manager.getMatrix("J"); const SparseMatrix & M = this->dof_manager.getMatrix("M"); const SparseMatrix & K = this->dof_manager.getMatrix("K"); + bool does_j_need_update = false; + does_j_need_update |= M.getRelease() != m_release; + does_j_need_update |= K.getRelease() != k_release; + if (this->dof_manager.hasMatrix("C")) { + const SparseMatrix & C = this->dof_manager.getMatrix("C"); + does_j_need_update |= C.getRelease() != c_release; + } + + if (!does_j_need_update) { + AKANTU_DEBUG_OUT(); + return; + } + + J.clear(); + Real c = this->getAccelerationCoefficient(type, delta_t); Real e = this->getDisplacementCoefficient(type, delta_t); if (!(e == 0.)) { // in explicit this coefficient is exactly 0. J.add(K, e); } J.add(M, c); - try { + m_release = M.getRelease(); + k_release = K.getRelease(); + + if (this->dof_manager.hasMatrix("C")) { Real d = this->getVelocityCoefficient(type, delta_t); const SparseMatrix & C = this->dof_manager.getMatrix("C"); J.add(C, d); - } catch (debug::Exception &) { + c_release = C.getRelease(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/integration_scheme/newmark-beta.hh b/src/model/integration_scheme/newmark-beta.hh index e2b191724..6272757e4 100644 --- a/src/model/integration_scheme/newmark-beta.hh +++ b/src/model/integration_scheme/newmark-beta.hh @@ -1,186 +1,195 @@ /** * @file newmark-beta.hh * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Tue Oct 05 2010 * @date last modification: Fri Oct 23 2015 * * @brief implementation of the newmark-@f$\beta@f$ integration scheme. This * implementation is taken from Méthodes numériques en mécanique des solides by * Alain Curnier \note{ISBN: 2-88074-247-1} * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NEWMARK_BETA_HH__ #define __AKANTU_NEWMARK_BETA_HH__ /* -------------------------------------------------------------------------- */ #include "integration_scheme_2nd_order.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * The three differentiate equations (dynamic and cinematic) are : * @f{eqnarray*}{ * M \ddot{u}_{n+1} + C \dot{u}_{n+1} + K u_{n+1} &=& q_{n+1} \\ * u_{n+1} &=& u_{n} + (1 - \alpha) \Delta t \dot{u}_{n} + \alpha \Delta t *\dot{u}_{n+1} + (1/2 - \alpha) \Delta t^2 \ddot{u}_n \\ * \dot{u}_{n+1} &=& \dot{u}_{n} + (1 - \beta) \Delta t \ddot{u}_{n} + \beta *\Delta t \ddot{u}_{n+1} * @f} * * Predictor: * @f{eqnarray*}{ * u^{0}_{n+1} &=& u_{n} + \Delta t \dot{u}_n + \frac{\Delta t^2}{2} *\ddot{u}_n \\ * \dot{u}^{0}_{n+1} &=& \dot{u}_{n} + \Delta t \ddot{u}_{n} \\ * \ddot{u}^{0}_{n+1} &=& \ddot{u}_{n} * @f} * * Solve : * @f[ (c M + d C + e K^i_{n+1}) w = = q_{n+1} - f^i_{n+1} - C \dot{u}^i_{n+1} *- M \ddot{u}^i_{n+1} @f] * * Corrector : * @f{eqnarray*}{ * \ddot{u}^{i+1}_{n+1} &=& \ddot{u}^{i}_{n+1} + c w \\ * \dot{u}^{i+1}_{n+1} &=& \dot{u}^{i}_{n+1} + d w \\ * u^{i+1}_{n+1} &=& u^{i}_{n+1} + e w * @f} * * c, d and e are parameters depending on the method used to solve the equations *@n * For acceleration : @f$ w = \delta \ddot{u}, e = \alpha \beta \Delta t^2, d = *\beta \Delta t, c = 1 @f$ @n * For velocity : @f$ w = \delta \dot{u}, e = 1/\beta \Delta t, d = *1, c = \alpha \Delta t @f$ @n * For displacement : @f$ w = \delta u, e = 1, d = *1/\alpha \Delta t, c = 1/\alpha \beta \Delta t^2 @f$ */ class NewmarkBeta : public IntegrationScheme2ndOrder { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NewmarkBeta(DOFManager & dof_manager, const ID & dof_id, Real alpha = 0., Real beta = 0.); ~NewmarkBeta(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void predictor(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs) const; void corrector(const SolutionType & type, Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const; void assembleJacobian(const SolutionType & type, Real delta_t); public: Real getAccelerationCoefficient(const SolutionType & type, Real delta_t) const; Real getVelocityCoefficient(const SolutionType & type, Real delta_t) const; Real getDisplacementCoefficient(const SolutionType & type, Real delta_t) const; private: template void allCorrector(Real delta_t, Array & u, Array & u_dot, Array & u_dot_dot, const Array & blocked_dofs, const Array & delta) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Beta, beta, Real); AKANTU_GET_MACRO(Alpha, alpha, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the @f$\beta@f$ parameter Real beta; /// the @f$\alpha@f$ parameter Real alpha; Real k; Real h; + + /// last release of M matrix + UInt m_release; + + /// last release of K matrix + UInt k_release; + + /// last release of C matrix + UInt c_release; }; /** * central difference method (explicit) * undamped stability condition : * @f$ \Delta t = \alpha \Delta t_{crit} = \frac{2}{\omega_{max}} \leq \min_{e} *\frac{l_e}{c_e} * */ class CentralDifference : public NewmarkBeta { public: CentralDifference(DOFManager & dof_manager, const ID & dof_id) : NewmarkBeta(dof_manager, dof_id, 0., 1. / 2.){}; }; //#include "integration_scheme/central_difference.hh" /// undamped trapezoidal rule (implicit) class TrapezoidalRule2 : public NewmarkBeta { public: TrapezoidalRule2(DOFManager & dof_manager, const ID & dof_id) : NewmarkBeta(dof_manager, dof_id, 1. / 2., 1. / 2.){}; }; /// Fox-Goodwin rule (implicit) class FoxGoodwin : public NewmarkBeta { public: FoxGoodwin(DOFManager & dof_manager, const ID & dof_id) : NewmarkBeta(dof_manager, dof_id, 1. / 6., 1. / 2.){}; }; /// Linear acceleration (implicit) class LinearAceleration : public NewmarkBeta { public: LinearAceleration(DOFManager & dof_manager, const ID & dof_id) : NewmarkBeta(dof_manager, dof_id, 1. / 3., 1. / 2.){}; }; /* -------------------------------------------------------------------------- */ __END_AKANTU__ #endif /* __AKANTU_NEWMARK_BETA_HH__ */ diff --git a/src/model/integration_scheme/pseudo_time.cc b/src/model/integration_scheme/pseudo_time.cc index f0c38ad92..c97fef22a 100644 --- a/src/model/integration_scheme/pseudo_time.cc +++ b/src/model/integration_scheme/pseudo_time.cc @@ -1,84 +1,88 @@ /** * @file pseudo_time.cc * * @author Nicolas Richart * * @date Wed Feb 17 09:49:10 2016 * * @brief Implementation of a really simple integration scheme * * @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 "pseudo_time.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ PseudoTime::PseudoTime(DOFManager & dof_manager, const ID & dof_id) - : IntegrationScheme(dof_manager, dof_id, 0) {} + : IntegrationScheme(dof_manager, dof_id, 0), k_release(0) {} /* -------------------------------------------------------------------------- */ -void PseudoTime::predictor(__attribute__((unused)) Real delta_t) {} +void PseudoTime::predictor(Real) {} /* -------------------------------------------------------------------------- */ -void PseudoTime::corrector(__attribute__((unused)) const SolutionType & type, - __attribute__((unused)) Real delta_t) { +void PseudoTime::corrector(const SolutionType &, Real) { Array & u = this->dof_manager.getDOFs(this->dof_id); const Array & delta = this->dof_manager.getSolution(this->dof_id); const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); UInt nb_nodes = u.getSize(); UInt nb_degree_of_freedom = u.getNbComponent() * nb_nodes; Array::scalar_iterator u_it = u.begin_reinterpret(nb_degree_of_freedom); Array::scalar_iterator u_end = u.end_reinterpret(nb_degree_of_freedom); Array::const_scalar_iterator delta_it = delta.begin_reinterpret(nb_degree_of_freedom); Array::const_scalar_iterator blocked_dofs_it = blocked_dofs.begin_reinterpret(nb_degree_of_freedom); for (; u_it != u_end; ++u_it, ++delta_it, ++blocked_dofs_it) { if (!(*blocked_dofs_it)) { *u_it += *delta_it; } } } /* -------------------------------------------------------------------------- */ -void PseudoTime::assembleJacobian(__attribute__((unused)) const SolutionType & type, - __attribute__((unused)) Real delta_t) { +void PseudoTime::assembleJacobian(const SolutionType &, Real) { SparseMatrix & J = this->dof_manager.getMatrix("J"); const SparseMatrix & K = this->dof_manager.getMatrix("K"); + if (K.getRelease() == k_release) + return; + + J.clear(); J.add(K); + + k_release = K.getRelease(); } /* -------------------------------------------------------------------------- */ -void PseudoTime::assembleResidual(__attribute__((unused)) bool is_lumped) {} +void PseudoTime::assembleResidual(bool) {} /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/integration_scheme/pseudo_time.hh b/src/model/integration_scheme/pseudo_time.hh index c409ab342..9be2d22d3 100644 --- a/src/model/integration_scheme/pseudo_time.hh +++ b/src/model/integration_scheme/pseudo_time.hh @@ -1,68 +1,70 @@ /** * @file pseudo_time.hh * * @author Nicolas Richart * * @date Wed Feb 17 09:46:05 2016 * * @brief Pseudo time integration scheme * * @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 "integration_scheme.hh" /* -------------------------------------------------------------------------- */ - #ifndef __AKANTU_PSEUDO_TIME_HH__ #define __AKANTU_PSEUDO_TIME_HH__ __BEGIN_AKANTU__ class PseudoTime : public IntegrationScheme { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PseudoTime(DOFManager & dof_manager, const ID & dof_id); virtual ~PseudoTime() {} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// generic interface of a predictor virtual void predictor(Real delta_t); /// generic interface of a corrector virtual void corrector(const SolutionType & type, Real delta_t); /// assemble the jacobian matrix - virtual void assembleJacobian(const SolutionType & type, - Real delta_t); + virtual void assembleJacobian(const SolutionType & type, Real delta_t); /// assemble the residual virtual void assembleResidual(bool is_lumped); + +protected: + /// last release of K matrix + UInt k_release; }; __END_AKANTU__ #endif /* __AKANTU_PSEUDO_TIME_HH__ */ diff --git a/src/model/non_linear_solver_linear.hh b/src/model/non_linear_solver_linear.hh index e76faf320..7d14507b9 100644 --- a/src/model/non_linear_solver_linear.hh +++ b/src/model/non_linear_solver_linear.hh @@ -1,76 +1,76 @@ /** * @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 "solver_mumps.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); /* ------------------------------------------------------------------------ */ /* 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_lumped.hh b/src/model/non_linear_solver_lumped.hh index f358fbb09..03ec67c44 100644 --- a/src/model/non_linear_solver_lumped.hh +++ b/src/model/non_linear_solver_lumped.hh @@ -1,82 +1,82 @@ /** * @file non_linear_solver_lumped.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 "solver_mumps.hh" +#include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NON_LINEAR_SOLVER_LUMPED_HH__ #define __AKANTU_NON_LINEAR_SOLVER_LUMPED_HH__ namespace akantu { class DOFManagerDefault; } __BEGIN_AKANTU__ class NonLinearSolverLumped : public NonLinearSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NonLinearSolverLumped(DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id = "non_linear_solver_lumped", UInt memory_id = 0); virtual ~NonLinearSolverLumped(); /* ------------------------------------------------------------------------ */ /* 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); static void solveLumped(const Array & A, Array & x, const Array & b, const Array & blocked_dofs, Real alpha); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: DOFManagerDefault & dof_manager; /// Coefficient to apply between x and A^{-1} b Real alpha; }; __END_AKANTU__ #endif /* __AKANTU_NON_LINEAR_SOLVER_LUMPED_HH__ */ diff --git a/src/model/non_linear_solver_newton_raphson.cc b/src/model/non_linear_solver_newton_raphson.cc index 213874d73..89ca4568b 100644 --- a/src/model/non_linear_solver_newton_raphson.cc +++ b/src/model/non_linear_solver_newton_raphson.cc @@ -1,184 +1,184 @@ /** * @file non_linear_solver_newton_raphson.cc * * @author Nicolas Richart * * @date Tue Aug 25 00:57:00 2015 * * @brief Implementation of the default NonLinearSolver * * @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_newton_raphson.hh" #include "dof_manager_default.hh" #include "solver_callback.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::NonLinearSolverNewtonRaphson( DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id, UInt memory_id) : NonLinearSolver(dof_manager, non_linear_solver_type, id, memory_id), dof_manager(dof_manager), solver(dof_manager, "J", id + ":sparse_solver", memory_id), n_iter(0), error(0.), converged(false) { this->supported_type.insert(_nls_newton_raphson_modified); this->supported_type.insert(_nls_newton_raphson); this->supported_type.insert(_nls_linear); this->checkIfTypeIsSupported(); this->registerParam("threshold", convergence_criteria, 1e-10, _pat_parsmod, "Threshold to consider results as converged"); this->registerParam("convergence_type", convergence_criteria_type, _scc_solution, _pat_parsmod, "Type of convergence criteria"); this->registerParam("max_iterations", max_iterations, 10, _pat_parsmod, "Max number of iterations"); this->registerParam("error", error, _pat_readable, "Last reached error"); this->registerParam("nb_iterations", n_iter, _pat_readable, "Last reached number of iterations"); this->registerParam("converged", converged, _pat_readable, "Did last solve converged"); this->registerParam("force_linear_recompute", force_linear_recompute, true, _pat_modifiable, "Force reassembly of the jacobian matrix"); } /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::~NonLinearSolverNewtonRaphson() {} /* ------------------------------------------------------------------------ */ void NonLinearSolverNewtonRaphson::solve(SolverCallback & solver_callback) { solver_callback.predictor(); solver_callback.assembleResidual(); if (this->non_linear_solver_type == _nls_newton_raphson_modified || (this->non_linear_solver_type == _nls_linear && this->force_linear_recompute)) { solver_callback.assembleJacobian(); this->force_linear_recompute = false; } this->n_iter = 0; this->converged = false; if (this->convergence_criteria_type == _scc_residual) { this->converged = this->testConvergence(this->dof_manager.getResidual()); if (this->converged) return; } do { if (this->non_linear_solver_type == _nls_newton_raphson) solver_callback.assembleJacobian(); this->solver.solve(); solver_callback.corrector(); // EventManager::sendEvent(NonLinearSolver::AfterSparseSolve(method)); if (this->convergence_criteria_type == _scc_residual) { solver_callback.assembleResidual(); this->converged = this->testConvergence(this->dof_manager.getResidual()); } else { this->converged = this->testConvergence(this->dof_manager.getGlobalSolution()); } if (this->convergence_criteria_type == _scc_solution && !this->converged) solver_callback.assembleResidual(); // this->dump(); this->n_iter++; AKANTU_DEBUG_INFO( "[" << this->convergence_criteria_type << "] Convergence iteration " << std::setw(std::log10(this->max_iterations)) << this->n_iter << ": error " << this->error << (this->converged ? " < " : " > ") << this->convergence_criteria); } while (!this->converged && this->n_iter < this->max_iterations); // this makes sure that you have correct strains and stresses after the // solveStep function (e.g., for dumping) if (this->convergence_criteria_type == _scc_solution) solver_callback.assembleResidual(); if (this->converged) { // EventManager::sendEvent( // SolidMechanicsModelEvent::AfterNonLinearSolverSolves(method)); } else if (this->n_iter == this->max_iterations) { AKANTU_CUSTOM_EXCEPTION(debug::NLSNotConvergedException( this->convergence_criteria, this->n_iter)); AKANTU_DEBUG_WARNING("[" << this->convergence_criteria_type << "] Convergence not reached after " << std::setw(std::log10(this->max_iterations)) << this->n_iter << " iteration" << (this->n_iter == 1 ? "" : "s") << "!"); } return; } /* -------------------------------------------------------------------------- */ bool NonLinearSolverNewtonRaphson::testConvergence(const Array & array) { AKANTU_DEBUG_IN(); const Array & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); UInt nb_degree_of_freedoms = array.getSize(); - Array::const_scalar_iterator arr_it = array.begin(); - Array::const_scalar_iterator bld_it = blocked_dofs.begin(); + auto arr_it = array.begin(); + auto bld_it = blocked_dofs.begin(); Real norm = 0.; for (UInt n = 0; n < nb_degree_of_freedoms; ++n, ++arr_it, ++bld_it) { bool is_local_node = this->dof_manager.isLocalOrMasterDOF(n); - if ((*bld_it == 0) && is_local_node) { + if ((! *bld_it) && is_local_node) { norm += *arr_it * *arr_it; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, _so_sum); norm = std::sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something went wrong in the solve phase"); this->error = norm; return (error < this->convergence_criteria); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/non_linear_solver_newton_raphson.hh b/src/model/non_linear_solver_newton_raphson.hh index f379802c4..d3e929d49 100644 --- a/src/model/non_linear_solver_newton_raphson.hh +++ b/src/model/non_linear_solver_newton_raphson.hh @@ -1,101 +1,101 @@ /** * @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 "solver_mumps.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); 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/model/time_step_solvers/time_step_solver.cc b/src/model/time_step_solvers/time_step_solver.cc index 67d0a9cfb..4df58122a 100644 --- a/src/model/time_step_solvers/time_step_solver.cc +++ b/src/model/time_step_solvers/time_step_solver.cc @@ -1,93 +1,93 @@ /** * @file time_step_solver.cc * * @author Nicolas Richart * * @date Mon Oct 12 16:56:43 2015 * * @brief Implementation of common part of TimeStepSolvers * * @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 "time_step_solver.hh" #include "non_linear_solver.hh" #include "dof_manager.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ TimeStepSolver::TimeStepSolver(DOFManager & dof_manager, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, const ID & id, UInt memory_id) : Memory(id, memory_id), SolverCallback(dof_manager), _dof_manager(dof_manager), type(type), time_step(0.), solver_callback(NULL), non_linear_solver(non_linear_solver) { this->registerSubRegistry("non_linear_solver", non_linear_solver); } /* -------------------------------------------------------------------------- */ TimeStepSolver::~TimeStepSolver() {} /* -------------------------------------------------------------------------- */ void TimeStepSolver::predictor() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); this->solver_callback->predictor(); } /* -------------------------------------------------------------------------- */ void TimeStepSolver::corrector() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); this->solver_callback->corrector(); } /* -------------------------------------------------------------------------- */ void TimeStepSolver::assembleJacobian() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); - this->_dof_manager.clearMatrix("J"); + // this->_dof_manager.clearMatrix("J"); this->solver_callback->assembleJacobian(); } /* -------------------------------------------------------------------------- */ void TimeStepSolver::assembleResidual() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); this->_dof_manager.clearResidual(); this->solver_callback->assembleResidual(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/solver/sparse_matrix.cc b/src/solver/sparse_matrix.cc index 4b9069d7c..ee229bde1 100644 --- a/src/solver/sparse_matrix.cc +++ b/src/solver/sparse_matrix.cc @@ -1,82 +1,80 @@ /** * @file sparse_matrix.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Mon Nov 16 2015 * * @brief implementation of the SparseMatrix class * * @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 /* -------------------------------------------------------------------------- */ +#include "dof_manager.hh" #include "sparse_matrix.hh" #include "static_communicator.hh" -#include "dof_manager.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(DOFManager & dof_manager, const MatrixType & matrix_type, const ID & id) : id(id), _dof_manager(dof_manager), matrix_type(matrix_type), - size(dof_manager.getSystemSize()), nb_non_zero(0), matrix_release(1) { + size(dof_manager.getSystemSize()), nb_non_zero(0) { AKANTU_DEBUG_IN(); - StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + const auto & comm = _dof_manager.getCommunicator(); this->nb_proc = comm.getNbProc(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(const SparseMatrix & matrix, const ID & id) - : id(id), _dof_manager(matrix._dof_manager), - matrix_type(matrix.matrix_type), size(matrix.size), - nb_proc(matrix.nb_proc), nb_non_zero(matrix.nb_non_zero), - matrix_release(1) {} + : SparseMatrix(matrix._dof_manager, matrix.matrix_type, id) { + nb_non_zero = matrix.nb_non_zero; +} /* -------------------------------------------------------------------------- */ SparseMatrix::~SparseMatrix() {} /* -------------------------------------------------------------------------- */ Array & operator*=(Array & vect, const SparseMatrix & mat) { Array tmp(vect.getSize(), vect.getNbComponent(), 0.); mat.matVecMul(vect, tmp); vect.copy(tmp); return vect; } /* -------------------------------------------------------------------------- */ void SparseMatrix::add(const SparseMatrix & B, Real alpha) { B.addMeTo(*this, alpha); } /* -------------------------------------------------------------------------- */ - __END_AKANTU__ diff --git a/src/solver/sparse_matrix.hh b/src/solver/sparse_matrix.hh index 5391595ef..760010a33 100644 --- a/src/solver/sparse_matrix.hh +++ b/src/solver/sparse_matrix.hh @@ -1,162 +1,161 @@ /** * @file sparse_matrix.hh * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Fri Oct 16 2015 * * @brief sparse matrix storage class (distributed assembled matrix) * This is a COO format (Coordinate List) * * @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" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SPARSE_MATRIX_HH__ #define __AKANTU_SPARSE_MATRIX_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class DOFManager; class TermsToAssemble; } namespace akantu { class SparseMatrix { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseMatrix(DOFManager & dof_manager, const MatrixType & matrix_type, const ID & id = "sparse_matrix"); SparseMatrix(const SparseMatrix & matrix, const ID & id = "sparse_matrix"); virtual ~SparseMatrix(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// remove the existing profile virtual void clearProfile(); /// set the matrix to 0 virtual void clear() = 0; /// add a non-zero element to the profile - virtual inline UInt addToProfile(UInt i, UInt j) = 0; + virtual UInt addToProfile(UInt i, UInt j) = 0; /// assemble a local matrix in the sparse one virtual void addToMatrix(UInt i, UInt j, Real value) = 0; /// save the profil in a file using the MatrixMarket file format virtual void saveProfile(__attribute__((unused)) const std::string &) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// save the matrix in a file using the MatrixMarket file format virtual void saveMatrix(__attribute__((unused)) const std::string &) const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// multiply the matrix by a coefficient virtual void mul(Real alpha) = 0; /// add matrices virtual void add(const SparseMatrix & matrix, Real alpha = 1.); /// Equivalent of *gemv in blas virtual void matVecMul(const Array & x, Array & y, Real alpha = 1., Real beta = 0.) const = 0; /// modify the matrix to "remove" the blocked dof virtual void applyBoundary(Real block_val = 1.) = 0; /// operator *= SparseMatrix & operator*=(Real alpha) { this->mul(alpha); return *this; } protected: /// This is the revert of add B += \alpha * *this; virtual void addMeTo(SparseMatrix & B, Real alpha) const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the values at potition i, j virtual inline Real operator()(__attribute__((unused)) UInt i, __attribute__((unused)) UInt j) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// return the values at potition i, j virtual inline Real & operator()(__attribute__((unused)) UInt i, __attribute__((unused)) UInt j) { AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_GET_MACRO(NbNonZero, nb_non_zero, UInt); AKANTU_GET_MACRO(Size, size, UInt); AKANTU_GET_MACRO(MatrixType, matrix_type, const MatrixType &); + virtual UInt getRelease() const = 0; + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: ID id; /// Underlying dof manager DOFManager & _dof_manager; /// sparce matrix type MatrixType matrix_type; /// Size of the matrix UInt size; /// number of processors UInt nb_proc; /// number of non zero element UInt nb_non_zero; - - /// matrix definition release - UInt matrix_release; }; Array & operator*=(Array & vect, const SparseMatrix & mat); } // akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix_inline_impl.cc" #endif /* __AKANTU_SPARSE_MATRIX_HH__ */ diff --git a/src/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc index e5c7fb56d..760148ccb 100644 --- a/src/solver/sparse_matrix_aij.cc +++ b/src/solver/sparse_matrix_aij.cc @@ -1,228 +1,228 @@ /** * @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 /* -------------------------------------------------------------------------- */ __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(0), value_release(0) {} + 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(0), value_release(0) {} + 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(); Array::const_scalar_iterator irn_val = irn.begin(); Array::const_scalar_iterator jcn_val = jcn.begin(); Array::scalar_iterator a_val = a.begin(); 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); 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; } } } ++irn_val; ++jcn_val; ++a_val; } this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); 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(); std::ofstream outfile; outfile.precision(std::numeric_limits::digits10); outfile.open(filename.c_str()); 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; for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn(i) << " " << this->jcn(i) << " " << this->a(i) << std::endl; } 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(); for (UInt k = 0; k < this->nb_non_zero; ++k, ++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 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) { 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)) { 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/sparse_matrix_aij.hh b/src/solver/sparse_matrix_aij.hh index dcf3f1b0e..6a4e43cf4 100644 --- a/src/solver/sparse_matrix_aij.hh +++ b/src/solver/sparse_matrix_aij.hh @@ -1,179 +1,179 @@ /** * @file sparse_matrix_aij.hh * * @author Nicolas Richart * * @date Mon Aug 17 21:22:57 2015 * * @brief AIJ implementation of the SparseMatrix (this the format used by Mumps) * * @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_array.hh" #include "aka_common.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SPARSE_MATRIX_AIJ_HH__ #define __AKANTU_SPARSE_MATRIX_AIJ_HH__ namespace akantu { class DOFManagerDefault; class TermsToAssemble; } __BEGIN_AKANTU__ class SparseMatrixAIJ : public SparseMatrix { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseMatrixAIJ(DOFManagerDefault & dof_manager, const MatrixType & matrix_type, const ID & id = "sparse_matrix_aij"); SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id = "sparse_matrix_aij"); virtual ~SparseMatrixAIJ(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// remove the existing profile inline void clearProfile(); /// add a non-zero element - virtual UInt addToProfile(UInt i, UInt j); + virtual inline UInt addToProfile(UInt i, UInt j); /// set the matrix to 0 virtual void clear(); /// assemble a local matrix in the sparse one - inline void addToMatrix(UInt i, UInt j, Real value); + virtual inline void addToMatrix(UInt i, UInt j, Real value); /// set the size of the matrix void resize(UInt size) { this->size = size; } /// modify the matrix to "remove" the blocked dof virtual void applyBoundary(Real block_val = 1.); /// save the profil in a file using the MatrixMarket file format virtual void saveProfile(const std::string & filename) const; /// save the matrix in a file using the MatrixMarket file format virtual void saveMatrix(const std::string & filename) const; /// copy assuming the profile are the same virtual void copyContent(const SparseMatrix & matrix); /// multiply the matrix by a scalar void mul(Real alpha); /// add matrix *this += B // virtual void add(const SparseMatrix & matrix, Real alpha); /// Equivalent of *gemv in blas virtual void matVecMul(const Array & x, Array & y, Real alpha = 1., Real beta = 0.) const; /* ------------------------------------------------------------------------ */ /// accessor to A_{ij} - if (i, j) not present it returns 0 inline Real operator()(UInt i, UInt j) const; /// accessor to A_{ij} - if (i, j) not present it fails, (i, j) should be /// first added to the profile inline Real & operator()(UInt i, UInt j); protected: /// This is the revert of add B += \alpha * *this; virtual void addMeTo(SparseMatrix & B, Real alpha) const; private: /// This is just to inline the addToMatrix function template void addMeToTemplated(MatrixType & B, Real alpha) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(IRN, irn, const Array &); AKANTU_GET_MACRO(JCN, jcn, const Array &); AKANTU_GET_MACRO(A, a, const Array &); /// The release changes at each call of a function that changes the profile, /// it in increasing but could overflow so it should be checked as /// (my_release != release) and not as (my_release < release) AKANTU_GET_MACRO(ProfileRelease, profile_release, UInt); AKANTU_GET_MACRO(ValueRelease, value_release, UInt); - + virtual UInt getRelease() const { return value_release; } protected: typedef std::pair KeyCOO; typedef unordered_map::type coordinate_list_map; /// get the pair corresponding to (i, j) inline KeyCOO key(UInt i, UInt j) const { if (this->matrix_type == _symmetric && (i > j)) return std::make_pair(j, i); return std::make_pair(i, j); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: DOFManagerDefault & dof_manager; /// row indexes Array irn; /// column indexes Array jcn; /// values : A[k] = Matrix[irn[k]][jcn[k]] Array a; /// Profile release UInt profile_release; /// Value release UInt value_release; /* * map for (i, j) -> k correspondence */ coordinate_list_map irn_jcn_k; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix_aij_inline_impl.cc" #endif /* __AKANTU_SPARSE_MATRIX_AIJ_HH__ */ diff --git a/src/solver/sparse_solver_mumps.cc b/src/solver/sparse_solver_mumps.cc index 9b2632a1c..435e0d4fe 100644 --- a/src/solver/sparse_solver_mumps.cc +++ b/src/solver/sparse_solver_mumps.cc @@ -1,429 +1,443 @@ /** * @file sparse_solver_mumps.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Tue Jan 19 2016 * * @brief implem of SparseSolverMumps class * * @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 * * @subsection Ctrl_param Control parameters * * ICNTL(1), * ICNTL(2), * ICNTL(3) : output streams for error, diagnostics, and global messages * * ICNTL(4) : verbose level : 0 no message - 4 all messages * * ICNTL(5) : type of matrix, 0 assembled, 1 elementary * * ICNTL(6) : control the permutation and scaling(default 7) see mumps doc for * more information * * ICNTL(7) : determine the pivot order (default 7) see mumps doc for more * information * * ICNTL(8) : describe the scaling method used * * ICNTL(9) : 1 solve A x = b, 0 solve At x = b * * ICNTL(10) : number of iterative refinement when NRHS = 1 * * ICNTL(11) : > 0 return statistics * * ICNTL(12) : only used for SYM = 2, ordering strategy * * ICNTL(13) : * * ICNTL(14) : percentage of increase of the estimated working space * * ICNTL(15-17) : not used * * ICNTL(18) : only used if ICNTL(5) = 0, 0 matrix centralized, 1 structure on * host and mumps give the mapping, 2 structure on host and distributed matrix * for facto, 3 distributed matrix * * ICNTL(19) : > 0, Shur complement returned * * ICNTL(20) : 0 rhs dense, 1 rhs sparse * * ICNTL(21) : 0 solution in rhs, 1 solution distributed in ISOL_loc and SOL_loc * allocated by user * * ICNTL(22) : 0 in-core, 1 out-of-core * * ICNTL(23) : maximum memory allocatable by mumps pre proc * * ICNTL(24) : controls the detection of "null pivot rows" * * ICNTL(25) : * * ICNTL(26) : * * ICNTL(27) : * * ICNTL(28) : 0 automatic choice, 1 sequential analysis, 2 parallel analysis * * ICNTL(29) : 0 automatic choice, 1 PT-Scotch, 2 ParMetis */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dof_manager_default.hh" #include "sparse_matrix_aij.hh" +#include "dof_synchronizer.hh" #if defined(AKANTU_USE_MPI) #include "static_communicator_mpi.hh" #include "mpi_type_wrapper.hh" #endif -#include "solver_mumps.hh" +#include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ // static std::ostream & operator <<(std::ostream & stream, const DMUMPS_STRUC_C // & _this) { // stream << "DMUMPS Data [" << std::endl; // stream << " + job : " << _this.job << std::endl; // stream << " + par : " << _this.par << std::endl; // stream << " + sym : " << _this.sym << std::endl; // stream << " + comm_fortran : " << _this.comm_fortran << std::endl; // stream << " + nz : " << _this.nz << std::endl; // stream << " + irn : " << _this.irn << std::endl; // stream << " + jcn : " << _this.jcn << std::endl; // stream << " + nz_loc : " << _this.nz_loc << std::endl; // stream << " + irn_loc : " << _this.irn_loc << std::endl; // stream << " + jcn_loc : " << _this.jcn_loc << std::endl; // stream << "]"; // return stream; // } __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseSolverMumps::SparseSolverMumps(DOFManagerDefault & dof_manager, const ID & matrix_id, const ID & id, const MemoryID & memory_id) : SparseSolver(dof_manager, matrix_id, id, memory_id), dof_manager(dof_manager), matrix(dof_manager.getMatrix(matrix_id)), master_rhs_solution(0, 1), is_initialized(false) { AKANTU_DEBUG_IN(); this->prank = communicator.whoAmI(); #ifdef AKANTU_USE_MPI this->parallel_method = _fully_distributed; #else // AKANTU_USE_MPI this->parallel_method = _not_parallel; #endif // AKANTU_USE_MPI this->mumps_data.par = 1; // The host is part of computations switch (this->parallel_method) { case _not_parallel: break; case _master_slave_distributed: this->mumps_data.par = 0; // The host is not part of the computations case _fully_distributed: #ifdef AKANTU_USE_MPI const StaticCommunicatorMPI & mpi_st_comm = dynamic_cast( communicator.getRealStaticCommunicator()); - + MPI_Comm mpi_comm = mpi_st_comm.getMPITypeWrapper().getMPICommunicator(); this->mumps_data.comm_fortran = - MPI_Comm_c2f(mpi_st_comm.getMPITypeWrapper().getMPICommunicator()); + MPI_Comm_c2f(mpi_comm); #else AKANTU_DEBUG_ERROR( "You cannot use parallel method to solve without activating MPI"); #endif break; } this->mumps_data.sym = 2 * (this->matrix.getMatrixType() == _symmetric); this->prank = communicator.whoAmI(); this->setOutputLevel(); this->mumps_data.job = _smj_initialize; // initialize dmumps_c(&this->mumps_data); this->setOutputLevel(); this->is_initialized = true; this->last_profile_release = this->matrix.getProfileRelease() - 1; this->last_value_release = this->matrix.getValueRelease() - 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseSolverMumps::~SparseSolverMumps() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::destroyInternalData() { if(this->is_initialized) { this->mumps_data.job = _smj_destroy; // destroy dmumps_c(&this->mumps_data); this->is_initialized = false; } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::checkInitialized() { if(!this->is_initialized) AKANTU_EXCEPTION("You cannot use an un/de-initialized mumps solver"); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::setOutputLevel() { // Output setup icntl(1) = 0; // error output icntl(2) = 0; // dignostics output icntl(3) = 0; // informations icntl(4) = 0; #if !defined(AKANTU_NDEBUG) DebugLevel dbg_lvl = debug::debugger.getDebugLevel(); if (AKANTU_DEBUG_TEST(dblDump)) { strcpy(this->mumps_data.write_problem, "mumps_matrix.mtx"); } icntl(1) = (dbg_lvl >= dblWarning) ? 6 : 0; icntl(3) = (dbg_lvl >= dblInfo) ? 6 : 0; icntl(2) = (dbg_lvl >= dblTrace) ? 6 : 0; if (dbg_lvl >= dblDump) { icntl(4) = 4; } else if (dbg_lvl >= dblTrace) { icntl(4) = 3; } else if (dbg_lvl >= dblInfo) { icntl(4) = 2; } else if (dbg_lvl >= dblWarning) { icntl(4) = 1; } #endif } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initMumpsData() { // Default Scaling icntl(8) = 77; // Assembled matrix icntl(5) = 0; /// Default centralized dense second member icntl(20) = 0; icntl(21) = 0; // automatic choice for analysis icntl(28) = 0; UInt size = this->matrix.getSize(); if (prank == 0) { this->master_rhs_solution.resize(size); } this->mumps_data.nz_alloc = 0; this->mumps_data.n = size; switch (this->parallel_method) { case _fully_distributed: icntl(18) = 3; // fully distributed this->mumps_data.nz_loc = matrix.getNbNonZero(); this->mumps_data.irn_loc = matrix.getIRN().storage(); this->mumps_data.jcn_loc = matrix.getJCN().storage(); break; case _not_parallel: case _master_slave_distributed: icntl(18) = 0; // centralized if (prank == 0) { this->mumps_data.nz = matrix.getNbNonZero(); this->mumps_data.irn = matrix.getIRN().storage(); this->mumps_data.jcn = matrix.getJCN().storage(); } else { this->mumps_data.nz = 0; this->mumps_data.irn = NULL; this->mumps_data.jcn = NULL; } break; default: AKANTU_DEBUG_ERROR("This case should not happen!!"); } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initialize() { AKANTU_DEBUG_IN(); this->analysis(); // icntl(14) = 80; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::analysis() { AKANTU_DEBUG_IN(); initMumpsData(); this->checkInitialized(); this->mumps_data.job = _smj_analyze; // analyze dmumps_c(&this->mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::factorize() { AKANTU_DEBUG_IN(); // should be null on slaves // this->mumps_data.rhs = this->master_rhs_solution.storage(); if (parallel_method == _fully_distributed) this->mumps_data.a_loc = this->matrix.getA().storage(); else { if (prank == 0) this->mumps_data.a = this->matrix.getA().storage(); } this->checkInitialized(); this->mumps_data.job = _smj_factorize; // factorize dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SparseSolverMumps::solve() { - AKANTU_DEBUG_IN(); - - this->setOutputLevel(); +void SparseSolverMumps::solve(Array & x, const Array & b) { + auto & synch = this->dof_manager.getSynchronizer(); + + if (this->prank == 0) { + this->master_rhs_solution.resize(this->dof_manager.getSystemSize()); + synch.gather(b, this->master_rhs_solution); + } else { + synch.gather(b); + } - this->dof_manager.applyBoundary(); + this->solveInternal(); - if (AKANTU_DEBUG_TEST(dblDump)) { - std::stringstream sstr; sstr << prank << ".mtx"; - this->matrix.saveMatrix("solver_mumps" + sstr.str()); + if (this->prank == 0) { + synch.scatter(x, this->master_rhs_solution); + } else { + synch.scatter(x); } +} +/* -------------------------------------------------------------------------- */ +void SparseSolverMumps::solve() { this->master_rhs_solution.copy(this->dof_manager.getGlobalResidual()); - // if (prank == 0) { - // matrix.getDOFSynchronizer().gather(this->rhs, 0, this->master_rhs_solution); - // } else { - // this->matrix.getDOFSynchronizer().gather(this->rhs, 0); - // } + + this->solveInternal(); + + this->dof_manager.setGlobalSolution(this->master_rhs_solution); + + this->dof_manager.splitSolutionPerDOFs(); +} + +/* -------------------------------------------------------------------------- */ +void SparseSolverMumps::solveInternal() { + AKANTU_DEBUG_IN(); + + this->setOutputLevel(); if (this->last_profile_release != this->matrix.getProfileRelease()) { this->analysis(); this->last_profile_release = this->matrix.getProfileRelease(); } + this->dof_manager.applyBoundary(this->matrix_id); + if (AKANTU_DEBUG_TEST(dblDump)) { + std::stringstream sstr; sstr << prank << ".mtx"; + this->matrix.saveMatrix("solver_mumps" + sstr.str()); + } + if (this->last_value_release != this->matrix.getValueRelease()) { this->factorize(); this->last_value_release = this->matrix.getValueRelease(); } if (prank == 0) { this->mumps_data.rhs = this->master_rhs_solution.storage(); } this->checkInitialized(); this->mumps_data.job = _smj_solve; // solve dmumps_c(&this->mumps_data); this->printError(); - this->dof_manager.setGlobalSolution(this->master_rhs_solution); - - // if (prank == 0) { - // matrix.getDOFSynchronizer().scatter(this->solution, 0, this->master_rhs_solution); - // } else { - // this->matrix.getDOFSynchronizer().gather(this->solution, 0); - // } - - this->dof_manager.splitSolutionPerDOFs(); - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::printError() { Vector _info_v(2); _info_v[0] = info(1); // to get errors _info_v[1] = -info(1); // to get warnings StaticCommunicator::getStaticCommunicator().allReduce(_info_v, _so_min); _info_v[1] = -_info_v[1]; if (_info_v[0] < 0) { // < 0 is an error switch (_info_v[0]) { case -10: AKANTU_DEBUG_ERROR("The matrix is singular"); break; case -9: { icntl(14) += 10; if (icntl(14) != 90) { // std::cout << "Dynamic memory increase of 10%" << std::endl; AKANTU_DEBUG_WARNING("MUMPS dynamic memory is insufficient it will be " "increased allowed to use 10% more"); // change releases to force a recompute this->last_value_release--; this->last_profile_release--; this->solve(); } else { AKANTU_DEBUG_ERROR("The MUMPS workarray is too small INFO(2)=" << info(2) << "No further increase possible"); break; } } default: AKANTU_DEBUG_ERROR("Error in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } else if (_info_v[1] > 0) { AKANTU_DEBUG_WARNING("Warning in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } __END_AKANTU__ diff --git a/src/solver/sparse_solver_mumps.hh b/src/solver/sparse_solver_mumps.hh index 887d99ce6..b30b0e8dc 100644 --- a/src/solver/sparse_solver_mumps.hh +++ b/src/solver/sparse_solver_mumps.hh @@ -1,156 +1,160 @@ /** * @file sparse_solver_mumps.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jan 19 2016 * * @brief Solver class implementation for the mumps solver * * @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 "sparse_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_MUMPS_HH__ #define __AKANTU_SOLVER_MUMPS_HH__ namespace akantu { class DOFManagerDefault; class SparseMatrixAIJ; } __BEGIN_AKANTU__ class SparseSolverMumps : public SparseSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseSolverMumps(DOFManagerDefault & dof_manager, const ID & matrix_id, const ID & id = "sparse_solver_mumps", const MemoryID & memory_id = 0); virtual ~SparseSolverMumps(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile and do the analysis part virtual void initialize(); - // void initializeSlave(SolverOptions & options = _solver_no_options); - /// analysis (symbolic facto + permutations) virtual void analysis(); /// factorize the matrix virtual void factorize(); /// solve the system + virtual void solve(Array & x, const Array & b); + + /// solve using residual and solution from the dof_manager virtual void solve(); private: /// print the error if any happened in mumps void printError(); + /// solve the system with master_rhs_solution as b and x + void solveInternal(); + /// set internal values; void initMumpsData(); /// set the level of verbosity of mumps based on the debug level of akantu void setOutputLevel(); protected: /// de-initialize the internal data virtual void destroyInternalData(); /// check if initialized and except if it is not the case void checkInitialized(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ private: /// access the control variable inline Int & icntl(UInt i) { return mumps_data.icntl[i - 1]; } /// access the results info inline Int & info(UInt i) { return mumps_data.info[i - 1]; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// DOFManager used by the Mumps implementation of the SparseSolver DOFManagerDefault & dof_manager; /// AIJ Matrix, usualy the jacobian matrix SparseMatrixAIJ & matrix; /// Full right hand side on the master processors and solution after solve Array master_rhs_solution; /// mumps data DMUMPS_STRUC_C mumps_data; /// Rank of the current process UInt prank; /// matrix release at last solve UInt last_profile_release; /// matrix release at last solve UInt last_value_release; /// check if the solver data are initialized bool is_initialized; /* ------------------------------------------------------------------------ */ /* Local types */ /* ------------------------------------------------------------------------ */ private: SolverParallelMethod parallel_method; // bool rhs_is_local; enum SolverMumpsJob { _smj_initialize = -1, _smj_analyze = 1, _smj_factorize = 2, _smj_solve = 3, _smj_analyze_factorize = 4, _smj_factorize_solve = 5, _smj_complete = 6, // analyze, factorize, solve _smj_destroy = -2 }; }; __END_AKANTU__ #endif /* __AKANTU_SOLVER_MUMPS_HH__ */ diff --git a/src/synchronizer/communication_buffer.hh b/src/synchronizer/communication_buffer.hh index 0ef7d3f7b..1a10d2467 100644 --- a/src/synchronizer/communication_buffer.hh +++ b/src/synchronizer/communication_buffer.hh @@ -1,192 +1,193 @@ /** * @file communication_buffer.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Sun Oct 19 2014 * * @brief Buffer for packing and unpacking data * * @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 "aka_array.hh" +#include "element.hh" #ifndef __AKANTU_COMMUNICATION_BUFFER_HH__ #define __AKANTU_COMMUNICATION_BUFFER_HH__ __BEGIN_AKANTU__ template class CommunicationBufferTemplated { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: explicit CommunicationBufferTemplated(UInt size) : buffer(size, 1) { ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); }; CommunicationBufferTemplated() : CommunicationBufferTemplated(0) {} CommunicationBufferTemplated(const CommunicationBufferTemplated & other) { buffer = other.buffer; ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); } CommunicationBufferTemplated & operator=(const CommunicationBufferTemplated & other) { if (this != &other) { buffer = other.buffer; ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); } return *this; } virtual ~CommunicationBufferTemplated() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// reset to "empty" inline void reset(); /// resize the internal buffer inline void resize(UInt size); /// clear buffer context inline void clear(); private: inline void packResize(UInt size); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: inline char * storage() { return buffer.storage(); }; /* ------------------------------------------------------------------------ */ /* Operators */ /* ------------------------------------------------------------------------ */ public: /// printing tool template inline std::string extractStream(UInt packet_size); /// packing data template inline CommunicationBufferTemplated & operator<< (const T & to_pack); template inline CommunicationBufferTemplated & operator<< (const Vector & to_pack); template inline CommunicationBufferTemplated & operator<< (const Matrix & to_pack); template inline CommunicationBufferTemplated & operator<< (const std::vector & to_pack); /// unpacking data template inline CommunicationBufferTemplated & operator>> (T & to_unpack); template inline CommunicationBufferTemplated & operator>> (Vector & to_unpack); template inline CommunicationBufferTemplated & operator>> (Matrix & to_unpack); template inline CommunicationBufferTemplated & operator>> (std::vector & to_unpack); inline CommunicationBufferTemplated & operator<< (const std::string & to_pack); inline CommunicationBufferTemplated & operator>> (std::string & to_unpack); private: template inline void packIterable (T & to_pack); template inline void unpackIterable (T & to_pack); /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: template static inline UInt sizeInBuffer(const T & data); template static inline UInt sizeInBuffer(const Vector & data); template static inline UInt sizeInBuffer(const Matrix & data); template static inline UInt sizeInBuffer(const std::vector & data); static inline UInt sizeInBuffer(const std::string & data); /// return the size in bytes of the stored values inline UInt getPackedSize(){ return ptr_pack - buffer.storage(); }; /// return the size in bytes of data left to be unpacked inline UInt getLeftToUnpack(){ return buffer.getSize() - (ptr_unpack - buffer.storage()); }; /// return the global size allocated inline UInt getSize(){ return buffer.getSize(); }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// current position for packing char * ptr_pack; /// current position for unpacking char * ptr_unpack; /// storing buffer Array buffer; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "communication_buffer_inline_impl.cc" #endif typedef CommunicationBufferTemplated CommunicationBuffer; typedef CommunicationBufferTemplated DynamicCommunicationBuffer; __END_AKANTU__ #endif /* __AKANTU_COMMUNICATION_BUFFER_HH__ */ diff --git a/src/synchronizer/communication_buffer_inline_impl.cc b/src/synchronizer/communication_buffer_inline_impl.cc index f0655f3be..73de024a3 100644 --- a/src/synchronizer/communication_buffer_inline_impl.cc +++ b/src/synchronizer/communication_buffer_inline_impl.cc @@ -1,315 +1,319 @@ /** * @file communication_buffer_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Thu Apr 14 2011 * @date last modification: Sun Oct 19 2014 * * @brief CommunicationBuffer inline implementation * * @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 . * */ /* -------------------------------------------------------------------------- */ template template inline UInt CommunicationBufferTemplated::sizeInBuffer(const T &) { return sizeof(T); } template template inline UInt CommunicationBufferTemplated::sizeInBuffer(const Vector & data) { UInt size = data.size() * sizeof(T); return size; } template template inline UInt CommunicationBufferTemplated::sizeInBuffer(const Matrix & data) { UInt size = data.size() * sizeof(T); return size; } template template inline UInt CommunicationBufferTemplated::sizeInBuffer( const std::vector & data) { UInt size = data.size() * sizeof(T) + sizeof(size_t); return size; } template inline UInt CommunicationBufferTemplated::sizeInBuffer( const std::string & data) { UInt size = data.size() * sizeof(std::string::value_type) + sizeof(size_t); return size; } /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::packResize(UInt size) { if (!is_static) { char * values = buffer.storage(); buffer.resize(buffer.getSize() + size); ptr_pack = buffer.storage() + (ptr_pack - values); ptr_unpack = buffer.storage() + (ptr_unpack - values); } } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator<<(const T & to_pack) { UInt size = sizeInBuffer(to_pack); packResize(size); - T * tmp = reinterpret_cast(ptr_pack); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.getSize()) >= (ptr_pack + size), "Packing too much data in the CommunicationBufferTemplated"); - *tmp = to_pack; + memcpy(ptr_pack, reinterpret_cast(&to_pack), size); ptr_pack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(T & to_unpack) { UInt size = sizeInBuffer(to_unpack); - T * tmp = reinterpret_cast(ptr_unpack); + + alignas(alignof(T)) std::array aligned_ptr; + memcpy(aligned_ptr.data(), ptr_unpack, size); + + T * tmp = reinterpret_cast(aligned_ptr.data()); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.getSize()) >= (ptr_unpack + size), "Unpacking too much data in the CommunicationBufferTemplated"); to_unpack = *tmp; + // memcpy(reinterpret_cast(&to_unpack), ptr_unpack, size); ptr_unpack += size; return *this; } /* -------------------------------------------------------------------------- */ /* Specialization */ /* -------------------------------------------------------------------------- */ /** * Vector */ /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator<<(const Vector & to_pack) { UInt size = sizeInBuffer(to_pack); packResize(size); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.getSize()) >= (ptr_pack + size), "Packing too much data in the CommunicationBufferTemplated"); memcpy(ptr_pack, to_pack.storage(), size); ptr_pack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(Vector & to_unpack) { UInt size = sizeInBuffer(to_unpack); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.getSize()) >= (ptr_unpack + size), "Unpacking too much data in the CommunicationBufferTemplated"); memcpy(to_unpack.storage(), ptr_unpack, size); ptr_unpack += size; return *this; } /** * Matrix */ /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator<<(const Matrix & to_pack) { UInt size = sizeInBuffer(to_pack); packResize(size); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.getSize()) >= (ptr_pack + size), "Packing too much data in the CommunicationBufferTemplated"); memcpy(ptr_pack, to_pack.storage(), size); ptr_pack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(Matrix & to_unpack) { UInt size = sizeInBuffer(to_unpack); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.getSize()) >= (ptr_unpack + size), "Unpacking too much data in the CommunicationBufferTemplated"); memcpy(to_unpack.storage(), ptr_unpack, size); ptr_unpack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline void CommunicationBufferTemplated::packIterable(T & to_pack) { operator<<(size_t(to_pack.size())); typename T::const_iterator it = to_pack.begin(); typename T::const_iterator end = to_pack.end(); for (; it != end; ++it) operator<<(*it); } /* -------------------------------------------------------------------------- */ template template inline void CommunicationBufferTemplated::unpackIterable(T & to_unpack) { size_t size; operator>>(size); to_unpack.resize(size); typename T::iterator it = to_unpack.begin(); typename T::iterator end = to_unpack.end(); for (; it != end; ++it) operator>>(*it); } /** * std::vector */ /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated:: operator<<(const std::vector & to_pack) { packIterable(to_pack); return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated:: operator>>(std::vector & to_unpack) { unpackIterable(to_unpack); return *this; } /** * std::string */ /* -------------------------------------------------------------------------- */ template inline CommunicationBufferTemplated & CommunicationBufferTemplated:: operator<<(const std::string & to_pack) { packIterable(to_pack); return *this; } /* -------------------------------------------------------------------------- */ template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(std::string & to_unpack) { unpackIterable(to_unpack); return *this; } /* -------------------------------------------------------------------------- */ template template inline std::string CommunicationBufferTemplated::extractStream(UInt block_size) { std::stringstream str; T * ptr = reinterpret_cast(buffer.storage()); UInt sz = buffer.getSize() / sizeof(T); UInt sz_block = block_size / sizeof(T); UInt n_block = 0; for (UInt i = 0; i < sz; ++i) { if (i % sz_block == 0) { str << std::endl << n_block << " "; ++n_block; } str << *ptr << " "; ++ptr; } return str.str(); } /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::resize(UInt size) { if (!is_static) { buffer.resize(0); } else { buffer.resize(size); } reset(); #ifndef AKANTU_NDEBUG clear(); #endif } /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::clear() { buffer.clear(); } /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::reset() { ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); } /* -------------------------------------------------------------------------- */ // template // inline CommunicationBufferTemplated & // CommunicationBufferTemplated::packMeshData (const MeshData & // to_pack, const ElementType & type) { // UInt size = to_pack.size(); // operator<<(size); // typename std::vector::iterator it = to_pack.begin(); // typename std::vector::iterator end = to_pack.end(); // for(;it != end; ++it) operator<<(*it); // return *this; //} diff --git a/src/synchronizer/dof_synchronizer.cc b/src/synchronizer/dof_synchronizer.cc index 7e17e8477..2bc478a4c 100644 --- a/src/synchronizer/dof_synchronizer.cc +++ b/src/synchronizer/dof_synchronizer.cc @@ -1,216 +1,220 @@ /** * @file dof_synchronizer.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Wed Oct 21 2015 * * @brief DOF synchronizing object implementation * * @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 "dof_synchronizer.hh" #include "dof_manager_default.hh" #include "mesh.hh" #include "node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /** * A DOFSynchronizer needs a mesh and the number of degrees of freedom * per node to be created. In the constructor computes the local and global dof * number for each dof. The member * proc_informations (std vector) is resized with the number of mpi * processes. Each entry in the vector is a PerProcInformations object * that contains the interactions of the current mpi process (prank) with the * mpi process corresponding to the position of that entry. Every * ProcInformations object contains one array with the dofs that have * to be sent to prank and a second one with dofs that willl be received form * prank. * This information is needed for the asychronous communications. The * constructor sets up this information. */ DOFSynchronizer::DOFSynchronizer(DOFManagerDefault & dof_manager, const ID & id, - MemoryID memory_id, const StaticCommunicator & comm) + MemoryID memory_id, + const StaticCommunicator & comm) : SynchronizerImpl(id, memory_id, comm), root(0), dof_manager(dof_manager), - slave_receive_dofs(0, 1, "dofs-to-receive-from-master"), + root_dofs(0, 1, "dofs-to-receive-from-master"), dof_changed(true) { std::vector dof_ids = dof_manager.getDOFIDs(); // Transfers nodes to global equation numbers in new schemes for (ID dof_id : dof_ids) { registerDOFs(dof_id); } this->initScatterGatherCommunicationScheme(); } /* -------------------------------------------------------------------------- */ DOFSynchronizer::~DOFSynchronizer() {} /* -------------------------------------------------------------------------- */ void DOFSynchronizer::registerDOFs(const ID & dof_id) { - if(this->nb_proc == 1) return; + if (this->nb_proc == 1) + return; typedef Communications::const_scheme_iterator const_scheme_iterator; const Array equation_numbers = dof_manager.getLocalEquationNumbers(dof_id); if (dof_manager.getSupportType(dof_id) == _dst_nodal) { const NodeSynchronizer & node_synchronizer = dof_manager.getMesh().getNodeSynchronizer(); const Array & associated_nodes = dof_manager.getDOFsAssociatedNodes(dof_id); const Communications & node_communications = node_synchronizer.getCommunications(); auto transcode_node_to_global_dof_scheme = [this, &associated_nodes, &equation_numbers](const_scheme_iterator it, const_scheme_iterator end, const CommunicationSendRecv & sr) -> void { for (; it != end; ++it) { auto & scheme = communications.createScheme(it->first, sr); const auto & node_scheme = it->second; for (auto & node : node_scheme) { auto an_begin = associated_nodes.begin(); auto an_it = an_begin; auto an_end = associated_nodes.end(); std::vector global_dofs_per_node; while ((an_it = std::find(an_it, an_end, node)) != an_end) { UInt pos = an_it - an_begin; UInt local_eq_num = equation_numbers(pos); - UInt global_eq_num = local_eq_num; - dof_manager.localToGlobalEquationNumber(global_eq_num); + UInt global_eq_num = + dof_manager.localToGlobalEquationNumber(local_eq_num); global_dofs_per_node.push_back(global_eq_num); ++an_it; } std::sort(global_dofs_per_node.begin(), global_dofs_per_node.end()); - std::transform( - global_dofs_per_node.begin(), global_dofs_per_node.end(), - global_dofs_per_node.begin(), [this](UInt g) -> UInt { - UInt l = dof_manager.globalToLocalEquationNumber(g); - return l; - }); + std::transform(global_dofs_per_node.begin(), + global_dofs_per_node.end(), + global_dofs_per_node.begin(), [this](UInt g) -> UInt { + UInt l = dof_manager.globalToLocalEquationNumber(g); + return l; + }); for (auto & leqnum : global_dofs_per_node) { scheme.push_back(leqnum); } } } }; auto nss_it = node_communications.begin_send_scheme(); auto nss_end = node_communications.end_send_scheme(); transcode_node_to_global_dof_scheme(nss_it, nss_end, _send); auto nrs_it = node_communications.begin_recv_scheme(); auto nrs_end = node_communications.end_recv_scheme(); transcode_node_to_global_dof_scheme(nrs_it, nrs_end, _recv); } dof_changed = true; } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::initScatterGatherCommunicationScheme() { AKANTU_DEBUG_IN(); if (this->nb_proc == 1) { + dof_changed = false; AKANTU_DEBUG_OUT(); return; } UInt nb_dofs = dof_manager.getLocalSystemSize(); - this->slave_receive_dofs.clear(); + this->root_dofs.clear(); this->master_receive_dofs.clear(); Array dofs_to_send; for (UInt n = 0; n < nb_dofs; ++n) { if (dof_manager.isLocalOrMasterDOF(n)) { - this->slave_receive_dofs.push_back(n); + auto & receive_per_proc = master_receive_dofs[this->root]; + UInt global_dof = dof_manager.localToGlobalEquationNumber(n); - UInt global_dof = n; - dof_manager.localToGlobalEquationNumber(global_dof); + root_dofs.push_back(n); + receive_per_proc.push_back(global_dof); dofs_to_send.push_back(global_dof); } } if (this->rank == UInt(this->root)) { Array nb_dof_per_proc(this->nb_proc); communicator.gather(dofs_to_send.getSize(), nb_dof_per_proc); std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == UInt(this->root)) continue; - Array & receive_per_proc = master_receive_dofs[p]; + auto & receive_per_proc = master_receive_dofs[p]; receive_per_proc.resize(nb_dof_per_proc(p)); - if(nb_dof_per_proc(p) != 0) + if (nb_dof_per_proc(p) != 0) requests.push_back(communicator.asyncReceive( - receive_per_proc, p, - Tag::genTag(p, 0, Tag::_GATHER_INITIALIZATION, this->hash_id))); + receive_per_proc, p, + Tag::genTag(p, 0, Tag::_GATHER_INITIALIZATION, this->hash_id))); } communicator.waitAll(requests); communicator.freeCommunicationRequest(requests); } else { communicator.gather(dofs_to_send.getSize(), this->root); AKANTU_DEBUG(dblDebug, "I have " << nb_dofs << " dofs (" << dofs_to_send.getSize() << " to send to master proc"); - if(dofs_to_send.getSize() != 0) - communicator.send( - dofs_to_send, this->root, - Tag::genTag(this->rank, 0, Tag::_GATHER_INITIALIZATION, this->hash_id)); + if (dofs_to_send.getSize() != 0) + communicator.send(dofs_to_send, this->root, + Tag::genTag(this->rank, 0, Tag::_GATHER_INITIALIZATION, + this->hash_id)); } dof_changed = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool DOFSynchronizer::hasChanged() { communicator.allReduce(dof_changed, _so_lor); return dof_changed; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/synchronizer/dof_synchronizer.hh b/src/synchronizer/dof_synchronizer.hh index 179b98f0a..624b125c4 100644 --- a/src/synchronizer/dof_synchronizer.hh +++ b/src/synchronizer/dof_synchronizer.hh @@ -1,110 +1,110 @@ /** * @file dof_synchronizer.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Tue Dec 08 2015 * * @brief Synchronize Array of DOFs * * @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_array.hh" #include "aka_common.hh" #include "synchronizer_impl.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Mesh; class DOFManagerDefault; } #ifndef __AKANTU_DOF_SYNCHRONIZER_HH__ #define __AKANTU_DOF_SYNCHRONIZER_HH__ namespace akantu { class DOFSynchronizer : public SynchronizerImpl { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFSynchronizer( DOFManagerDefault & dof_manager, const ID & id = "dof_synchronizer", MemoryID memory_id = 0, const StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator()); virtual ~DOFSynchronizer(); virtual void registerDOFs(const ID & dof_id); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template /// Gather the DOF value on the root proccessor void gather(const Array & to_gather, Array & gathered); /// Gather the DOF value on the root proccessor template void gather(const Array & to_gather); /// Scatter a DOF Array form root to all processors template void scatter(Array & scattered, const Array & to_scatter); /// Scatter a DOF Array form root to all processors template void scatter(Array & scattered); template void synchronize(Array & vector) const; template