diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 9f942f591..449c1b3c3 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,518 +1,519 @@ /** * @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 SpacialDirection { _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 }; /// Define the node/dof type enum NodeType : Int { _nt_pure_gost = -3, _nt_master = -2, _nt_normal = -1 }; /// 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 }; } #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/common/aka_common_inline_impl.cc b/src/common/aka_common_inline_impl.cc index b73e5da3f..3a0023066 100644 --- a/src/common/aka_common_inline_impl.cc +++ b/src/common/aka_common_inline_impl.cc @@ -1,455 +1,458 @@ /** * @file aka_common_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief inline implementations of common akantu type descriptions * * @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 * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_COMMON_INLINE_IMPL_CC__ #define __AKANTU_AKA_COMMON_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type) { switch (type) { case _not_ghost: stream << "not_ghost"; break; case _ghost: stream << "ghost"; break; case _casper: stream << "Casper the friendly ghost"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for TimeStepSolverType inline std::ostream & operator<<(std::ostream & stream, const TimeStepSolverType & type) { switch (type) { case _tsst_static: stream << "static"; break; case _tsst_dynamic: stream << "dynamic"; break; case _tsst_dynamic_lumped: stream << "dynamic_lumped"; break; case _tsst_not_defined: stream << "not defined time step solver"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard input stream operator for TimeStepSolverType inline std::istream & operator>>(std::istream & stream, TimeStepSolverType & type) { std::string str; stream >> str; if (str == "static") type = _tsst_static; else if (str == "dynamic") type = _tsst_dynamic; else if (str == "dynamic_lumped") type = _tsst_dynamic_lumped; else { AKANTU_ERROR("The type " << str << " is not a recognized TimeStepSolverType"); stream.setstate(std::ios::failbit); } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for NonLinearSolverType inline std::ostream & operator<<(std::ostream & stream, const NonLinearSolverType & type) { switch (type) { case _nls_linear: stream << "linear"; break; case _nls_newton_raphson: stream << "newton_raphson"; break; case _nls_newton_raphson_modified: stream << "newton_raphson_modified"; break; case _nls_lumped: stream << "lumped"; break; case _nls_auto: stream << "auto"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard input stream operator for NonLinearSolverType inline std::istream & operator>>(std::istream & stream, NonLinearSolverType & type) { std::string str; stream >> str; if (str == "linear") type = _nls_linear; else if (str == "newton_raphson") type = _nls_newton_raphson; else if (str == "newton_raphson_modified") type = _nls_newton_raphson_modified; else if (str == "lumped") type = _nls_lumped; else if (str == "auto") type = _nls_auto; else type = _nls_auto; return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for IntegrationSchemeType inline std::ostream & operator<<(std::ostream & stream, const IntegrationSchemeType & type) { switch (type) { case _ist_pseudo_time: stream << "pseudo_time"; break; case _ist_forward_euler: stream << "forward_euler"; break; case _ist_trapezoidal_rule_1: stream << "trapezoidal_rule_1"; break; case _ist_backward_euler: stream << "backward_euler"; break; case _ist_central_difference: stream << "central_difference"; break; case _ist_fox_goodwin: stream << "fox_goodwin"; break; case _ist_trapezoidal_rule_2: stream << "trapezoidal_rule_2"; break; case _ist_linear_acceleration: stream << "linear_acceleration"; break; case _ist_newmark_beta: stream << "newmark_beta"; break; case _ist_generalized_trapezoidal: stream << "generalized_trapezoidal"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard input stream operator for IntegrationSchemeType inline std::istream & operator>>(std::istream & stream, IntegrationSchemeType & type) { std::string str; stream >> str; if (str == "pseudo_time") type = _ist_pseudo_time; else if (str == "forward_euler") type = _ist_forward_euler; else if (str == "trapezoidal_rule_1") type = _ist_trapezoidal_rule_1; else if (str == "backward_euler") type = _ist_backward_euler; else if (str == "central_difference") type = _ist_central_difference; else if (str == "fox_goodwin") type = _ist_fox_goodwin; else if (str == "trapezoidal_rule_2") type = _ist_trapezoidal_rule_2; else if (str == "linear_acceleration") type = _ist_linear_acceleration; else if (str == "newmark_beta") type = _ist_newmark_beta; else if (str == "generalized_trapezoidal") type = _ist_generalized_trapezoidal; else { AKANTU_ERROR("The type " << str << " is not a recognized IntegrationSchemeType"); stream.setstate(std::ios::failbit); } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type) { switch (type) { case _gst_whatever: stream << "_gst_whatever"; break; + case _gst_ask_nodes: + stream << "_gst_ask_nodes"; + break; case _gst_update: stream << "_gst_update"; break; case _gst_size: stream << "_gst_size"; break; case _gst_smm_mass: stream << "_gst_smm_mass"; break; case _gst_smm_for_gradu: stream << "_gst_smm_for_gradu"; break; case _gst_smm_boundary: stream << "_gst_smm_boundary"; break; case _gst_smm_uv: stream << "_gst_smm_uv"; break; case _gst_smm_res: stream << "_gst_smm_res"; break; case _gst_smm_init_mat: stream << "_gst_smm_init_mat"; break; case _gst_smm_stress: stream << "_gst_smm_stress"; break; case _gst_smmc_facets: stream << "_gst_smmc_facets"; break; case _gst_smmc_facets_conn: stream << "_gst_smmc_facets_conn"; break; case _gst_smmc_facets_stress: stream << "_gst_smmc_facets_stress"; break; case _gst_smmc_damage: stream << "_gst_smmc_damage"; break; case _gst_giu_global_conn: stream << "_gst_giu_global_conn"; break; case _gst_ce_groups: stream << "_gst_ce_groups"; break; case _gst_gm_clusters: stream << "_gst_gm_clusters"; break; case _gst_htm_temperature: stream << "_gst_htm_temperature"; break; case _gst_htm_gradient_temperature: stream << "_gst_htm_gradient_temperature"; break; case _gst_htm_phi: stream << "_gst_htm_phi"; break; case _gst_htm_gradient_phi: stream << "_gst_htm_gradient_phi"; break; case _gst_mnl_for_average: stream << "_gst_mnl_for_average"; break; case _gst_mnl_weight: stream << "_gst_mnl_weight"; break; case _gst_nh_criterion: stream << "_gst_nh_criterion"; break; case _gst_test: stream << "_gst_test"; break; case _gst_user_1: stream << "_gst_user_1"; break; case _gst_user_2: stream << "_gst_user_2"; break; case _gst_material_id: stream << "_gst_material_id"; break; case _gst_for_dump: stream << "_gst_for_dump"; break; case _gst_cf_nodal: stream << "_gst_cf_nodal"; break; case _gst_cf_incr: stream << "_gst_cf_incr"; break; case _gst_solver_solution: stream << "_gst_solver_solution"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SolveConvergenceCriteria inline std::ostream & operator<<(std::ostream & stream, const SolveConvergenceCriteria & criteria) { switch (criteria) { case _scc_residual: stream << "_scc_residual"; break; case _scc_solution: stream << "_scc_solution"; break; case _scc_residual_mass_wgh: stream << "_scc_residual_mass_wgh"; break; } return stream; } inline std::istream & operator>>(std::istream & stream, SolveConvergenceCriteria & criteria) { std::string str; stream >> str; if (str == "residual") criteria = _scc_residual; else if (str == "solution") criteria = _scc_solution; else if (str == "residual_mass_wgh") criteria = _scc_residual_mass_wgh; else { stream.setstate(std::ios::failbit); } return stream; } /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str) { std::string lstr = str; std::transform(lstr.begin(), lstr.end(), lstr.begin(), (int (*)(int))tolower); return lstr; } /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim) { std::string trimed = to_trim; // left trim trimed.erase(trimed.begin(), std::find_if(trimed.begin(), trimed.end(), std::not1(std::ptr_fun(isspace)))); // right trim trimed.erase(std::find_if(trimed.rbegin(), trimed.rend(), std::not1(std::ptr_fun(isspace))) .base(), trimed.end()); return trimed; } } // akantu #include namespace akantu { /* -------------------------------------------------------------------------- */ template std::string printMemorySize(UInt size) { Real real_size = size * sizeof(T); UInt mult = 0; if (real_size != 0) mult = (std::log(real_size) / std::log(2)) / 10; std::stringstream sstr; real_size /= Real(1 << (10 * mult)); sstr << std::setprecision(2) << std::fixed << real_size; std::string size_prefix; switch (mult) { case 0: sstr << ""; break; case 1: sstr << "Ki"; break; case 2: sstr << "Mi"; break; case 3: sstr << "Gi"; break; // I started on this type of machines // (32bit computers) (Nicolas) case 4: sstr << "Ti"; break; case 5: sstr << "Pi"; break; case 6: sstr << "Ei"; break; // theoritical limit of RAM of the current // computers in 2014 (64bit computers) (Nicolas) case 7: sstr << "Zi"; break; case 8: sstr << "Yi"; break; default: AKANTU_ERROR( "The programmer in 2014 didn't thought so far (even wikipedia does not " "go further)." << " You have at least 1024 times more than a yobibit of RAM!!!" << " Just add the prefix corresponding in this switch case."); } sstr << "Byte"; return sstr.str(); } } // akantu #endif /* __AKANTU_AKA_COMMON_INLINE_IMPL_CC__ */ diff --git a/src/mesh_utils/global_ids_updater.cc b/src/mesh_utils/global_ids_updater.cc index 1fddc36c1..897e02066 100644 --- a/src/mesh_utils/global_ids_updater.cc +++ b/src/mesh_utils/global_ids_updater.cc @@ -1,110 +1,112 @@ /** * @file global_ids_updater.cc * * @author Marco Vocialta * * @date creation: Fri Apr 13 2012 * @date last modification: Fri Dec 08 2017 * * @brief Functions of the GlobalIdsUpdater * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "global_ids_updater.hh" #include "element_synchronizer.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { UInt GlobalIdsUpdater::updateGlobalIDs(UInt local_nb_new_nodes) { UInt total_nb_new_nodes = this->updateGlobalIDsLocally(local_nb_new_nodes); this->synchronizeGlobalIDs(); return total_nb_new_nodes; } UInt GlobalIdsUpdater::updateGlobalIDsLocally(UInt local_nb_new_nodes) { const auto & comm = mesh.getCommunicator(); Int rank = comm.whoAmI(); Int nb_proc = comm.getNbProc(); if (nb_proc == 1) return local_nb_new_nodes; /// resize global ids array Array & nodes_global_ids = mesh.getGlobalNodesIds(); UInt old_nb_nodes = mesh.getNbNodes() - local_nb_new_nodes; nodes_global_ids.resize(mesh.getNbNodes(), -1); /// compute the number of global nodes based on the number of old nodes Matrix local_master_nodes(2, nb_proc, 0); for (UInt n = 0; n < old_nb_nodes; ++n) if (mesh.isLocalOrMasterNode(n)) ++local_master_nodes(0, rank); /// compute amount of local or master doubled nodes for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) if (mesh.isLocalOrMasterNode(n)) ++local_master_nodes(1, rank); comm.allGather(local_master_nodes); local_master_nodes = local_master_nodes.transpose(); UInt old_global_nodes = std::accumulate(local_master_nodes(0).storage(), local_master_nodes(0).storage() + nb_proc, 0); /// update global number of nodes UInt total_nb_new_nodes = std::accumulate(local_master_nodes(1).storage(), local_master_nodes(1).storage() + nb_proc, 0); if (total_nb_new_nodes == 0) return 0; /// set global ids of local and master nodes UInt starting_index = std::accumulate(local_master_nodes(1).storage(), local_master_nodes(1).storage() + rank, old_global_nodes); for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) { if (mesh.isLocalOrMasterNode(n)) { nodes_global_ids(n) = starting_index; ++starting_index; + } else { + nodes_global_ids(n) = -1; } } MeshAccessor mesh_accessor(mesh); mesh_accessor.setNbGlobalNodes(old_global_nodes + total_nb_new_nodes); return total_nb_new_nodes; } void GlobalIdsUpdater::synchronizeGlobalIDs() { this->synchronizer.slaveReductionOnce(*this, _gst_giu_global_conn); this->synchronizer.synchronizeOnce(*this, _gst_giu_global_conn); } } // akantu diff --git a/src/mesh_utils/global_ids_updater_inline_impl.cc b/src/mesh_utils/global_ids_updater_inline_impl.cc index e4ce2707b..461a07ba6 100644 --- a/src/mesh_utils/global_ids_updater_inline_impl.cc +++ b/src/mesh_utils/global_ids_updater_inline_impl.cc @@ -1,121 +1,127 @@ /** * @file global_ids_updater_inline_impl.cc * * @author Marco Vocialta * * @date creation: Fri Oct 02 2015 * @date last modification: Sun Aug 13 2017 * * @brief Implementation of the inline functions of GlobalIdsUpdater * * @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 "global_ids_updater.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ #define __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt GlobalIdsUpdater::getNbData(const Array & elements, const SynchronizationTag & tag) const { UInt size = 0; if (tag == _gst_giu_global_conn) size += Mesh::getNbNodesPerElementList(elements) * sizeof(UInt); return size; } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag == _gst_giu_global_conn) packUnpackGlobalConnectivity(buffer, elements); } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag == _gst_giu_global_conn) packUnpackGlobalConnectivity(buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void GlobalIdsUpdater::packUnpackGlobalConnectivity( CommunicationBuffer & buffer, const Array & elements) const { AKANTU_DEBUG_IN(); ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; Array::const_vector_iterator conn_begin; UInt nb_nodes_per_elem = 0; auto & global_nodes_ids = mesh.getGlobalNodesIds(); - for(auto & el : elements) { + for (auto & el : elements) { if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; const Array & connectivity = mesh.getConnectivity(current_element_type, current_ghost_type); nb_nodes_per_elem = connectivity.getNbComponent(); conn_begin = connectivity.begin(nb_nodes_per_elem); } /// get element connectivity const Vector current_conn = conn_begin[el.element]; /// loop on all connectivity nodes for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = current_conn(n); if (pack_mode) { - if(mesh.isPureGhostNode(node)) { + if (mesh.isPureGhostNode(node)) { buffer << UInt(-1); } else { UInt index = global_nodes_ids(node); buffer << index; } } else { UInt index; buffer >> index; - if (global_nodes_ids(node) == UInt(-1)) + if (global_nodes_ids(node) == UInt(-1) and mesh.isSlaveNode(node)) { global_nodes_ids(node) = index; + } else { + AKANTU_DEBUG_ASSERT( + index == UInt(-1) or mesh.isPureGhostNode(node) or index == global_nodes_ids(node), + "Two processors disagree on the global number of a node " + << index << " against " << global_nodes_ids(node)); + } } } } AKANTU_DEBUG_OUT(); } } // akantu #endif /* __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ */ diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 29484a9fb..39bee7f7e 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,917 +1,920 @@ /** * @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 "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_type(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_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 = std::make_unique( *this, this->id + ":dof_synchronizer", this->memory_id); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() = default; /* -------------------------------------------------------------------------- */ 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.size() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.size() == 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) { this->dofs[dof_id] = std::make_unique(dof_id); return *this->dofs[dof_id]; } -/* -------------------------------------------------------------------------- */ -class GlobalDOFInfoDataAccessor : public DataAccessor { -public: - using size_type = - typename std::unordered_map>::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 { - auto it = dofs_per_node.find(node); - AKANTU_DEBUG_ASSERT(it != dofs_per_node.end(), - "Did not find the node " << node); - return it->second[nth_dof]; - } - - UInt getNbData(const Array & nodes, - const SynchronizationTag & tag) const override { - if (tag == _gst_size) { - return nodes.size() * 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; - } - - void packData(CommunicationBuffer & buffer, const Array & nodes, - const SynchronizationTag & tag) const override { - for (auto node : nodes) { - auto it = dofs_per_node.find(node); - if (tag == _gst_size) { - if (it != dofs_per_node.end()) { - buffer << size_type(it->second.size()); - } else { - buffer << size_type(0); - } - } else if (tag == _gst_update) { - if (it != dofs_per_node.end()) - buffer << it->second; - } - } - } - - void unpackData(CommunicationBuffer & buffer, const Array & nodes, - const SynchronizationTag & tag) override { - for (auto node : nodes) { - auto it = dofs_per_node.find(node); - if (tag == _gst_size) { - size_type size; - buffer >> size; - if (size != 0) - dofs_per_node[node].resize(size); - } else if (tag == _gst_update) { - if (it != dofs_per_node.end()) - buffer >> it->second; - } - } - } - -protected: - std::unordered_map> dofs_per_node; -}; - /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs) { // 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; }); } // 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, 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.); 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(); DOFData & 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 Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & 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 Array & 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.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.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 Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt size = A.size(); 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.size(); 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.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, [&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, + const DOFManagerDefault & dof_manager, + const Mesh & mesh) + : dof_data(dof_data) { + for (auto && pair : + zip(dof_data.local_equation_number, dof_data.associated_nodes)) { + UInt node, dof; + std::tie(dof, node) = pair; + if (mesh.isMasterNode(node)) { + dofs_per_node[node].push_back( + dof_manager.localToGlobalEquationNumber(dof)); + } + } + } + + UInt getNthDOFForNode(UInt nth_dof, UInt node) const { + auto it = dofs_per_node.find(node); + AKANTU_DEBUG_ASSERT(it != dofs_per_node.end(), + "Did not find the node " << node); + return it->second[nth_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); + } + + 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; + } + + void packData(CommunicationBuffer & buffer, const Array & nodes, + const SynchronizationTag & tag) const override { + if (tag == _gst_ask_nodes) { + for (auto & node : nodes) { + for (auto & dof : dofs_per_node.at(node)) { + buffer << dof; + } + } + } + } + + void unpackData(CommunicationBuffer & buffer, const Array & nodes, + const SynchronizationTag & tag) override { + if (tag == _gst_ask_nodes) { + auto nb_dofs_per_node = dof_data.dof->getNbComponent(); + for (auto & node : nodes) { + auto & dofs = dofs_per_node[node]; + for (auto _[[gnu::unused]] : arange(nb_dofs_per_node)) { + UInt dof; + buffer >> dof; + dofs.push_back(dof); + } + } + } + } + +protected: + std::unordered_map> dofs_per_node; + DOFManagerDefault::DOFDataDefault & dof_data; +}; + /* -------------------------------------------------------------------------- */ 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; - GlobalDOFInfoDataAccessor data_accessor; - // resize all relevant arrays this->residual.resize(this->local_system_size, 0.); this->dofs_type.resize(this->local_system_size); 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); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); // 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); } // update per dof info for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number.push_back(local_eq_num); bool is_local_dof = true; // determine the dof type switch (support_type) { case _dst_nodal: { UInt node = getNode(d / dof_data.dof->getNbComponent()); this->dofs_type(local_eq_num) = this->mesh->getNodeType(node); dof_data.associated_nodes.push_back(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; + this->global_equation_number(local_eq_num) = UInt(-1); } } // synchronize the global numbering for slaves if (support_type == _dst_nodal && this->synchronizer) { auto nb_dofs_per_node = dof_data.dof->getNbComponent(); + GlobalDOFInfoDataAccessor data_accessor(dof_data, *this, *this->mesh); auto & node_synchronizer = this->mesh->getNodeSynchronizer(); - node_synchronizer.synchronizeOnce(data_accessor, _gst_size); - node_synchronizer.synchronizeOnce(data_accessor, _gst_update); + node_synchronizer.synchronizeOnce(data_accessor, _gst_ask_nodes); std::vector counters(nb_new_local_dofs / nb_dofs_per_node); for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; if (not this->isSlaveDOF(local_eq_num)) continue; UInt node = d / nb_dofs_per_node; UInt dof_count = counters[node]++; node = getNode(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; } } + + + this->first_global_dof_id += std::accumulate( + nb_dofs_per_proc.begin() + prank + 1, nb_dofs_per_proc.end(), 0); + } /* -------------------------------------------------------------------------- */ // 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 { 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 { 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 e8efa3b50..6f96f03c8 100644 --- a/src/model/dof_manager_default.hh +++ b/src/model/dof_manager_default.hh @@ -1,361 +1,362 @@ /** * @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, const Array & array_to_assemble, Real scale_factor = 1.) override; /// Assemble an array to the global lumped matrix array void assembleToLumpedMatrix(const ID & dof_id, const 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); 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, const std::function & getNode); 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 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: 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_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{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__ */