diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index d93aa8087..10e138763 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,563 +1,563 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ using MemoryID = UInt; // using Surface = std::string; // using SurfacePair= std::pair; // using SurfacePairList = std::list; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; #define AKANTU_PP_ENUM(s, data, i, elem) \ BOOST_PP_TUPLE_REM() \ elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data))) } // namespace akantu #if (defined(__GNUC__) || defined(__GNUG__)) #define AKA_GCC_VERSION \ (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) #if AKA_GCC_VERSION < 60000 #define AKANTU_ENUM_HASH(type_name) \ namespace std { \ template <> struct hash<::akantu::type_name> { \ using argument_type = ::akantu::type_name; \ size_t operator()(const argument_type & e) const noexcept { \ auto ue = underlying_type_t(e); \ return uh(ue); \ } \ \ private: \ const hash> uh{}; \ }; \ } #else #define AKANTU_ENUM_HASH(type_name) #endif // AKA_GCC_VERSION #endif // GNU #include "aka_element_classes_info.hh" namespace akantu { #define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) #define AKANTU_PP_TYPE_TO_STR(s, data, elem) \ ({BOOST_PP_CAT(data::_, elem), BOOST_PP_STRINGIZE(elem)}) #define AKANTU_PP_STR_TO_TYPE(s, data, elem) \ ({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data::_, elem)}) #define AKANTU_ENUM_DECLARE(type_name, list) \ enum class type_name { \ BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \ }; #define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ } \ AKANTU_ENUM_HASH(type_name) \ namespace aka { \ inline std::string to_string(const ::akantu::type_name & type) { \ static std::unordered_map<::akantu::type_name, std::string> convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, \ ::akantu::type_name, list))}; \ return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << aka::to_string(type); \ return stream; \ } #define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ inline std::istream & operator>>(std::istream & stream, type_name & type) { \ std::string str; \ stream >> str; \ static std::unordered_map convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, type_name, list))}; \ type = convert.at(str); \ return stream; \ } /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #ifndef SWIG // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { _model, _solid_mechanics_model, _solid_mechanics_model_cohesive, _heat_transfer_model, _structural_mechanics_model, _embedded_model }; #endif /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent _nls_lumped, ///< Case of lumped mass or equivalent matrix _nls_auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_ask_nodes, _gst_size, //--- SolidMechanicsModel tags --- _gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _gst_smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _gst_smm_uv, ///< synchronization of the nodal velocities and displacement _gst_smm_res, ///< synchronization of the nodal residual _gst_smm_init_mat, ///< synchronization of the data to initialize materials _gst_smm_stress, ///< synchronization of the stresses to compute the internal /// forces _gst_smmc_facets, ///< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, ///< synchronization of facet global connectivity _gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet /// synch _gst_smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _gst_giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _gst_ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gst_gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _gst_htm_temperature, ///< synchronization of the nodal temperature _gst_htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _gst_htm_phi, ///< synchronization of the nodal level set value phi _gst_htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _gst_mnl_for_average, ///< synchronization of data to average in non local /// material _gst_mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _gst_nh_criterion, // --- General tags --- _gst_test, ///< Test tag _gst_user_1, ///< tag for user simulations _gst_user_2, ///< tag for user simulations _gst_material_id, ///< synchronization of the material ids _gst_for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _gst_cf_nodal, ///< synchronization of disp, velo, and current position _gst_cf_incr, ///< synchronization of increment // --- Solver tags --- _gst_solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, - _local_master_mask = 0x22, //~(_master & _periodic_mask) + _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } inline std::ostream & operator<<(std::ostream & stream, const NodeFlag & flag) { using under = std::underlying_type_t; stream << under(flag); return stream; } } // namespace akantu #ifndef SWIG AKANTU_ENUM_HASH(GhostType) #endif namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ struct TensorTrait {}; /* -------------------------------------------------------------------------- */ template using is_tensor = std::is_base_of; /* -------------------------------------------------------------------------- */ template using is_scalar = std::is_arithmetic; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ // #if defined(__INTEL_COMPILER) // /// remark #981: operands are evaluated in unspecified order // #pragma warning(disable : 981) // /// remark #383: value copied to temporary, reference to temporary used // #pragma warning(disable : 383) // #endif // defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ } // namespace akantu #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic number * is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/model/dof_manager.cc b/src/model/dof_manager.cc index 795e3182b..7c43b65af 100644 --- a/src/model/dof_manager.cc +++ b/src/model/dof_manager.cc @@ -1,604 +1,605 @@ /** * @file dof_manager.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the common parts of the DOFManagers * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager.hh" #include "communicator.hh" #include "element_group.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), communicator(Communicator::getStaticCommunicator()) {} /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(Mesh & mesh, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), mesh(&mesh), local_system_size(0), pure_local_system_size(0), system_size(0), communicator(mesh.getCommunicator()) { this->mesh->registerEventHandler(*this, _ehp_dof_manager); } /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() = default; /* -------------------------------------------------------------------------- */ // void DOFManager::getEquationsNumbers(const ID &, Array &) { // AKANTU_TO_IMPLEMENT(); // } /* -------------------------------------------------------------------------- */ std::vector DOFManager::getDOFIDs() const { std::vector keys; for (const auto & dof_data : this->dofs) keys.push_back(dof_data.first); return keys; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_element; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_vect.size() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); Array::const_matrix_iterator elem_it = elementary_vect.begin(nb_degree_of_freedom, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++elem_it) { UInt element = el; if (filter_it != nullptr) { // conn_it = conn_begin + *filter_it; element = *filter_it; } // const Vector & conn = *conn_it; const Matrix & elemental_val = *elem_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_node = connectivity(element, n) * nb_degree_of_freedom; Vector assemble(array_assembeled.storage() + offset_node, nb_degree_of_freedom); Vector elem_val = elemental_val(n); assemble.aXplusY(elem_val, scale_factor); } if (filter_it != nullptr) ++filter_it; // else // ++conn_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.clear(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToResidual(dof_id, array_localy_assembeled, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.clear(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToLumpedMatrix(dof_id, array_localy_assembeled, lumped_mtx, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = *pair.second; this->assembleMatMulVectToResidual(dof_id, A_id, *dof_data.dof, scale_factor); } } /* -------------------------------------------------------------------------- */ DOFManager::DOFData::DOFData(const ID & dof_id) : support_type(_dst_generic), group_support("__mesh__"), dof(nullptr), blocked_dofs(nullptr), increment(nullptr), previous(nullptr), solution(0, 1, dof_id + ":solution"), local_equation_number(0, 1, dof_id + ":local_equation_number") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData::~DOFData() = default; /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManager::getNewDOFData(const ID & dof_id) { auto it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } std::unique_ptr dofs_storage = std::make_unique(dof_id); this->dofs[dof_id] = std::move(dofs_storage); return *dofs_storage; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { DOFData & dofs_storage = this->getDOFData(dof_id); dofs_storage.dof = &dofs_array; UInt nb_local_dofs = 0; UInt nb_pure_local = 0; const DOFSupportType & support_type = dofs_storage.support_type; switch (support_type) { case _dst_nodal: { UInt nb_nodes = 0; const ID & group = dofs_storage.group_support; NodeGroup * node_group = nullptr; if (group == "__mesh__") { nb_nodes = this->mesh->getNbNodes(); } else { node_group = &this->mesh->getElementGroup(group).getNodeGroup(); nb_nodes = node_group->size(); } nb_local_dofs = nb_nodes; AKANTU_DEBUG_ASSERT( dofs_array.size() == nb_local_dofs, "The array of dof is too shot to be associated to nodes."); for (UInt n = 0; n < nb_nodes; ++n) { UInt node = n; if (node_group) node = node_group->getNodes()(n); nb_pure_local += this->mesh->isLocalOrMasterNode(node) ? 1 : 0; + nb_local_dofs -= this->mesh->isPeriodicSlave(node) ? 1 : 0; } nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); break; } case _dst_generic: { nb_local_dofs = nb_pure_local = dofs_array.size() * dofs_array.getNbComponent(); break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; communicator.allReduce(nb_pure_local, SynchronizerOperation::_sum); this->system_size += nb_pure_local; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { DOFData & dofs_storage = this->getNewDOFData(dof_id); dofs_storage.support_type = support_type; this->registerDOFsInternal(dof_id, dofs_array); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group) { DOFData & dofs_storage = this->getNewDOFData(dof_id); dofs_storage.support_type = _dst_nodal; dofs_storage.group_support = support_group; this->registerDOFsInternal(dof_id, dofs_array); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsPrevious(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.previous != nullptr) { AKANTU_EXCEPTION("The previous dofs array for " << dof_id << " has already been registered"); } dof.previous = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsIncrement(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.increment != nullptr) { AKANTU_EXCEPTION("The dofs increment array for " << dof_id << " has already been registered"); } dof.increment = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative) { DOFData & dof = this->getDOFData(dof_id); std::vector *> & derivatives = dof.dof_derivatives; if (derivatives.size() < order) { derivatives.resize(order, nullptr); } else { if (derivatives[order - 1] != nullptr) { AKANTU_EXCEPTION("The dof derivatives of order " << order << " already been registered for this dof (" << dof_id << ")"); } } derivatives[order - 1] = &dofs_derivative; } /* -------------------------------------------------------------------------- */ void DOFManager::registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs) { DOFData & dof = this->getDOFData(dof_id); if (dof.blocked_dofs != nullptr) { AKANTU_EXCEPTION("The blocked dofs array for " << dof_id << " has already been registered"); } dof.blocked_dofs = &blocked_dofs; } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { auto it = this->dofs.begin(); auto end = this->dofs.end(); for (; it != end; ++it) { DOFData & dof_data = *it->second; dof_data.solution.resize(dof_data.dof->size() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(it->first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix) { SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } SparseMatrix & ret = *matrix; this->matrices[matrix_id] = std::move(matrix); return ret; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix Array & DOFManager::getNewLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it != this->lumped_matrices.end()) { AKANTU_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } auto mtx = std::make_unique>(this->local_system_size, 1, matrix_id); this->lumped_matrices[matrix_id] = std::move(mtx); return *this->lumped_matrices[matrix_id]; } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::registerNonLinearSolver( const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver) { NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it != this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " already exists in " << this->id); } NonLinearSolver & ret = *non_linear_solver; this->non_linear_solvers[non_linear_solver_id] = std::move(non_linear_solver); return ret; } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::registerTimeStepSolver( const ID & time_step_solver_id, std::unique_ptr & time_step_solver) { TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it != this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " already exists in " << this->id); } TimeStepSolver & ret = *time_step_solver; this->time_step_solvers[time_step_solver_id] = std::move(time_step_solver); return ret; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::getMatrix(const ID & id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it == this->matrices.end()) { AKANTU_SILENT_EXCEPTION("The matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasMatrix(const ID & id) const { ID mtx_id = this->id + ":mtx:" + id; auto it = this->matrices.find(mtx_id); return it != this->matrices.end(); } /* -------------------------------------------------------------------------- */ Array & DOFManager::getLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const Array & DOFManager::getLumpedMatrix(const ID & id) const { ID matrix_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasLumpedMatrix(const ID & id) const { ID mtx_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(mtx_id); return it != this->lumped_matrices.end(); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::getNonLinearSolver(const ID & id) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it == this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasNonLinearSolver(const ID & id) const { ID solver_id = this->id + ":nls:" + id; auto it = this->non_linear_solvers.find(solver_id); return it != this->non_linear_solvers.end(); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::getTimeStepSolver(const ID & id) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it == this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasTimeStepSolver(const ID & solver_id) const { ID time_step_solver_id = this->id + ":tss:" + solver_id; auto it = this->time_step_solvers.find(time_step_solver_id); return it != this->time_step_solvers.end(); } /* -------------------------------------------------------------------------- */ void DOFManager::savePreviousDOFs(const ID & dofs_id) { this->getPreviousDOFs(dofs_id).copy(this->getDOFs(dofs_id)); } /* -------------------------------------------------------------------------- */ /* Mesh Events */ /* -------------------------------------------------------------------------- */ std::pair DOFManager::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { auto & dof_data = this->getDOFData(dof_id); UInt nb_new_local_dofs = 0; UInt nb_new_pure_local = 0; nb_new_local_dofs = nodes_list.size(); for (const auto & node : nodes_list) { nb_new_pure_local += this->mesh->isLocalOrMasterNode(node) ? 1 : 0; } const auto & dof_array = *dof_data.dof; nb_new_pure_local *= dof_array.getNbComponent(); nb_new_local_dofs *= dof_array.getNbComponent(); this->pure_local_system_size += nb_new_pure_local; this->local_system_size += nb_new_local_dofs; UInt nb_new_global = nb_new_pure_local; communicator.allReduce(nb_new_global, SynchronizerOperation::_sum); this->system_size += nb_new_global; dof_data.solution.resize(dof_data.solution.size() + nb_new_local_dofs); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesAdded(const Array & nodes_list, const NewNodesEvent &) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) continue; const auto & group = dof_data.group_support; if (group == "__mesh__") { this->updateNodalDOFs(dof_id, nodes_list); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup(); Array new_nodes_list; for (const auto & node : nodes_list) { if (node_group.find(node) != UInt(-1)) new_nodes_list.push_back(node); } this->updateNodalDOFs(dof_id, new_nodes_list); } } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array &, const NewElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved(const Array &, const ElementTypeMapArray &, const RemovedElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) {} /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 9ad40d7d7..e1b7ad401 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,953 +1,980 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Thu Feb 08 2018 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager_default.hh" #include "communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver_default.hh" #include "periodic_node_synchronizer.hh" #include "sparse_matrix_aij.hh" #include "terms_to_assemble.hh" #include "time_step_solver_default.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace 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_flag(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), global_equation_number(0, 1, "global_equation_number"), 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_flag(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), global_equation_number(0, 1, "global_equation_number"), first_global_dof_id(0), synchronizer(nullptr) { if (this->mesh->isDistributed()) this->synchronizer = std::make_unique( *this, this->id + ":dof_synchronizer", this->memory_id); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() = default; /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::makeConsistentForPeriodicity(const ID & dof_id, Array & array) { auto & dof_data = this->getDOFDataTyped(dof_id); if (dof_data.support_type != _dst_nodal) return; if (not mesh->isPeriodic()) return; this->mesh->getPeriodicNodeSynchronizer() .reduceSynchronizeWithPBCSlaves(array); } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); auto & dof_data = this->getDOFDataTyped(dof_id); AKANTU_DEBUG_ASSERT(dof_data.local_equation_number.size() == array_to_assemble.size() * array_to_assemble.getNbComponent(), "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); if (dof_data.support_type == _dst_nodal and mesh->isPeriodic()) { for (auto && data : zip(dof_data.local_equation_number, dof_data.associated_nodes, make_view(array_to_assemble))) { auto && equ_num = std::get<0>(data); auto && node = std::get<1>(data); auto && arr = std::get<2>(data); if (not this->mesh->isPeriodicSlave(node)) { global_array(equ_num) += scale_factor * (arr); } } } else { for (auto && data : zip(dof_data.local_equation_number, make_view(array_to_assemble))) { auto && equ_num = std::get<0>(data); auto && arr = std::get<1>(data); global_array(equ_num) += scale_factor * (arr); } } 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) { this->dofs[dof_id] = std::make_unique(dof_id); return *this->dofs[dof_id]; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs) { // auto prank = this->communicator.whoAmI(); // auto psize = this->communicator.getNbProc(); // access the relevant data to update auto & dof_data = this->getDOFDataTyped(dof_id); const auto & support_type = dof_data.support_type; const auto & group = dof_data.group_support; - if (support_type == _dst_nodal and group != "__mesh__") { - auto & support_nodes = - this->mesh->getElementGroup(group).getNodeGroup().getNodes(); - this->updateDOFsData( - dof_data, nb_dofs, nb_pure_local_dofs, - [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); - } else { - - this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs, - [](UInt node) -> UInt { return node; }); + switch(support_type) { + case _dst_nodal: + if (group != "__mesh__") { + auto & support_nodes = + this->mesh->getElementGroup(group).getNodeGroup().getNodes(); + this->updateDOFsData( + dof_data, nb_dofs, nb_pure_local_dofs, support_nodes.size(), + [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); + } else { + this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs, mesh->getNbNodes(), + [](UInt node) -> UInt { return node; }); + } + break; + case _dst_generic: + this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs); + break; } // 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; std::unique_ptr sm = std::make_unique(*this, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, 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); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, 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; std::unique_ptr nls; switch (type) { #if defined(AKANTU_IMPLICIT) case _nls_newton_raphson: case _nls_newton_raphson_modified: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } case _nls_linear: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } #endif case _nls_lumped: { nls = std::make_unique( *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"); } return this->registerNonLinearSolver(non_linear_solver_id, nls); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique( *this, type, non_linear_solver, time_step_solver_id, this->memory_id); return this->registerTimeStepSolver(time_step_solver_id, 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() { auto it = this->dofs.begin(); auto 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.size(); 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::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, Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); this->makeConsistentForPeriodicity(dof_id, array_to_assemble); this->assembleToGlobalArray(dof_id, array_to_assemble, this->residual, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); this->makeConsistentForPeriodicity(dof_id, array_to_assemble); 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.); Array tmp_residual(this->residual.size(), 1, 0.); A.matVecMul(data_cache, tmp_residual, scale_factor, 1.); this->residual += tmp_residual; } /* -------------------------------------------------------------------------- */ 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); auto A_it = A.begin(); auto A_end = A.end(); auto x_it = data_cache.begin(); auto 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(); auto & dof_data = this->getDOFData(dof_id); AKANTU_DEBUG_ASSERT(dof_data.support_type == _dst_nodal, "This function applies only on Nodal dofs"); this->addToProfile(matrix_id, dof_id, type, ghost_type); const auto & equation_number = this->getLocalEquationNumbers(dof_id); auto & A = this->getMatrix(matrix_id); UInt nb_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { if (dof_data.group_support != "__mesh__") { const auto & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); nb_element = group_elements.size(); filter_it = group_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } } AKANTU_DEBUG_ASSERT(elementary_mat.size() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = dof_data.dof->getNbComponent(); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); auto conn_begin = connectivity.begin(nb_nodes_per_element); auto conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); std::transform(element_eq_nb.begin(), element_eq_nb.end(), element_eq_nb.begin(), [&](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.size()); else this->addUnsymmetricElementalMatrixToSymmetric(A, *el_mat_it, element_eq_nb, A.size()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb, A.size()); } 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) { UInt gi = this->localToGlobalEquationNumber(equation_number_m(term.i())); UInt gj = this->localToGlobalEquationNumber(equation_number_n(term.j())); A.add(gi, gj, term); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const DOFData & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) return; auto mat_dof = std::make_pair(matrix_id, dof_id); auto type_pair = std::make_pair(type, ghost_type); auto prof_it = this->matrix_profiled_dofs.find(mat_dof); if (prof_it != this->matrix_profiled_dofs.end() && std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) != prof_it->second.end()) return; UInt nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); const auto & equation_number = this->getLocalEquationNumbers(dof_id); auto & A = this->getMatrix(matrix_id); auto size = A.size(); auto 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.size(); UInt * ge_it = nullptr; if (dof_data.group_support != "__mesh__") { const auto & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); ge_it = group_elements.storage(); nb_elements = group_elements.size(); } UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; Vector element_eq_nb(size_mat); for (UInt e = 0; e < nb_elements; ++e) { if (ge_it) cit = cbegin + *ge_it; this->extractElementEquationNumber( equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (ge_it) ++ge_it; else ++cit; for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = element_eq_nb(i); if (c_irn < size) { for (UInt j = 0; j < size_mat; ++j) { UInt c_jcn = element_eq_nb(j); if (c_jcn < size) { A.add(c_irn, c_jcn); } } } } } this->matrix_profiled_dofs[mat_dof].push_back(type_pair); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::applyBoundary(const ID & matrix_id) { SparseMatrixAIJ & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getValueRelease()) { auto are_equal = std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (not are_equal) J.applyBoundary(); previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getValueRelease(); } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getGlobalResidual() { if (this->synchronizer) { if (not this->global_residual) { this->global_residual = std::make_unique>(0, 1, "global_residual"); } if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->global_residual->resize(this->system_size); this->synchronizer->gather(this->residual, *this->global_residual); } else { this->synchronizer->gather(this->residual); } return *this->global_residual; } else { return this->residual; } } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getResidual() const { return this->residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::setGlobalSolution(const Array & solution) { if (this->synchronizer) { if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->synchronizer->scatter(this->global_solution, solution); } else { this->synchronizer->scatter(this->global_solution); } } else { AKANTU_DEBUG_ASSERT(solution.size() == this->global_solution.size(), "Sequential call to this function needs the solution " "to be the same size as the global_solution"); this->global_solution.copy(solution); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { DOFManager::onNodesAdded(nodes_list, event); if (this->synchronizer) this->synchronizer->onNodesAdded(nodes_list); } /* -------------------------------------------------------------------------- */ std::pair DOFManagerDefault::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { UInt nb_new_local_dofs, nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = DOFManager::updateNodalDOFs(dof_id, nodes_list); auto & dof_data = this->getDOFDataTyped(dof_id); - updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, + updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, nodes_list.size(), [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor(DOFManagerDefault::DOFDataDefault & dof_data, DOFManagerDefault & dof_manager) : dof_data(dof_data), dof_manager(dof_manager) { for (auto && pair : zip(dof_data.local_equation_number, dof_data.associated_nodes)) { UInt node, dof; std::tie(dof, node) = pair; dofs_per_node[node].push_back(dof); } } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_ask_nodes) { return nodes.size() * dof_data.dof->getNbComponent() * sizeof(UInt); } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_ask_nodes) { for (auto & node : nodes) { auto & dofs = dofs_per_node.at(node); for (auto & dof : dofs) { buffer << dof_manager.global_equation_number(dof); } } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { if (tag == _gst_ask_nodes) { for (auto & node : nodes) { auto & dofs = dofs_per_node[node]; for (auto dof : dofs) { UInt global_dof; buffer >> global_dof; dof_manager.global_equation_number(dof) = global_dof; dof_manager.global_to_local_mapping[global_dof] = dof; } } } } protected: std::unordered_map> dofs_per_node; DOFManagerDefault::DOFDataDefault & dof_data; DOFManagerDefault & dof_manager; }; /* -------------------------------------------------------------------------- */ -void DOFManagerDefault::updateDOFsData( - DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, - const std::function & getNode) { - auto prank = this->communicator.whoAmI(); - auto psize = this->communicator.getNbProc(); - - // access the relevant data to update - const auto & support_type = dof_data.support_type; - +void DOFManagerDefault::resizeGlobalArrays() { // resize all relevant arrays this->residual.resize(this->local_system_size, 0.); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_solution.resize(this->local_system_size, 0.); this->global_blocked_dofs.resize(this->local_system_size, true); this->previous_global_blocked_dofs.resize(this->local_system_size, true); this->global_equation_number.resize(this->local_system_size, -1); for (auto & lumped_matrix : lumped_matrices) lumped_matrix.second->resize(this->local_system_size); - auto first_dof_pos = dof_data.local_equation_number.size(); - dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + - nb_new_local_dofs); + matrix_profiled_dofs.clear(); + for (auto & matrix : matrices) { + matrix.second->clearProfile(); + } +} + +/* -------------------------------------------------------------------------- */ +auto DOFManagerDefault::computeFirstDOFIDs(UInt nb_new_local_dofs, + UInt nb_new_pure_local) { + auto prank = this->communicator.whoAmI(); + auto psize = this->communicator.getNbProc(); // determine the first local/global dof id to use Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = nb_new_pure_local; this->communicator.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_new_local_dofs; - if (support_type == _dst_nodal) { - dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + - nb_new_local_dofs); - } + auto first_global_dof_id = this->first_global_dof_id; + auto first_local_dof_id = this->local_system_size - nb_new_local_dofs; + + this->first_global_dof_id += std::accumulate( + nb_dofs_per_proc.begin() + prank, nb_dofs_per_proc.end(), 0); + + return std::make_pair(first_local_dof_id, first_global_dof_id); +} +/* -------------------------------------------------------------------------- */ +void DOFManagerDefault::updateDOFsData( + DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, + UInt nb_node, const std::function & getNode) { + auto nb_local_dofs_added = nb_node * dof_data.dof->getNbComponent(); + + resizeGlobalArrays(); + + auto first_dof_pos = dof_data.local_equation_number.size(); + dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + + nb_local_dofs_added); + dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + + nb_local_dofs_added); std::unordered_map, UInt> masters_dofs; // update per dof info - auto local_eq_num = first_dof_id; - for (auto d : arange(nb_new_local_dofs)) { - auto is_periodic_slave = false; - auto is_periodic_master = false; - auto is_local_dof = true; - auto dof_flag = NodeFlag::_normal; - - // determine the dof type - switch (support_type) { - case _dst_nodal: { - auto node = getNode(d / dof_data.dof->getNbComponent()); - dof_flag = this->mesh->getNodeFlag(node); - - dof_data.associated_nodes.push_back(node); - is_local_dof = this->mesh->isLocalOrMasterNode(node); - is_periodic_slave = this->mesh->isPeriodicSlave(node); - is_periodic_master = this->mesh->isPeriodicMaster(node); - break; - } - case _dst_generic: { - break; - } - default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } - } + UInt local_eq_num, first_global_dof_id; + std::tie(local_eq_num, first_global_dof_id) = + computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); + for (auto d : arange(nb_local_dofs_added)) { + auto node = getNode(d / dof_data.dof->getNbComponent()); + auto dof_flag = this->mesh->getNodeFlag(node); + + dof_data.associated_nodes.push_back(node); + auto is_local_dof = this->mesh->isLocalOrMasterNode(node); + auto is_periodic_slave = this->mesh->isPeriodicSlave(node); + auto is_periodic_master = this->mesh->isPeriodicMaster(node); if (is_periodic_slave) { dof_data.local_equation_number.push_back(UInt(-1)); continue; } // update equation numbers this->dofs_flag(local_eq_num) = dof_flag; dof_data.local_equation_number.push_back(local_eq_num); if (is_local_dof) { - this->global_equation_number(local_eq_num) = this->first_global_dof_id; - this->global_to_local_mapping[this->first_global_dof_id] = local_eq_num; - ++this->first_global_dof_id; + this->global_equation_number(local_eq_num) = first_global_dof_id; + this->global_to_local_mapping[first_global_dof_id] = local_eq_num; + ++first_global_dof_id; } else { this->global_equation_number(local_eq_num) = UInt(-1); } if (is_periodic_master) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof = d % dof_data.dof->getNbComponent(); masters_dofs.insert( std::make_pair(std::make_pair(node, dof), local_eq_num)); } ++local_eq_num; } // correct periodic slave equation numbers - if (support_type == _dst_nodal and this->mesh->isPeriodic()) { + if (this->mesh->isPeriodic()) { auto assoc_begin = dof_data.associated_nodes.begin(); - for (auto d : arange(nb_new_local_dofs)) { + for (auto d : arange(nb_local_dofs_added)) { auto node = dof_data.associated_nodes(first_dof_pos + d); if (not this->mesh->isPeriodicSlave(node)) continue; auto master_node = this->mesh->getPeriodicMaster(node); auto dof = d % dof_data.dof->getNbComponent(); dof_data.local_equation_number(first_dof_pos + d) = masters_dofs[std::make_pair(master_node, dof)]; } } // synchronize the global numbering for slaves - if (support_type == _dst_nodal && this->synchronizer) { + if (this->synchronizer) { GlobalDOFInfoDataAccessor data_accessor(dof_data, *this); auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, _gst_ask_nodes); } +} - this->first_global_dof_id += std::accumulate( - nb_dofs_per_proc.begin() + prank + 1, nb_dofs_per_proc.end(), 0); +/* -------------------------------------------------------------------------- */ +void DOFManagerDefault::updateDOFsData(DOFDataDefault & dof_data, + UInt nb_new_local_dofs, + UInt nb_new_pure_local) { + resizeGlobalArrays(); - matrix_profiled_dofs.clear(); - for (auto & matrix : matrices) { - matrix.second->clearProfile(); + dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + + nb_new_local_dofs); + + UInt first_local_dof_id, first_global_dof_id; + std::tie(first_local_dof_id, first_global_dof_id) = + computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); + + // update per dof info + for (auto _ [[gnu::unused]] : arange(nb_new_local_dofs)) { + // update equation numbers + this->dofs_flag(first_local_dof_id) = NodeFlag::_normal; + ; + dof_data.local_equation_number.push_back(first_local_dof_id); + + this->global_equation_number(first_local_dof_id) = first_global_dof_id; + this->global_to_local_mapping[first_global_dof_id] = first_local_dof_id; + + ++first_global_dof_id; + ++first_local_dof_id; } } /* -------------------------------------------------------------------------- */ // register in factory static bool default_dof_manager_is_registered[[gnu::unused]] = DefaultDOFManagerFactory::getInstance().registerAllocator( - "default", - [](const ID & id, - const MemoryID & mem_id) -> std::unique_ptr { + "default", [](const ID & id, + const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(id, mem_id); }); static bool dof_manager_is_registered[[gnu::unused]] = DOFManagerFactory::getInstance().registerAllocator( - "default", - [](Mesh & mesh, const ID & id, - const MemoryID & mem_id) -> std::unique_ptr { + "default", [](Mesh & mesh, const ID & id, + const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(mesh, id, mem_id); }); } // namespace akantu diff --git a/src/model/dof_manager_default.hh b/src/model/dof_manager_default.hh index f6bda2457..2ed32c812 100644 --- a/src/model/dof_manager_default.hh +++ b/src/model/dof_manager_default.hh @@ -1,360 +1,364 @@ /** * @file dof_manager_default.hh * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Jan 31 2018 * * @brief Default implementation of the dof manager * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_DEFAULT_HH__ #define __AKANTU_DOF_MANAGER_DEFAULT_HH__ namespace akantu { class SparseMatrixAIJ; class NonLinearSolverDefault; class TimeStepSolverDefault; class DOFSynchronizer; } // namespace akantu 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); ~DOFManagerDefault() override; protected: struct DOFDataDefault : public DOFData { explicit DOFDataDefault(const ID & dof_id); /// associated node for _dst_nodal dofs only Array associated_nodes; }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ private: void registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs); public: /// register an array of degree of freedom void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) override; /// the dof as an implied type of _dst_nodal and is defined only on a subset /// of nodes void registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support) override; /// Assemble an array to the global residual array - void assembleToResidual(const ID & dof_id, - Array & array_to_assemble, + void assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor = 1.) override; /// Assemble an array to the global lumped matrix array void assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor = 1.) override; /** * 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 **/ 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) override; /// multiply a vector by a matrix and assemble the result to the residual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) override; /// multiply a vector by a lumped matrix and assemble the result to the /// residual void assembleLumpedMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) override; /// assemble coupling terms between to dofs void assemblePreassembledMatrix(const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) override; 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); template - void makeConsistentForPeriodicity(const ID & dof_id, - Array & array); + void makeConsistentForPeriodicity(const ID & dof_id, Array & array); public: /// clear the residual void clearResidual() override; /// sets the matrix to 0 void clearMatrix(const ID & mtx) override; /// sets the lumped matrix to 0 void clearLumpedMatrix(const ID & mtx) override; /// update the global dofs vector virtual void updateGlobalBlockedDofs(); /// apply boundary conditions to jacobian matrix virtual void applyBoundary(const ID & matrix_id = "J"); // void getEquationsNumbers(const ID & dof_id, // Array & equation_numbers) override; 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 void getSolutionPerDOFs(const ID & dof_id, Array & solution_array) override; /// 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 void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) override; 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); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ protected: std::pair updateNodalDOFs(const ID & dof_id, const Array & nodes_list) override; private: void updateDOFsData(DOFDataDefault & dof_data, UInt nb_new_local_dofs, - UInt nb_new_pure_local, + UInt nb_new_pure_local, UInt nb_nodes, const std::function & getNode); + void updateDOFsData(DOFDataDefault & dof_data, UInt nb_new_local_dofs, + UInt nb_new_pure_local); + + void resizeGlobalArrays(); + auto computeFirstDOFIDs(UInt nb_new_local_dofs, UInt nb_new_pure_local); + public: /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get an instance of a new SparseMatrix SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) override; /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) override; /// Get the reference of an existing matrix SparseMatrixAIJ & getMatrix(const ID & matrix_id); /* ------------------------------------------------------------------------ */ /* Non Linear Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver NonLinearSolver & getNewNonLinearSolver( const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type) override; /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver TimeStepSolver & getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) override; /* ------------------------------------------------------------------------ */ /// Get the solution array AKANTU_GET_MACRO_NOT_CONST(GlobalSolution, global_solution, Array &); /// Set the global solution array void setGlobalSolution(const Array & solution); /// Get the global residual array across processors (SMP call) const Array & getGlobalResidual(); /// Get the residual array const Array & getResidual() const; /// Get the blocked dofs array AKANTU_GET_MACRO(GlobalBlockedDOFs, global_blocked_dofs, const Array &); /// Get the blocked dofs array AKANTU_GET_MACRO(PreviousGlobalBlockedDOFs, previous_global_blocked_dofs, const Array &); /// Get the location type of a given dof inline bool isLocalOrMasterDOF(UInt local_dof_num); /// Answer to the question is a dof a slave dof ? inline bool isSlaveDOF(UInt local_dof_num); /// get the equation numbers (in local numbering) corresponding to a dof ID inline const Array & getLocalEquationNumbers(const ID & dof_id) const; /// tells if the dof manager knows about a global dof bool hasGlobalEquationNumber(UInt global) const; /// return the local index of the global equation number inline UInt globalToLocalEquationNumber(UInt global) const; /// converts local equation numbers to global equation numbers; inline UInt localToGlobalEquationNumber(UInt local) const; /// get the array of dof types (use only if you know what you do...) inline NodeFlag getDOFFlag(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: DOFData & getNewDOFData(const ID & dof_id) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend class GlobalDOFInfoDataAccessor; // using AIJMatrixMap = std::map>; // using DefaultNonLinearSolversMap = // std::map>; // using DefaultTimeStepSolversMap = // std::map>; using DOFToMatrixProfile = std::map, std::vector>>; /// 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 std::unique_ptr> 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_flag; /// 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{0}; /// equation number in global numbering Array global_equation_number; using equation_numbers_map = std::unordered_map; /// dual information of global_equation_number equation_numbers_map global_to_local_mapping; /// accumulator to know what would be the next global id to use UInt first_global_dof_id{0}; /// synchronizer to maintain coherency in dof fields std::unique_ptr synchronizer; }; } // namespace akantu #include "dof_manager_default_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_DEFAULT_HH__ */