diff --git a/.clang-tidy b/.clang-tidy index ce13071c2..d973c6bbf 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -1,4 +1,4 @@ -Checks: '-modernize-use-trailing-return-type, modernize-use-*, performance-*, mpi-*, openmp-*, readability-*, cppcoreguidelines-*' +Checks: 'modernize-use-*, -modernize-use-trailing-return-type*, performance-*, mpi-*, openmp-*, readability-*, cppcoreguidelines-*, -cppcoreguidelines-pro-bounds-array-to-pointer-decay*, -cppcoreguidelines-macro-usage*' AnalyzeTemporaryDtors: false HeaderFilterRegex: 'src/.*' FormatStyle: file diff --git a/src/common/aka_array.cc b/src/common/aka_array.cc index f02204b9c..52739c70f 100644 --- a/src/common/aka_array.cc +++ b/src/common/aka_array.cc @@ -1,100 +1,97 @@ /** * @file aka_array.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of akantu::Array * * * 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 #include /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /* Functions ArrayBase */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template <> UInt Array::find(const Real & elem) const { AKANTU_DEBUG_IN(); Real epsilon = std::numeric_limits::epsilon(); auto it = std::find_if(begin(), end(), [&elem, &epsilon](auto && a) { return std::abs(a - elem) <= epsilon; }); AKANTU_DEBUG_OUT(); return (it != end()) ? end() - it : UInt(-1); } /* -------------------------------------------------------------------------- */ template <> -Array & Array::operator*=(__attribute__((unused)) - const ElementType & alpha) { +Array & +Array::operator*=(const ElementType & /*alpha*/) { AKANTU_TO_IMPLEMENT(); return *this; } template <> -Array & Array:: -operator-=(__attribute__((unused)) const Array & vect) { +Array & +Array::operator-=(__attribute__((unused)) + const Array & /*vect*/) { AKANTU_TO_IMPLEMENT(); return *this; } template <> -Array & Array:: -operator+=(__attribute__((unused)) const Array & vect) { +Array & +Array::operator+=(const Array & /*vect*/) { AKANTU_TO_IMPLEMENT(); return *this; } -template <> -Array & Array::operator*=(__attribute__((unused)) - const char & alpha) { +template <> Array & Array::operator*=(const char & /*alpha*/) { AKANTU_TO_IMPLEMENT(); return *this; } template <> -Array & Array::operator-=(__attribute__((unused)) - const Array & vect) { +Array & Array::operator-=(const Array & /*vect*/) { AKANTU_TO_IMPLEMENT(); return *this; } template <> -Array & Array::operator+=(__attribute__((unused)) - const Array & vect) { +Array & Array::operator+=(const Array & /*vect*/) { AKANTU_TO_IMPLEMENT(); return *this; } } // namespace akantu diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index be17e197d..be2fe7c77 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,642 +1,643 @@ /** * @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 * * * 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #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 #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Constants */ /* -------------------------------------------------------------------------- */ namespace { [[gnu::unused]] constexpr UInt _all_dimensions{ std::numeric_limits::max()}; #ifdef AKANTU_NDEBUG [[gnu::unused]] constexpr Real REAL_INIT_VALUE{0.}; #else [[gnu::unused]] constexpr Real REAL_INIT_VALUE{ std::numeric_limits::quiet_NaN()}; #endif } // namespace /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; using MemoryID = UInt; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "aka_enum_macros.hh" /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* 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 }; #if !defined(DOXYGEN) // 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_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_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 }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_NON_LINEAR_SOLVER_TYPES \ (linear) \ (newton_raphson) \ (newton_raphson_modified) \ (lumped) \ (gmres) \ (bfgs) \ (cg) \ (auto) // clang-format on AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) #else /// Type of non linear resolution available in akantu enum class NonLinearSolverType { _linear, ///< No non linear convergence loop _newton_raphson, ///< Regular Newton-Raphson _newton_raphson_modified, ///< Newton-Raphson with initial tangent _lumped, ///< Case of lumped mass or equivalent matrix _gmres, _bfgs, _cg, _auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_TIME_STEP_SOLVER_TYPE \ (static) \ (dynamic) \ (dynamic_lumped) \ (not_defined) // clang-format on AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) #else /// Type of time stepping solver enum class TimeStepSolverType { _static, ///< Static solution _dynamic, ///< Dynamic solver _dynamic_lumped, ///< Dynamic solver with lumped mass _not_defined, ///< For not defined cases }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_INTEGRATION_SCHEME_TYPE \ (pseudo_time) \ (forward_euler) \ (trapezoidal_rule_1) \ (backward_euler) \ (central_difference) \ (fox_goodwin) \ (trapezoidal_rule_2) \ (linear_acceleration) \ (newmark_beta) \ (generalized_trapezoidal) // clang-format on AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) #else /// Type of integration scheme enum class IntegrationSchemeType { _pseudo_time, ///< Pseudo Time _forward_euler, ///< GeneralizedTrapezoidal(0) _trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _backward_euler, ///< GeneralizedTrapezoidal(1) _central_difference, ///< NewmarkBeta(0, 1/2) _fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ (residual) \ (solution) \ (residual_mass_wgh) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) #else /// enum SolveConvergenceCriteria different convergence criteria enum class SolveConvergenceCriteria { _residual, ///< Use residual to test the convergence _solution, ///< Use solution to test the convergence _residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; #endif /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum MatrixType 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 }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_SYNCHRONIZATION_TAG \ (whatever) \ (update) \ (ask_nodes) \ (size) \ (smm_mass) \ (smm_for_gradu) \ (smm_boundary) \ (smm_uv) \ (smm_res) \ (smm_init_mat) \ (smm_stress) \ (smmc_facets) \ (smmc_facets_conn) \ (smmc_facets_stress) \ (smmc_damage) \ (giu_global_conn) \ (ce_groups) \ (gm_clusters) \ (htm_temperature) \ (htm_gradient_temperature) \ (htm_phi) \ (htm_gradient_phi) \ (mnl_for_average) \ (mnl_weight) \ (nh_criterion) \ (test) \ (user_1) \ (user_2) \ (material_id) \ (for_dump) \ (cf_nodal) \ (cf_incr) \ (solver_solution) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) #else /// @enum SynchronizationTag type of synchronizations enum class SynchronizationTag { //--- Generic tags --- _whatever, _update, _ask_nodes, _size, //--- SolidMechanicsModel tags --- _smm_mass, ///< synchronization of the SolidMechanicsModel.mass _smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _smm_uv, ///< synchronization of the nodal velocities and displacement _smm_res, ///< synchronization of the nodal residual _smm_init_mat, ///< synchronization of the data to initialize materials _smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _smmc_facets, ///< synchronization of facet data to setup facet synch _smmc_facets_conn, ///< synchronization of facet global connectivity _smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _htm_temperature, ///< synchronization of the nodal temperature _htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _htm_phi, ///< synchronization of the nodal level set value phi _htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _mnl_for_average, ///< synchronization of data to average in non local /// material _mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _nh_criterion, // --- General tags --- _test, ///< Test tag _user_1, ///< tag for user simulations _user_2, ///< tag for user simulations _material_id, ///< synchronization of the material ids _for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _cf_nodal, ///< synchronization of disp, velo, and current position _cf_incr, ///< synchronization of increment // --- Solver tags --- _solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; #endif /// @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 = 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)); } std::ostream & operator<<(std::ostream & stream, NodeFlag flag); } // namespace akantu AKANTU_ENUM_HASH(GhostType) 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; +namespace { + constexpr ghost_type_t ghost_types{_casper}; +} /// standard output stream operator for GhostType -inline std::ostream & operator<<(std::ostream & stream, GhostType type); +// inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT ' ' #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ #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_DEREF_PTR(name, ptr) \ inline decltype(auto) get##name() const { \ if (not ptr) { \ AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ } \ return (*ptr); \ } #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); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ struct TensorTrait {}; struct TensorProxyTrait {}; } // namespace akantu /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ namespace aka { /* ------------------------------------------------------------------------ */ template using is_tensor = std::is_base_of; template using is_tensor_proxy = std::is_base_of; /* ------------------------------------------------------------------------ */ template using is_scalar = std::is_arithmetic; /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> bool is_of_type(T && t) { return ( dynamic_cast>::value, std::add_const_t, R>>>(&t) != nullptr); } /* -------------------------------------------------------------------------- */ template bool is_of_type(std::unique_ptr & t) { return ( dynamic_cast::value, std::add_const_t, R>>>( t.get()) != nullptr); } /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { static_assert( disjunction< std::is_base_of, std::decay_t>, // down-cast std::is_base_of, std::decay_t> // up-cast >::value, "Type T and R are not valid for a as_type conversion"); return dynamic_cast>::value, std::add_const_t, R>>>(t); } /* -------------------------------------------------------------------------- */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { return &as_type(*t); } /* -------------------------------------------------------------------------- */ template decltype(auto) as_type(const std::shared_ptr & t) { return std::dynamic_pointer_cast(t); } } // namespace aka #include "aka_common_inline_impl.hh" #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(); #define AKANTU_CURRENT_FUNCTION \ (std::string(__func__) + "():" + std::to_string(__LINE__)) } // namespace akantu /* -------------------------------------------------------------------------- */ #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.hh b/src/common/aka_common_inline_impl.hh index 78e807f84..5403ad538 100644 --- a/src/common/aka_common_inline_impl.hh +++ b/src/common/aka_common_inline_impl.hh @@ -1,149 +1,149 @@ /** * @file aka_common_inline_impl.hh * * @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 * * * 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 "aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /// standard output stream operator for GhostType -inline std::ostream & operator<<(std::ostream & stream, GhostType type) { +inline std::ostream & operator<<(std::ostream & stream, const 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; } /* -------------------------------------------------------------------------- */ 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; } namespace { template inline std::string trim_p(const std::string & to_trim, pred && p) { std::string trimed = to_trim; auto && not_ = [&](auto && a) { return not p(a); }; // left trim trimed.erase(trimed.begin(), std::find_if(trimed.begin(), trimed.end(), not_)); // right trim trimed.erase(std::find_if(trimed.rbegin(), trimed.rend(), not_).base(), trimed.end()); return trimed; } } // namespace /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim) { return trim_p(to_trim, [&](auto && a) { return std::isspace(a); }); } inline std::string trim(const std::string & to_trim, char c) { return trim_p(to_trim, [&c](auto && a) { return (a == c); }); } /* -------------------------------------------------------------------------- */ 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(); } } // namespace akantu diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index 26f59ea60..c75d92d1d 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,356 +1,363 @@ /** * @file aka_error.cc * * @author Nicolas Richart * * @date creation: Mon Sep 06 2010 * @date last modification: Sun Dec 03 2017 * * @brief handling of errors * * * 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 "aka_error.hh" #include "aka_common.hh" #include "aka_config.hh" #include "aka_iterators.hh" +#include "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ #include #include #include #if (defined(READLINK_COMMAND) || defined(ADDR2LINE_COMMAND)) && \ (!defined(_WIN32)) #include #include #endif #include #include #include #include #include #include #include #include #include #ifdef AKANTU_USE_MPI #include #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace debug { // static void printBacktraceAndExit(int) { std::terminate(); } // /* ------------------------------------------------------------------------ */ // void initSignalHandler() { std::signal(SIGSEGV, &printBacktraceAndExit); } /* ------------------------------------------------------------------------ */ std::string demangle(const char * symbol) { int status; std::string result; char * demangled_name; if ((demangled_name = abi::__cxa_demangle(symbol, nullptr, nullptr, &status)) != nullptr) { result = demangled_name; free(demangled_name); } else { result = symbol; } return result; } /* ------------------------------------------------------------------------ */ #if (defined(READLINK_COMMAND) || defined(ADDR2LINK_COMMAND)) && \ (!defined(_WIN32)) std::string exec(const std::string & cmd) { FILE * pipe = popen(cmd.c_str(), "r"); if (!pipe) return ""; char buffer[1024]; std::string result = ""; while (!feof(pipe)) { if (fgets(buffer, 128, pipe) != nullptr) result += buffer; } result = result.substr(0, result.size() - 1); pclose(pipe); return result; } #endif auto getBacktrace() -> std::vector { std::vector backtrace_lines; #if not defined(_WIN32) #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) std::string me = ""; char buf[1024]; /* The manpage says it won't null terminate. Let's zero the buffer. */ memset(buf, 0, sizeof(buf)); /* Note we use sizeof(buf)-1 since we may need an extra char for NUL. */ if (readlink("/proc/self/exe", buf, sizeof(buf) - 1)) me = std::string(buf); std::ifstream inmaps; inmaps.open("/proc/self/maps"); std::map addr_map; std::string line; while (inmaps.good()) { std::getline(inmaps, line); std::stringstream sstr(line); size_t first = line.find('-'); std::stringstream sstra(line.substr(0, first)); size_t addr; sstra >> std::hex >> addr; std::string lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; if (lib != "" && addr_map.find(lib) == addr_map.end()) { addr_map[lib] = addr; } } if (me != "") addr_map[me] = 0; #endif /// \todo for windows this part could be coded using CaptureStackBackTrace /// and SymFromAddr const size_t max_depth = 100; size_t stack_depth; void * stack_addrs[max_depth]; char ** stack_strings; size_t i; stack_depth = backtrace(stack_addrs, max_depth); stack_strings = backtrace_symbols(stack_addrs, stack_depth); /// -1 to remove the call to the printBacktrace function for (i = 1; i < stack_depth; i++) { std::string bt_line(stack_strings[i]); size_t first, second; if ((first = bt_line.find('(')) != std::string::npos && (second = bt_line.find('+')) != std::string::npos) { std::string location = bt_line.substr(0, first); #if defined(READLINK_COMMAND) std::string location_cmd = std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + std::string(" -f ") + location; location = exec(location_cmd); #endif std::string call = demangle(bt_line.substr(first + 1, second - first - 1).c_str()); size_t f = bt_line.find('['); size_t s = bt_line.find(']'); std::string address = bt_line.substr(f + 1, s - f - 1); std::stringstream sstra(address); size_t addr; sstra >> std::hex >> addr; std::string trace = location + " [" + call + "]"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) auto it = addr_map.find(location); if (it != addr_map.end()) { std::stringstream syscom; syscom << BOOST_PP_STRINGIZE(ADDR2LINE_COMMAND) << " 0x" << std::hex << (addr - it->second) << " -i -e " << location; std::string line = exec(syscom.str()); trace += " (" + line + ")"; } else { #endif std::stringstream sstr_addr; sstr_addr << std::hex << addr; trace += " (0x" + sstr_addr.str() + ")"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) } #endif backtrace_lines.push_back(trace); } else { backtrace_lines.push_back(bt_line); } } free(stack_strings); #endif return backtrace_lines; } /* ------------------------------------------------------------------------ */ void printBacktrace(const std::vector & backtrace) { auto w = size_t(std::floor(std::log10(double(backtrace.size()))) + 1); std::cerr << "BACKTRACE : " << backtrace.size() << " stack frames.\n"; for (auto && data : enumerate(backtrace)) std::cerr << " [" << std::setw(w) << (std::get<0>(data) + 1) << "] " << std::get<1>(data) << "\n"; std::cerr << "END BACKTRACE" << std::endl; } /* ------------------------------------------------------------------------ */ namespace { void terminate_handler() { auto eptr = std::current_exception(); auto t = abi::__cxa_current_exception_type(); auto name = t ? demangle(t->name()) : std::string("unknown"); try { if (eptr) std::rethrow_exception(eptr); else { printBacktrace(); std::cerr << AKANTU_LOCATION << "!! Execution terminated for unknown reasons !!" << std::endl; } } catch (Exception & e) { printBacktrace(e.backtrace()); std::cerr << "!! Uncaught akantu::Exception of type " << name << " !!\nwhat(): \"" << e.what() << "\"" << std::endl; } catch (std::exception & e) { std::cerr << "!! Uncaught exception of type " << name << " !!\nwhat(): \"" << e.what() << "\"" << std::endl; } catch (...) { std::cerr << "!! Something strange of type \"" << name << "\" was thrown.... !!" << std::endl; } + + if (debugger.printBacktrace()) { + std::cerr << "Random generator seed: " << RandomGenerator::seed() + << std::endl; + printBacktrace(); + } } } // namespace /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() { cout = &std::cerr; level = dblWarning; parallel_context = ""; file_open = false; print_backtrace = false; //initSignalHandler(); std::set_terminate(terminate_handler); } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if (file_open) { dynamic_cast(cout)->close(); delete cout; } } /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { if (status != 0) std::terminate(); std::exit(0); } /*------------------------------------------------------------------------- */ void Debugger::throwException(const std::string & info, const std::string & file, unsigned int line, __attribute__((unused)) bool silent, __attribute__((unused)) const std::string & location, const std::string & module) const noexcept(false) { #if !defined(AKANTU_NDEBUG) if (not silent) { printMessage("###", dblWarning, info + " " + location, module); } #endif debug::Exception ex(info, file, line); ex.setModule(module); throw ex; } /* ------------------------------------------------------------------------ */ void Debugger::printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info, const std::string & module) const { if (testLevel(level, module)) { double timestamp = std::chrono::duration_cast>( std::chrono::system_clock::now().time_since_epoch()) .count(); *(cout) << parallel_context << "{" << (size_t)timestamp << "} " << prefix << " " << info << std::endl; } } /* ------------------------------------------------------------------------ */ void Debugger::setDebugLevel(const DebugLevel & level) { this->level = level; } /* ------------------------------------------------------------------------ */ const DebugLevel & Debugger::getDebugLevel() const { return this->level; } /* ------------------------------------------------------------------------ */ void Debugger::setLogFile(const std::string & filename) { if (file_open) { dynamic_cast(cout)->close(); delete cout; } auto * fileout = new std::ofstream(filename.c_str()); file_open = true; cout = fileout; } std::ostream & Debugger::getOutputStream() { return *cout; } /* ------------------------------------------------------------------------ */ void Debugger::setParallelContext(int rank, int size) { std::stringstream sstr; UInt pad = std::ceil(std::log10(size)); sstr << "<" << getpid() << ">[R" << std::setfill(' ') << std::right << std::setw(pad) << rank << "|S" << size << "] "; parallel_context = sstr.str(); } void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } const DebugLevel & getDebugLevel() { return debugger.getDebugLevel(); } /* -------------------------------------------------------------------------- */ void exit(int status) { debugger.exit(status); } } // namespace debug } // namespace akantu diff --git a/src/common/aka_error.hh b/src/common/aka_error.hh index 0ab2e17af..3badb946f 100644 --- a/src/common/aka_error.hh +++ b/src/common/aka_error.hh @@ -1,405 +1,415 @@ /** * @file aka_error.hh * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Tue Feb 20 2018 * * @brief error management and internal exceptions * * * 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 #include #include #include #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ERROR_HH__ #define __AKANTU_ERROR_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ enum DebugLevel { dbl0 = 0, dblError = 0, dblAssert = 0, dbl1 = 1, dblException = 1, dblCritical = 1, dbl2 = 2, dblMajor = 2, dbl3 = 3, dblCall = 3, dblSecondary = 3, dblHead = 3, dbl4 = 4, dblWarning = 4, dbl5 = 5, dblInfo = 5, dbl6 = 6, dblIn = 6, dblOut = 6, dbl7 = 7, dbl8 = 8, dblTrace = 8, dbl9 = 9, dblAccessory = 9, dbl10 = 10, dblDebug = 42, dbl100 = 100, dblDump = 100, dblTest = 1337 }; /* -------------------------------------------------------------------------- */ #define AKANTU_LOCATION \ - "(" << __func__ << "(): " << __FILE__ << ":" << __LINE__ << ")" + "(" << std::string(__func__) << "(): " << std::string(__FILE__) << ":" \ + << std::to_string(__LINE__) \ + << ")" // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay) /* -------------------------------------------------------------------------- */ namespace debug { void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel(); void initSignalHandler(); std::string demangle(const char * symbol); template std::string demangle() { return demangle(typeid(T).name()); } template std::string demangle(const T & t) { return demangle(typeid(t).name()); } auto exec(const std::string & cmd) -> std::string; auto getBacktrace() -> std::vector; void printBacktrace(const std::vector & backtrace = getBacktrace()); void exit(int status) __attribute__((noreturn)); /* ------------------------------------------------------------------------ */ /// exception class that can be thrown by akantu class Exception : public std::exception { /* ---------------------------------------------------------------------- */ /* Constructors/Destructors */ /* ---------------------------------------------------------------------- */ protected: explicit Exception(std::string info = "") : _info(std::move(info)), _file("") {} public: //! full constructor Exception(std::string info, std::string file, unsigned int line) : _info(std::move(info)), _file(std::move(file)), _line(line) {} //! destructor ~Exception() noexcept override = default; /* ---------------------------------------------------------------------- */ /* Methods */ /* ---------------------------------------------------------------------- */ public: const char * what() const noexcept override { return _info.c_str(); } virtual const std::string info() const noexcept { std::stringstream stream; stream << debug::demangle(typeid(*this).name()) << " : " << _info << " [" << _file << ":" << _line << "]"; return stream.str(); } public: void setInfo(const std::string & info) { _info = info; } void setFile(const std::string & file) { _file = file; } void setLine(unsigned int line) { _line = line; } void setModule(const std::string & module) { _module = module; } void setBacktrace(const std::vector & backtrace) { backtrace_ = backtrace; } decltype(auto) backtrace() const { return backtrace_; } /* ---------------------------------------------------------------------- */ /* Class Members */ /* ---------------------------------------------------------------------- */ protected: /// exception description and additionals std::string _info; private: /// file it is thrown from std::string _file; /// line it is thrown from unsigned int _line{0}; /// module in which exception was raised std::string _module{"core"}; std::vector backtrace_; }; class CriticalError : public Exception {}; class AssertException : public Exception {}; class NotImplementedException : public Exception {}; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Exception & _this) { stream << _this.what(); return stream; } /* -------------------------------------------------------------------------- */ class Debugger { public: Debugger(); virtual ~Debugger(); Debugger(const Debugger &) = default; Debugger & operator=(const Debugger &) = default; void exit(int status) __attribute__((noreturn)); void throwException(const std::string & info, const std::string & file, unsigned int line, bool, const std::string &, const std::string & module) const noexcept(false) __attribute__((noreturn)); /*----------------------------------------------------------------------- */ template void throwCustomException(const Except & ex, const std::string & info, const std::string & file, unsigned int line, const std::string & module) const noexcept(false) __attribute__((noreturn)); /*----------------------------------------------------------------------- */ template void throwCustomException(const Except & ex, const std::string & file, unsigned int line, const std::string & module_) const noexcept(false) __attribute__((noreturn)); void printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info, const std::string & module_) const; void setOutStream(std::ostream & out) { cout = &out; } std::ostream & getOutStream() { return *cout; } public: void setParallelContext(int rank, int size); void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel() const; void setLogFile(const std::string & filename); std::ostream & getOutputStream(); inline bool testLevel(const DebugLevel & level, const std::string & module = "core") const { auto level_reached = (this->level >= (level)); auto correct_module = (level <= dblCritical) or (modules_to_debug.size() == 0) or (modules_to_debug.find(module) != modules_to_debug.end()); return level_reached and correct_module; } void printBacktrace(bool on_off) { this->print_backtrace = on_off; } bool printBacktrace() { return this->print_backtrace; } void addModuleToDebug(const std::string & id) { modules_to_debug.insert(id); } void removeModuleToDebug(const std::string & id) { auto it = modules_to_debug.find(id); if (it != modules_to_debug.end()) modules_to_debug.erase(it); } void listModules() { for (auto & module_ : modules_to_debug) { (*cout) << module_ << std::endl; } } private: std::string parallel_context; std::ostream * cout; bool file_open; DebugLevel level; bool print_backtrace; std::set modules_to_debug; }; extern Debugger debugger; } // namespace debug /* -------------------------------------------------------------------------- */ #define AKANTU_STRINGIZE_(str) #str #define AKANTU_STRINGIZE(str) AKANTU_STRINGIZE_(str) /* -------------------------------------------------------------------------- */ #define AKANTU_DEBUG_MODULE AKANTU_STRINGIZE(AKANTU_MODULE) /* -------------------------------------------------------------------------- */ #define AKANTU_STRINGSTREAM_IN(_str, _sstr) \ ; \ do { \ std::stringstream _dbg_s_info; \ _dbg_s_info << _sstr; \ _str = _dbg_s_info.str(); \ } while (false) /* -------------------------------------------------------------------------- */ #define AKANTU_EXCEPTION(info) AKANTU_EXCEPTION_(info, false) #define AKANTU_SILENT_EXCEPTION(info) AKANTU_EXCEPTION_(info, true) #define AKANTU_EXCEPTION_(info, silent) \ do { \ std::stringstream _dbg_str; \ _dbg_str << info; \ std::stringstream _dbg_loc; \ _dbg_loc << AKANTU_LOCATION; \ ::akantu::debug::debugger.throwException(_dbg_str.str(), __FILE__, \ __LINE__, silent, _dbg_loc.str(), \ AKANTU_DEBUG_MODULE); \ } while (false) #define AKANTU_CUSTOM_EXCEPTION_INFO(ex, info) \ do { \ std::stringstream _dbg_str; \ _dbg_str << info; \ ::akantu::debug::debugger.throwCustomException( \ ex, _dbg_str.str(), __FILE__, __LINE__, AKANTU_DEBUG_MODULE); \ } while (false) #define AKANTU_CUSTOM_EXCEPTION(ex) \ do { \ ::akantu::debug::debugger.throwCustomException(ex, __FILE__, __LINE__, \ AKANTU_DEBUG_MODULE); \ } while (false) /* -------------------------------------------------------------------------- */ #ifdef AKANTU_NDEBUG #define AKANTU_DEBUG_TEST(level) (false) #define AKANTU_DEBUG_LEVEL_IS_TEST() \ (::akantu::debug::debugger.testLevel(dblTest, AKANTU_DEBUG_MODULE)) #define AKANTU_DEBUG(level, info) #define AKANTU_DEBUG_(pref, level, info) #define AKANTU_DEBUG_IN() #define AKANTU_DEBUG_OUT() #define AKANTU_DEBUG_INFO(info) #define AKANTU_DEBUG_WARNING(info) #define AKANTU_DEBUG_TRACE(info) #define AKANTU_DEBUG_ASSERT(test, info) #define AKANTU_ERROR(info) \ AKANTU_CUSTOM_EXCEPTION_INFO(::akantu::debug::CriticalError(), info) /* -------------------------------------------------------------------------- */ #else #define AKANTU_DEBUG(level, info) AKANTU_DEBUG_(" ", level, info) #define AKANTU_DEBUG_(pref, level, info) \ do { \ std::string _dbg_str; \ AKANTU_STRINGSTREAM_IN(_dbg_str, info << " " << AKANTU_LOCATION); \ ::akantu::debug::debugger.printMessage(pref, level, _dbg_str, \ AKANTU_DEBUG_MODULE); \ } while (false) #define AKANTU_DEBUG_TEST(level) \ (::akantu::debug::debugger.testLevel(level, AKANTU_DEBUG_MODULE)) #define AKANTU_DEBUG_LEVEL_IS_TEST() \ (::akantu::debug::debugger.testLevel(dblTest)) #define AKANTU_DEBUG_IN() \ - AKANTU_DEBUG_("==>", ::akantu::dblIn, __func__ << "()") + AKANTU_DEBUG_( \ + "==>", ::akantu::dblIn, \ + __func__ \ + << "()") // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay) #define AKANTU_DEBUG_OUT() \ - AKANTU_DEBUG_("<==", ::akantu::dblOut, __func__ << "()") + AKANTU_DEBUG_( \ + "<==", ::akantu::dblOut, \ + __func__ \ + << "()") // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay) #define AKANTU_DEBUG_INFO(info) AKANTU_DEBUG_("---", ::akantu::dblInfo, info) #define AKANTU_DEBUG_WARNING(info) \ AKANTU_DEBUG_("/!\\", ::akantu::dblWarning, info) #define AKANTU_DEBUG_TRACE(info) AKANTU_DEBUG_(">>>", ::akantu::dblTrace, info) #define AKANTU_DEBUG_ASSERT(test, info) \ do { \ if (not(test)) \ AKANTU_CUSTOM_EXCEPTION_INFO(::akantu::debug::AssertException(), \ "assert [" << #test << "] " << info); \ } while (false) #define AKANTU_ERROR(info) \ do { \ AKANTU_DEBUG_("!!! ", ::akantu::dblError, info); \ AKANTU_CUSTOM_EXCEPTION_INFO(::akantu::debug::CriticalError(), info); \ } while (false) #endif // AKANTU_NDEBUG #define AKANTU_TO_IMPLEMENT() \ - AKANTU_CUSTOM_EXCEPTION_INFO(::akantu::debug::NotImplementedException(), \ - __func__ << " : not implemented yet !") + AKANTU_CUSTOM_EXCEPTION_INFO( \ + ::akantu::debug::NotImplementedException(), \ + __func__ \ + << " : not implemented yet !") // NOLINT(cppcoreguidelines-pro-bounds-array-to-pointer-decay) /* -------------------------------------------------------------------------- */ namespace debug { /* ------------------------------------------------------------------------ */ template void Debugger::throwCustomException(const Except & ex, const std::string & info, const std::string & file, unsigned int line, const std::string & module_) const noexcept(false) { auto & nc_ex = const_cast(ex); nc_ex.setInfo(info); nc_ex.setFile(file); nc_ex.setLine(line); nc_ex.setModule(module_); if (::akantu::debug::debugger.printBacktrace()) nc_ex.setBacktrace(::akantu::debug::getBacktrace()); throw ex; } /* ------------------------------------------------------------------------ */ template void Debugger::throwCustomException(const Except & ex, const std::string & file, unsigned int line, const std::string & module_) const noexcept(false) { auto & nc_ex = const_cast(ex); nc_ex.setFile(file); nc_ex.setLine(line); nc_ex.setModule(module_); if (::akantu::debug::debugger.printBacktrace()) nc_ex.setBacktrace(::akantu::debug::getBacktrace()); throw ex; } } // namespace debug } // namespace akantu #endif /* __AKANTU_ERROR_HH__ */ diff --git a/src/common/aka_extern.cc b/src/common/aka_extern.cc index d7d555a72..3f38bec5c 100644 --- a/src/common/aka_extern.cc +++ b/src/common/aka_extern.cc @@ -1,107 +1,103 @@ /** * @file aka_extern.cc * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Tue Feb 20 2018 * * @brief initialisation of all global variables * to insure the order of creation * * * 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 "aka_array.hh" #include "aka_common.hh" #include "aka_math.hh" #include "aka_named_argument.hh" #include "aka_random_generator.hh" #include "communication_tag.hh" #include "cppargparse.hh" #include "parser.hh" #include "solid_mechanics_model.hh" #if defined(AKANTU_COHESIVE_ELEMENT) #include "solid_mechanics_model_cohesive.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #endif namespace akantu { /* -------------------------------------------------------------------------- */ /* error.hpp variables */ /* -------------------------------------------------------------------------- */ namespace debug { /** \todo write function to get this * values from the environment or a config file */ /// standard output for debug messages std::ostream * _akantu_debug_cout = &std::cerr; /// standard output for normal messages std::ostream & _akantu_cout = std::cout; /// parallel context used in debug messages std::string _parallel_context = ""; Debugger debugger; #if defined(AKANTU_DEBUG_TOOLS) DebugElementManager element_manager; #endif } // namespace debug -/* -------------------------------------------------------------------------- */ -/// list of ghost iterable types -ghost_type_t ghost_types(_casper); - /* -------------------------------------------------------------------------- */ /// Paser for commandline arguments ::cppargparse::ArgumentParser static_argparser; /// Parser containing the information parsed by the input file given to initFull Parser static_parser; bool Parser::permissive_parser = false; /* -------------------------------------------------------------------------- */ Real Math::tolerance = 1e2 * std::numeric_limits::epsilon(); /* -------------------------------------------------------------------------- */ const UInt _all_dimensions [[gnu::unused]] = UInt(-1); /* -------------------------------------------------------------------------- */ const Array empty_filter(0, 1, "empty_filter"); /* -------------------------------------------------------------------------- */ template <> long int RandomGenerator::_seed = 5489u; template <> std::default_random_engine RandomGenerator::generator(5489u); /* -------------------------------------------------------------------------- */ int Tag::max_tag = 0; /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/common/aka_safe_enum.hh b/src/common/aka_safe_enum.hh index 267d5b719..630f14dc7 100644 --- a/src/common/aka_safe_enum.hh +++ b/src/common/aka_safe_enum.hh @@ -1,85 +1,94 @@ /** * @file aka_safe_enum.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Tue Nov 07 2017 * * @brief Safe enums type (see More C++ Idioms/Type Safe Enum on Wikibooks * http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Type_Safe_Enum) * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_SAFE_ENUM_HH__ #define __AKANTU_AKA_SAFE_ENUM_HH__ namespace akantu { /// Safe enumerated type template class safe_enum : public def { using type = typename def::type; public: - explicit safe_enum(type v = def::_end_) : val(v) {} - safe_enum(safe_enum && other) = default; - safe_enum & operator=(safe_enum && other) = default; + constexpr explicit safe_enum(type v = def::_end_) : val(v) {} - inner underlying() const { return val; } + constexpr inner underlying() const { return val; } - bool operator==(const safe_enum & s) const { return this->val == s.val; } - bool operator!=(const safe_enum & s) const { return this->val != s.val; } - bool operator<(const safe_enum & s) const { return this->val < s.val; } - bool operator<=(const safe_enum & s) const { return this->val <= s.val; } - bool operator>(const safe_enum & s) const { return this->val > s.val; } - bool operator>=(const safe_enum & s) const { return this->val >= s.val; } + constexpr bool operator==(const safe_enum & s) const { + return this->val == s.val; + } + constexpr bool operator!=(const safe_enum & s) const { + return this->val != s.val; + } + constexpr bool operator<(const safe_enum & s) const { + return this->val < s.val; + } + constexpr bool operator<=(const safe_enum & s) const { + return this->val <= s.val; + } + constexpr bool operator>(const safe_enum & s) const { + return this->val > s.val; + } + constexpr bool operator>=(const safe_enum & s) const { + return this->val >= s.val; + } - operator inner() { return val; }; + constexpr operator inner() { return val; }; public: // Works only if enumerations are contiguous. - class iterator { + class const_iterator { public: - explicit iterator(type v) : it(v) {} - iterator & operator++() { + constexpr explicit const_iterator(type v) : it(v) {} + constexpr const_iterator & operator++() { ++it; return *this; } - safe_enum operator*() { return safe_enum(static_cast(it)); } - bool operator!=(iterator const & it) { return it.it != this->it; } + constexpr safe_enum operator*() { return safe_enum(static_cast(it)); } + constexpr bool operator!=(const_iterator const & it) { return it.it != this->it; } private: int it; }; - static iterator begin() { return iterator(def::_begin_); } + constexpr auto begin() const { return const_iterator(def::_begin_); } + constexpr auto end() const { return const_iterator(def::_end_); } - static iterator end() { return iterator(def::_end_); } - -protected: +private: inner val; }; } // namespace akantu #endif /* __AKANTU_AKA_SAFE_ENUM_HH__ */ diff --git a/src/fe_engine/shape_functions.cc b/src/fe_engine/shape_functions.cc index 735ace8b7..b8224e57e 100644 --- a/src/fe_engine/shape_functions.cc +++ b/src/fe_engine/shape_functions.cc @@ -1,234 +1,236 @@ /** * @file shape_functions.cc * * @author Nicolas Richart * * @date creation: Wed Aug 09 2017 * @date last modification: Wed Oct 11 2017 * * @brief implementation of th shape functions interface * * * Copyright (©) 2016-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 "shape_functions.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ShapeFunctions::ShapeFunctions(const Mesh & mesh, UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), shapes("shapes_generic", id, memory_id), shapes_derivatives("shapes_derivatives_generic", id, memory_id), mesh(mesh), _spatial_dimension(spatial_dimension) {} /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::initElementalFieldInterpolationFromIntegrationPoints( const Array & interpolation_points_coordinates, ElementTypeMapArray & interpolation_points_coordinates_matrices, ElementTypeMapArray & quad_points_coordinates_inv_matrices, const Array & quadrature_points_coordinates, const GhostType & ghost_type, const Array & element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); UInt nb_element_filter; if (element_filter == empty_filter) nb_element_filter = nb_element; else nb_element_filter = element_filter.size(); UInt nb_quad_per_element = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_interpolation_points_per_elem = interpolation_points_coordinates.size() / nb_element; AKANTU_DEBUG_ASSERT(interpolation_points_coordinates.size() % nb_element == 0, "Number of interpolation points should be a multiple of " "total number of elements"); if (!quad_points_coordinates_inv_matrices.exists(type, ghost_type)) quad_points_coordinates_inv_matrices.alloc( nb_element_filter, nb_quad_per_element * nb_quad_per_element, type, ghost_type); else quad_points_coordinates_inv_matrices(type, ghost_type) .resize(nb_element_filter); if (!interpolation_points_coordinates_matrices.exists(type, ghost_type)) interpolation_points_coordinates_matrices.alloc( nb_element_filter, nb_interpolation_points_per_elem * nb_quad_per_element, type, ghost_type); else interpolation_points_coordinates_matrices(type, ghost_type) .resize(nb_element_filter); Array & quad_inv_mat = quad_points_coordinates_inv_matrices(type, ghost_type); Array & interp_points_mat = interpolation_points_coordinates_matrices(type, ghost_type); Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); Array::const_matrix_iterator quad_coords_it = quadrature_points_coordinates.begin_reinterpret( spatial_dimension, nb_quad_per_element, nb_element_filter); Array::const_matrix_iterator points_coords_begin = interpolation_points_coordinates.begin_reinterpret( spatial_dimension, nb_interpolation_points_per_elem, nb_element); Array::matrix_iterator inv_quad_coord_it = quad_inv_mat.begin(nb_quad_per_element, nb_quad_per_element); Array::matrix_iterator int_points_mat_it = interp_points_mat.begin( nb_interpolation_points_per_elem, nb_quad_per_element); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element_filter; ++el, ++inv_quad_coord_it, ++int_points_mat_it, ++quad_coords_it) { /// matrix containing the quadrature points coordinates const Matrix & quad_coords = *quad_coords_it; /// matrix to store the matrix inversion result Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /// insert the quad coordinates in a matrix compatible with the /// interpolation buildElementalFieldInterpolationMatrix(quad_coords, quad_coord_matrix); /// invert the interpolation matrix inv_quad_coord_matrix.inverse(quad_coord_matrix); /// matrix containing the interpolation points coordinates const Matrix & points_coords = points_coords_begin[element_filter(el)]; /// matrix to store the interpolation points coordinates /// compatible with these functions Matrix & inv_points_coord_matrix = *int_points_mat_it; /// insert the quad coordinates in a matrix compatible with the /// interpolation buildElementalFieldInterpolationMatrix(points_coords, inv_points_coord_matrix); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ShapeFunctions::initElementalFieldInterpolationFromIntegrationPoints( const ElementTypeMapArray & interpolation_points_coordinates, ElementTypeMapArray & interpolation_points_coordinates_matrices, ElementTypeMapArray & quad_points_coordinates_inv_matrices, const ElementTypeMapArray & quadrature_points_coordinates, const ElementTypeMapArray * element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); for (auto ghost_type : ghost_types) { auto types_iterable = mesh.elementTypes(spatial_dimension, ghost_type); if (element_filter) { types_iterable = element_filter->elementTypes(spatial_dimension, ghost_type); } for (auto type : types_iterable) { UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array * elem_filter; if (element_filter) elem_filter = &((*element_filter)(type, ghost_type)); else elem_filter = &(empty_filter); #define AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS(type) \ this->initElementalFieldInterpolationFromIntegrationPoints( \ interpolation_points_coordinates(type, ghost_type), \ interpolation_points_coordinates_matrices, \ quad_points_coordinates_inv_matrices, \ quadrature_points_coordinates(type, ghost_type), ghost_type, \ *elem_filter) AKANTU_BOOST_REGULAR_ELEMENT_SWITCH( AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS); #undef AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ShapeFunctions::interpolateElementalFieldFromIntegrationPoints( const ElementTypeMapArray & field, const ElementTypeMapArray & interpolation_points_coordinates_matrices, const ElementTypeMapArray & quad_points_coordinates_inv_matrices, ElementTypeMapArray & result, const GhostType & ghost_type, const ElementTypeMapArray * element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); auto types_iterable = mesh.elementTypes(spatial_dimension, ghost_type); - if (element_filter) { + if (element_filter != nullptr) { types_iterable = element_filter->elementTypes(spatial_dimension, ghost_type); } for (auto type : types_iterable) { UInt nb_element = mesh.getNbElement(type, ghost_type); - if (nb_element == 0) + if (nb_element == 0) { continue; - + } + const Array * elem_filter; - if (element_filter) + if (element_filter != nullptr) { elem_filter = &((*element_filter)(type, ghost_type)); - else + } else { elem_filter = &(empty_filter); + } #define AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS(type) \ interpolateElementalFieldFromIntegrationPoints( \ field(type, ghost_type), \ interpolation_points_coordinates_matrices(type, ghost_type), \ quad_points_coordinates_inv_matrices(type, ghost_type), result, \ ghost_type, *elem_filter) AKANTU_BOOST_REGULAR_ELEMENT_SWITCH( AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS); #undef AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/fe_engine/shape_functions_inline_impl.hh b/src/fe_engine/shape_functions_inline_impl.hh index 7160fa79a..aa724c4a4 100644 --- a/src/fe_engine/shape_functions_inline_impl.hh +++ b/src/fe_engine/shape_functions_inline_impl.hh @@ -1,400 +1,400 @@ /** * @file shape_functions_inline_impl.hh * * @author Guillaume Anciaux * @author Fabian Barras * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Oct 27 2010 * @date last modification: Tue Feb 20 2018 * * @brief ShapeFunctions inline implementation * * * 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 "fe_engine.hh" #include "shape_functions.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_FUNCTIONS_INLINE_IMPL_HH__ #define __AKANTU_SHAPE_FUNCTIONS_INLINE_IMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ inline const Array & ShapeFunctions::getShapes(const ElementType & el_type, const GhostType & ghost_type) const { return shapes(FEEngine::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ inline const Array & ShapeFunctions::getShapesDerivatives(const ElementType & el_type, const GhostType & ghost_type) const { return shapes_derivatives(FEEngine::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ inline UInt ShapeFunctions::getShapeSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt shape_size = 0; #define GET_SHAPE_SIZE(type) shape_size = ElementClass::getShapeSize() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SIZE); // , #undef GET_SHAPE_SIZE AKANTU_DEBUG_OUT(); return shape_size; } /* -------------------------------------------------------------------------- */ inline UInt ShapeFunctions::getShapeDerivativesSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt shape_derivatives_size = 0; #define GET_SHAPE_DERIVATIVES_SIZE(type) \ shape_derivatives_size = ElementClass::getShapeDerivativesSize() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_SIZE); // , #undef GET_SHAPE_DERIVATIVES_SIZE AKANTU_DEBUG_OUT(); return shape_derivatives_size; } /* -------------------------------------------------------------------------- */ template void ShapeFunctions::setIntegrationPointsByType(const Matrix & points, const GhostType & ghost_type) { if (not this->integration_points.exists(type, ghost_type)) this->integration_points(type, ghost_type).shallowCopy(points); } /* -------------------------------------------------------------------------- */ inline void ShapeFunctions::buildInterpolationMatrix(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { switch (integration_order) { case 1: { for (UInt i = 0; i < coordinates.cols(); ++i) coordMatrix(i, 0) = 1; break; } case 2: { UInt nb_quadrature_points = coordMatrix.cols(); for (UInt i = 0; i < coordinates.cols(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(j - 1, i); } break; } default: { AKANTU_TO_IMPLEMENT(); break; } } } /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::buildElementalFieldInterpolationMatrix( const Matrix &, Matrix &, UInt) const { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_segment_2>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_segment_3>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_triangle_3>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_triangle_6>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_tetrahedron_4>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_tetrahedron_10>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /** * @todo Write a more efficient interpolation for quadrangles by * dropping unnecessary quadrature points * */ /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_quadrangle_4>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { if (integration_order != ElementClassProperty<_quadrangle_4>::polynomial_degree) { AKANTU_TO_IMPLEMENT(); } else { for (UInt i = 0; i < coordinates.cols(); ++i) { Real x = coordinates(0, i); Real y = coordinates(1, i); coordMatrix(i, 0) = 1; coordMatrix(i, 1) = x; coordMatrix(i, 2) = y; coordMatrix(i, 3) = x * y; } } } /* -------------------------------------------------------------------------- */ template <> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_quadrangle_8>( const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { if (integration_order != ElementClassProperty<_quadrangle_8>::polynomial_degree) { AKANTU_TO_IMPLEMENT(); } else { for (UInt i = 0; i < coordinates.cols(); ++i) { // UInt j = 0; Real x = coordinates(0, i); Real y = coordinates(1, i); coordMatrix(i, 0) = 1; coordMatrix(i, 1) = x; coordMatrix(i, 2) = y; coordMatrix(i, 3) = x * y; // for (UInt e = 0; e <= 2; ++e) { // for (UInt n = 0; n <= 2; ++n) { // coordMatrix(i, j) = std::pow(x, e) * std::pow(y, n); // ++j; // } // } } } } /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::interpolateElementalFieldFromIntegrationPoints( const Array & field, const Array & interpolation_points_coordinates_matrices, const Array & quad_points_coordinates_inv_matrices, ElementTypeMapArray & result, const GhostType & ghost_type, const Array & element_filter) const { AKANTU_DEBUG_IN(); auto nb_element = this->mesh.getNbElement(type, ghost_type); auto nb_quad_per_element = GaussIntegrationElement::getNbQuadraturePoints(); auto nb_interpolation_points_per_elem = interpolation_points_coordinates_matrices.getNbComponent() / nb_quad_per_element; - if (!result.exists(type, ghost_type)) + if (not result.exists(type, ghost_type)) result.alloc(nb_element * nb_interpolation_points_per_elem, field.getNbComponent(), type, ghost_type); if (element_filter != empty_filter) nb_element = element_filter.size(); Matrix coefficients(nb_quad_per_element, field.getNbComponent()); auto & result_vec = result(type, ghost_type); auto field_it = field.begin_reinterpret(field.getNbComponent(), nb_quad_per_element, nb_element); auto interpolation_points_coordinates_it = interpolation_points_coordinates_matrices.begin( nb_interpolation_points_per_elem, nb_quad_per_element); auto result_begin = result_vec.begin_reinterpret( field.getNbComponent(), nb_interpolation_points_per_elem, result_vec.size() / nb_interpolation_points_per_elem); auto inv_quad_coord_it = quad_points_coordinates_inv_matrices.begin( nb_quad_per_element, nb_quad_per_element); /// loop over the elements of the current filter and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { /** * matrix containing the inversion of the quadrature points' * coordinates */ const auto & inv_quad_coord_matrix = *inv_quad_coord_it; /** * multiply it by the field values over quadrature points to get * the interpolation coefficients */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// matrix containing the points' coordinates const auto & coord = *interpolation_points_coordinates_it; /// multiply the coordinates matrix by the coefficients matrix and store the /// result Matrix res(result_begin[element_filter(el)]); res.mul(coefficients, coord); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::interpolateElementalFieldOnIntegrationPoints( const Array & u_el, Array & uq, const GhostType & ghost_type, const Array & shapes, const Array & filter_elements) const { auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_nodes_per_element = ElementClass::getShapeSize(); auto nb_points = shapes.size() / mesh.getNbElement(type, ghost_type); auto nb_degree_of_freedom = u_el.getNbComponent() / nb_nodes_per_element; Array::const_matrix_iterator N_it; Array * filtered_N = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filtered_N = new Array(0, shapes.getNbComponent()); FEEngine::filterElementalData(mesh, shapes, *filtered_N, type, ghost_type, filter_elements); N_it = filtered_N->begin_reinterpret(nb_nodes_per_element, nb_points, nb_element); } else { N_it = shapes.begin_reinterpret(nb_nodes_per_element, nb_points, nb_element); } uq.resize(nb_element * nb_points); auto u_it = u_el.begin(nb_degree_of_freedom, nb_nodes_per_element); auto inter_u_it = uq.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element); for (UInt el = 0; el < nb_element; ++el, ++N_it, ++u_it, ++inter_u_it) { const auto & u = *u_it; const auto & N = *N_it; auto & inter_u = *inter_u_it; inter_u.template mul(u, N); } delete filtered_N; } /* -------------------------------------------------------------------------- */ template void ShapeFunctions::gradientElementalFieldOnIntegrationPoints( const Array & u_el, Array & out_nablauq, const GhostType & ghost_type, const Array & shapes_derivatives, const Array & filter_elements) const { AKANTU_DEBUG_IN(); auto nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); auto nb_points = integration_points(type, ghost_type).cols(); auto element_dimension = ElementClass::getNaturalSpaceDimension(); auto nb_degree_of_freedom = u_el.getNbComponent() / nb_nodes_per_element; auto nb_element = mesh.getNbElement(type, ghost_type); Array::const_matrix_iterator B_it; Array * filtered_B = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filtered_B = new Array(0, shapes_derivatives.getNbComponent()); FEEngine::filterElementalData(mesh, shapes_derivatives, *filtered_B, type, ghost_type, filter_elements); B_it = filtered_B->begin(element_dimension, nb_nodes_per_element); } else { B_it = shapes_derivatives.begin(element_dimension, nb_nodes_per_element); } out_nablauq.resize(nb_element * nb_points); auto u_it = u_el.begin(nb_degree_of_freedom, nb_nodes_per_element); auto nabla_u_it = out_nablauq.begin(nb_degree_of_freedom, element_dimension); for (UInt el = 0; el < nb_element; ++el, ++u_it) { const auto & u = *u_it; for (UInt q = 0; q < nb_points; ++q, ++B_it, ++nabla_u_it) { const auto & B = *B_it; auto & nabla_u = *nabla_u_it; nabla_u.template mul(u, B); } } delete filtered_B; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_SHAPE_FUNCTIONS_INLINE_IMPL_HH__ */ diff --git a/src/mesh/element_type_map_tmpl.hh b/src/mesh/element_type_map_tmpl.hh index 88230ffb4..f4eee8d7b 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,790 +1,790 @@ /** * @file element_type_map_tmpl.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Tue Feb 20 2018 * * @brief implementation of template functions of the ElementTypeMap and * ElementTypeMapArray classes * * * 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 "aka_static_if.hh" #include "element_type_map.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "element_type_conversion.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template inline std::string ElementTypeMap::printType(const SupportType & type, const GhostType & ghost_type) { std::stringstream sstr; sstr << "(" << ghost_type << ":" << type << ")"; return sstr.str(); } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::exists( const SupportType & type, const GhostType & ghost_type) const { return this->getData(ghost_type).find(type) != this->getData(ghost_type).end(); } /* -------------------------------------------------------------------------- */ template inline const Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMap::printType(type, ghost_type) << " in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); return it->second; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) { return this->getData(ghost_type)[type]; } /* -------------------------------------------------------------------------- */ template template inline Stored & ElementTypeMap:: operator()(U && insertee, const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it != this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("Element of type " << ElementTypeMap::printType(type, ghost_type) << " already in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); } else { auto & data = this->getData(ghost_type); const auto & res = data.insert(std::make_pair(type, std::forward(insertee))); it = res.first; } return it->second; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ template inline const typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) const { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ /// Works only if stored is a pointer to a class with a printself method template void ElementTypeMap::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> [" << std::endl; - for (auto gt : ghost_types) { + for (auto && gt : ghost_types) { const DataMap & data = getData(gt); for (auto & pair : data) { stream << space << space << ElementTypeMap::printType(pair.first, gt) << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template ElementTypeMap::ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ template ElementTypeMap::~ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ /* ElementTypeMapArray */ /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::copy( const ElementTypeMapArray & other) { for (auto ghost_type : ghost_types) { for (auto type : this->elementTypes(_all_dimensions, ghost_type, _ek_not_defined)) { const auto & array_to_copy = other(type, ghost_type); auto & array = this->alloc(0, array_to_copy.getNbComponent(), type, ghost_type); array.copy(array_to_copy); } } } /* -------------------------------------------------------------------------- */ template ElementTypeMapArray::ElementTypeMapArray( const ElementTypeMapArray & other) : parent(), Memory(other.id + "_copy", other.memory_id), name(other.name + "_copy") { this->copy(other); } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray::alloc( UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value) { std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; Array * tmp; auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { auto id = this->id + ":" + std::to_string(type) + ghost_id; tmp = &(Memory::alloc(id, size, nb_component, default_value)); this->getData(ghost_type)[type] = tmp; } else { AKANTU_DEBUG_INFO( "The vector " << this->id << this->printType(type, ghost_type) << " already exists, it is resized instead of allocated."); tmp = it->second; it->second->resize(size); } return *tmp; } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value) { this->alloc(size, nb_component, type, _not_ghost, default_value); this->alloc(size, nb_component, type, _ghost, default_value); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::free() { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & pair : data) { dealloc(pair.second->getID()); } data.clear(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::clear() { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->clear(); } } } /* -------------------------------------------------------------------------- */ template template inline void ElementTypeMapArray::set(const ST & value) { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->set(value); } } } /* -------------------------------------------------------------------------- */ template inline const Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this const ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class(\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class (\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect) { auto it = this->getData(ghost_type).find(type); if (AKANTU_DEBUG_TEST(dblWarning) && it != this->getData(ghost_type).end() && it->second != &vect) { AKANTU_DEBUG_WARNING( "The Array " << this->printType(type, ghost_type) << " is already registred, this call can lead to a memory leak."); } this->getData(ghost_type)[type] = &(const_cast &>(vect)); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::onElementsRemoved( const ElementTypeMapArray & new_numbering) { for (auto gt : ghost_types) { for (auto & type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { auto support_type = convertType(type); if (this->exists(support_type, gt)) { const auto & renumbering = new_numbering(type, gt); if (renumbering.size() == 0) continue; auto & vect = this->operator()(support_type, gt); auto nb_component = vect.getNbComponent(); Array tmp(renumbering.size(), nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.size(); ++i) { UInt new_i = renumbering(i); if (new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i * nb_component, nb_component * sizeof(T)); ++new_size; } } tmp.resize(new_size); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; for (UInt g = _not_ghost; g <= _ghost; ++g) { auto gt = (GhostType)g; const DataMap & data = this->getData(gt); typename DataMap::const_iterator it; for (it = data.begin(); it != data.end(); ++it) { stream << space << space << ElementTypeMapArray::printType(it->first, gt) << " [" << std::endl; it->second->printself(stream, indent + 3); stream << space << space << " ]" << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /* SupportType Iterator */ /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek) : list_begin(list_begin), list_end(list_end), dim(dim), kind(ek) {} /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( const type_iterator & it) : list_begin(it.list_begin), list_end(it.list_end), dim(it.dim), kind(it.kind) {} /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator:: operator=(const type_iterator & it) { if (this != &it) { list_begin = it.list_begin; list_end = it.list_end; dim = it.dim; kind = it.kind; } return *this; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() const { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator::operator++() { ++list_begin; while ((list_begin != list_end) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(list_begin->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(list_begin->first))))) ++list_begin; return *this; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator ElementTypeMap::type_iterator::operator++(int) { type_iterator tmp(*this); operator++(); return tmp; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator==(const type_iterator & other) const { return this->list_begin == other.list_begin; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator!=(const type_iterator & other) const { return this->list_begin != other.list_begin; } /* -------------------------------------------------------------------------- */ template auto ElementTypeMap::ElementTypesIteratorHelper::begin() -> iterator { auto b = container.get().getData(ghost_type).begin(); auto e = container.get().getData(ghost_type).end(); // loop until the first valid type while ((b != e) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(b->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(b->first))))) ++b; return iterator(b, e, dim, kind); } template auto ElementTypeMap::ElementTypesIteratorHelper::end() -> iterator { auto e = container.get().getData(ghost_type).end(); return iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ template auto ElementTypeMap::elementTypesImpl( UInt dim, GhostType ghost_type, ElementKind kind) const -> ElementTypesIteratorHelper { return ElementTypesIteratorHelper(*this, dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ template template auto ElementTypeMap::elementTypesImpl( const use_named_args_t & unused, pack &&... _pack) const -> ElementTypesIteratorHelper { return ElementTypesIteratorHelper(*this, unused, _pack...); } /* -------------------------------------------------------------------------- */ template inline auto ElementTypeMap::firstType( UInt dim, GhostType ghost_type, ElementKind kind) const -> type_iterator { return elementTypes(dim, ghost_type, kind).begin(); } /* -------------------------------------------------------------------------- */ template inline auto ElementTypeMap::lastType( UInt dim, GhostType ghost_type, ElementKind kind) const -> type_iterator { typename DataMap::const_iterator e; e = getData(ghost_type).end(); return typename ElementTypeMap::type_iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const ElementTypeMap & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ class ElementTypeMapArrayInitializer { protected: using CompFunc = std::function; public: ElementTypeMapArrayInitializer( const CompFunc & comp_func, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined) : comp_func(comp_func), spatial_dimension(spatial_dimension), ghost_type(ghost_type), element_kind(element_kind) {} const GhostType & ghostType() const { return ghost_type; } virtual UInt nbComponent(const ElementType & type) const { return comp_func(type, ghostType()); } virtual bool isNodal() const { return false; } protected: CompFunc comp_func; UInt spatial_dimension; GhostType ghost_type; ElementKind element_kind; }; /* -------------------------------------------------------------------------- */ class MeshElementTypeMapArrayInitializer : public ElementTypeMapArrayInitializer { using CompFunc = ElementTypeMapArrayInitializer::CompFunc; public: MeshElementTypeMapArrayInitializer( const Mesh & mesh, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false) : MeshElementTypeMapArrayInitializer( mesh, [nb_component](const ElementType &, const GhostType &) -> UInt { return nb_component; }, spatial_dimension, ghost_type, element_kind, with_nb_element, with_nb_nodes_per_element) {} MeshElementTypeMapArrayInitializer( const Mesh & mesh, const CompFunc & comp_func, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false) : ElementTypeMapArrayInitializer(comp_func, spatial_dimension, ghost_type, element_kind), mesh(mesh), with_nb_element(with_nb_element), with_nb_nodes_per_element(with_nb_nodes_per_element) {} decltype(auto) elementTypes() const { return mesh.elementTypes(this->spatial_dimension, this->ghost_type, this->element_kind); } virtual UInt size(const ElementType & type) const { if (with_nb_element) return mesh.getNbElement(type, this->ghost_type); return 0; } UInt nbComponent(const ElementType & type) const override { UInt res = ElementTypeMapArrayInitializer::nbComponent(type); if (with_nb_nodes_per_element) return (res * mesh.getNbNodesPerElement(type)); return res; } bool isNodal() const override { return with_nb_nodes_per_element; } protected: const Mesh & mesh; bool with_nb_element; bool with_nb_nodes_per_element; }; /* -------------------------------------------------------------------------- */ class FEEngineElementTypeMapArrayInitializer : public MeshElementTypeMapArrayInitializer { public: FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, const ElementTypeMapArrayInitializer::CompFunc & nb_component, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); UInt size(const ElementType & type) const override; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; ElementTypesIteratorHelper elementTypes() const; protected: const FEEngine & fe_engine; }; /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const Func & f, const T & default_value, bool do_not_default) { this->is_nodal = f.isNodal(); auto ghost_type = f.ghostType(); for (auto & type : f.elementTypes()) { if (not this->exists(type, ghost_type)) if (do_not_default) { auto & array = this->alloc(0, f.nbComponent(type), type, ghost_type); array.resize(f.size(type)); } else { this->alloc(f.size(type), f.nbComponent(type), type, ghost_type, default_value); } else { auto & array = this->operator()(type, ghost_type); if (not do_not_default) array.resize(f.size(type), default_value); else array.resize(f.size(type)); } } } /* -------------------------------------------------------------------------- */ /** * All parameters are named optionals * \param _nb_component a functor giving the number of components per * (ElementType, GhostType) pair or a scalar giving a unique number of * components * regardless of type * \param _spatial_dimension a filter for the elements of a specific dimension * \param _element_kind filter with element kind (_ek_regular, _ek_structural, * ...) * \param _with_nb_element allocate the arrays with the number of elements for * each * type in the mesh * \param _with_nb_nodes_per_element multiply the number of components by the * number of nodes per element * \param _default_value default inital value * \param _do_not_default do not initialize the allocated arrays * \param _ghost_type filter a type of ghost */ template template void ElementTypeMapArray::initialize(const Mesh & mesh, pack &&... _pack) { GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); bool all_ghost_types = requested_ghost_type == _casper; for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; auto functor = MeshElementTypeMapArrayInitializer( mesh, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, mesh.getSpatialDimension()), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined), OPTIONAL_NAMED_ARG(with_nb_element, false), OPTIONAL_NAMED_ARG(with_nb_nodes_per_element, false)); this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ /** * All parameters are named optionals * \param _nb_component a functor giving the number of components per * (ElementType, GhostType) pair or a scalar giving a unique number of * components * regardless of type * \param _spatial_dimension a filter for the elements of a specific dimension * \param _element_kind filter with element kind (_ek_regular, _ek_structural, * ...) * \param _default_value default inital value * \param _do_not_default do not initialize the allocated arrays * \param _ghost_type filter a specific ghost type * \param _all_ghost_types get all ghost types */ template template void ElementTypeMapArray::initialize(const FEEngine & fe_engine, pack &&... _pack) { bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, true); GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _not_ghost); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; auto functor = FEEngineElementTypeMapArrayInitializer( fe_engine, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, UInt(-2)), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)); this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template inline T & ElementTypeMapArray:: operator()(const Element & element, UInt component) { return this->operator()(element.type, element.ghost_type)(element.element, component); } /* -------------------------------------------------------------------------- */ template inline const T & ElementTypeMapArray:: operator()(const Element & element, UInt component) const { return this->operator()(element.type, element.ghost_type)(element.element, component); } /* -------------------------------------------------------------------------- */ template UInt ElementTypeMapArray::sizeImpl( UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & kind) const { UInt size = 0; for (auto && type : this->elementTypes(spatial_dimension, ghost_type, kind)) { size += this->operator()(type, ghost_type).size(); } return size; } /* -------------------------------------------------------------------------- */ template template UInt ElementTypeMapArray::size(pack &&... _pack) const { UInt size = 0; bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, true); GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _not_ghost); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; size += sizeImpl(OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)); } return size; } } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ */ diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index eb97e5afa..20288a86b 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,627 +1,659 @@ /** * @file mesh.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief class handling meshes * * * 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 "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" #include "mesh_global_data_updater.hh" #include "mesh_io.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "mesh_utils_distribution.hh" #include "node_synchronizer.hh" #include "periodic_node_synchronizer.hh" /* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + #ifdef AKANTU_USE_IOHELPER #include "dumper_field.hh" #include "dumper_internal_material_field.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator) : Memory(id, memory_id), GroupManager(*this, id + ":group_manager", memory_id), MeshData("mesh_data", id, memory_id), connectivities("connectivities", id, memory_id), ghosts_counters("ghosts_counters", id, memory_id), normals("normals", id, memory_id), spatial_dimension(spatial_dimension), size(spatial_dimension, 0.), bbox(spatial_dimension), bbox_local(spatial_dimension), communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); this->nodes_flags = std::make_shared>(0, 1, NodeFlag::_normal, id + ":nodes_flags"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id, memory_id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const { barycenter.resize(getNbElement(type, ghost_type)); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type}, std::get<1>(data)); } } class FacetGlobalConnectivityAccessor : public DataAccessor { public: FacetGlobalConnectivityAccessor(Mesh & mesh) : global_connectivity("global_connectivity", "facet_connectivity_synchronizer") { global_connectivity.initialize( mesh, _spatial_dimension = _all_dimensions, _with_nb_element = true, _with_nb_nodes_per_element = true, _element_kind = _ek_regular); mesh.getGlobalConnectivity(global_connectivity); } UInt getNbData(const Array & elements, const SynchronizationTag & tag) const { UInt size = 0; if (tag == SynchronizationTag::_smmc_facets_conn) { UInt nb_nodes = Mesh::getNbNodesPerElementList(elements); size += nb_nodes * sizeof(UInt); } return size; } void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag == SynchronizationTag::_smmc_facets_conn) { for (const auto & element : elements) { auto & conns = global_connectivity(element.type, element.ghost_type); for (auto n : arange(conns.getNbComponent())) { buffer << conns(element.element, n); } } } } void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag == SynchronizationTag::_smmc_facets_conn) { for (const auto & element : elements) { auto & conns = global_connectivity(element.type, element.ghost_type); for (auto n : arange(conns.getNbComponent())) { buffer >> conns(element.element, n); } } } } AKANTU_GET_MACRO(GlobalConnectivity, (global_connectivity), decltype(auto)); protected: ElementTypeMapArray global_connectivity; }; /* -------------------------------------------------------------------------- */ Mesh & Mesh::initMeshFacets(const ID & id) { AKANTU_DEBUG_IN(); if (mesh_facets) { AKANTU_DEBUG_OUT(); return *mesh_facets; } mesh_facets = std::make_unique(spatial_dimension, this->nodes, getID() + ":" + id, getMemoryID()); mesh_facets->mesh_parent = this; mesh_facets->is_mesh_facets = true; mesh_facets->nodes_flags = this->nodes_flags; mesh_facets->nodes_global_ids = this->nodes_global_ids; MeshUtils::buildAllFacets(*this, *mesh_facets, 0); if (mesh.isDistributed()) { mesh_facets->is_distributed = true; mesh_facets->element_synchronizer = std::make_unique( *mesh_facets, mesh.getElementSynchronizer()); FacetGlobalConnectivityAccessor data_accessor(*mesh_facets); /// communicate mesh_facets->element_synchronizer->synchronizeOnce( data_accessor, SynchronizationTag::_smmc_facets_conn); /// flip facets MeshUtils::flipFacets(*mesh_facets, data_accessor.getGlobalConnectivity(), _ghost); } /// transfers the the mesh physical names to the mesh facets if (not this->hasData("physical_names")) { AKANTU_DEBUG_OUT(); return *mesh_facets; } auto & mesh_phys_data = this->getData("physical_names"); auto & phys_data = mesh_facets->getData("physical_names"); phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); ElementTypeMapArray barycenters(getID(), "temporary_barycenters"); barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); for (auto && ghost_type : ghost_types) { for (auto && type : barycenters.elementTypes(spatial_dimension - 1, ghost_type)) { mesh_facets->getBarycenters(barycenters(type, ghost_type), type, ghost_type); } } for_each_element( mesh, [&](auto && element) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); auto norm_barycenter = barycenter.norm(); auto tolerance = Math::getTolerance(); if (norm_barycenter > tolerance) tolerance *= norm_barycenter; const auto & element_to_facet = mesh_facets->getElementToSubelement( element.type, element.ghost_type); Vector barycenter_facet(spatial_dimension); auto range = enumerate(make_view( barycenters(element.type, element.ghost_type), spatial_dimension)); #ifndef AKANTU_NDEBUG auto min_dist = std::numeric_limits::max(); #endif // this is a spacial search coded the most inefficient way. auto facet = std::find_if(range.begin(), range.end(), [&](auto && data) { auto facet = std::get<0>(data); if (element_to_facet(facet)[1] == ElementNull) return false; auto norm_distance = barycenter.distance(std::get<1>(data)); #ifndef AKANTU_NDEBUG min_dist = std::min(min_dist, norm_distance); #endif return (norm_distance < tolerance); }); if (facet == range.end()) { AKANTU_DEBUG_INFO("The element " << element << " did not find its associated facet in the " "mesh_facets! Try to decrease math tolerance. " "The closest element was at a distance of " << min_dist); return; } // set physical name phys_data(Element{element.type, UInt(std::get<0>(*facet)), element.ghost_type}) = mesh_phys_data(element); }, _spatial_dimension = spatial_dimension - 1); mesh_facets->createGroupsFromMeshData("physical_names"); AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { AKANTU_DEBUG_ASSERT(not is_distributed, "You cannot read a mesh that is already distributed"); MeshIO mesh_io; mesh_io.read(filename, *this, mesh_io_type); auto types = this->elementTypes(spatial_dimension, _not_ghost, _ek_not_defined); auto it = types.begin(); auto last = types.end(); if (it == last) { AKANTU_DEBUG_WARNING( "The mesh contained in the file " << filename << " does not seem to be of the good dimension." << " No element of dimension " << spatial_dimension << " were read."); } this->makeReady(); } /* -------------------------------------------------------------------------- */ void Mesh::write(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO mesh_io; mesh_io.write(filename, *this, mesh_io_type); } /* -------------------------------------------------------------------------- */ void Mesh::makeReady() { this->nb_global_nodes = this->nodes->size(); this->computeBoundingBox(); this->nodes_flags->resize(nodes->size(), NodeFlag::_normal); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Mesh [" << std::endl; stream << space << " + id : " << getID() << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent + 2); stream << space << " + connectivities [" << std::endl; connectivities.printself(stream, indent + 2); stream << space << " ]" << std::endl; GroupManager::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox() { AKANTU_DEBUG_IN(); bbox_local.reset(); for (auto & pos : make_view(*nodes, spatial_dimension)) { // if(!isPureGhostNode(i)) bbox_local += pos; } if (this->is_distributed) { bbox = bbox_local.allSum(*communicator); } else { bbox = bbox_local; } size = bbox.size(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : global_connectivity.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_not_defined, _ghost_type = ghost_type)) { if (not connectivities.exists(type, ghost_type)) continue; auto & local_conn = connectivities(type, ghost_type); auto & g_connectivity = global_connectivity(type, ghost_type); UInt nb_nodes = local_conn.size() * local_conn.getNbComponent(); std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [&](UInt l) -> UInt { return this->getNodeGlobalId(l); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name, const std::string & group_name) { if (group_name == "all") return this->getDumper(dumper_name); else return element_groups[group_name]->getDumper(dumper_name); } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & arrays) { ElementTypeMap nb_data_per_elem; for (auto type : arrays.elementTypes()) { UInt nb_elements = this->getNbElement(type); auto & array = arrays(type); nb_data_per_elem(type) = array.getNbComponent() * array.size(); nb_data_per_elem(type) /= nb_elements; } return nb_data_per_elem; } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array); template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array); /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind) { std::shared_ptr field; ElementTypeMapArray * internal = nullptr; try { internal = &(this->getData(field_id)); } catch (...) { return nullptr; } ElementTypeMap nb_data_per_elem = this->getNbDataPerElem(*internal); field = this->createElementalField( *internal, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); return field; } template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); #endif /* -------------------------------------------------------------------------- */ void Mesh::distributeImpl( Communicator & communicator, std::function edge_weight_function [[gnu::unused]], std::function vertex_weight_function [[gnu::unused]]) { AKANTU_DEBUG_ASSERT(is_distributed == false, "This mesh is already distribute"); this->communicator = &communicator; this->element_synchronizer = std::make_unique( *this, this->getID() + ":element_synchronizer", this->getMemoryID(), true); this->node_synchronizer = std::make_unique( *this, this->getID() + ":node_synchronizer", this->getMemoryID(), true); Int psize = this->communicator->getNbProc(); if (psize > 1) { #ifdef AKANTU_USE_SCOTCH Int prank = this->communicator->whoAmI(); if (prank == 0) { MeshPartitionScotch partition(*this, spatial_dimension); partition.partitionate(psize, edge_weight_function, vertex_weight_function); MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition); } else { MeshUtilsDistribution::distributeMeshCentralized(*this, 0); } #else if (psize > 1) { AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool"); } #endif } // if (psize > 1) this->is_distributed = true; this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const Array & node_list, Array & elements) { for (const auto & node : node_list) for (const auto & element : *nodes_to_elements[node]) elements.push_back(element); } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements() { Element e; UInt nb_nodes = nodes->size(); for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) this->nodes_to_elements[n]->clear(); else this->nodes_to_elements[n] = std::make_unique>(); } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); Array::const_iterator> conn_it = connectivities(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) nodes_to_elements[conn(n)]->insert(e); } } } } /* -------------------------------------------------------------------------- */ std::tuple Mesh::updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { if (global_data_updater) return this->global_data_updater->updateData(nodes_event, elements_event); else { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } } /* -------------------------------------------------------------------------- */ void Mesh::registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->global_data_updater = std::move(global_data_updater); } /* -------------------------------------------------------------------------- */ void Mesh::eraseElements(const Array & elements) { ElementTypeMap last_element; RemovedElementsEvent event(*this, "new_numbering", AKANTU_CURRENT_FUNCTION); auto & remove_list = event.getList(); auto & new_numbering = event.getNewNumbering(); for (auto && el : elements) { if (el.ghost_type != _not_ghost) { auto & count = ghosts_counters(el); --count; - if (count > 0) + if (count > 0) { continue; + } } remove_list.push_back(el); - if (not last_element.exists(el.type, el.ghost_type)) { - UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); - last_element(nb_element - 1, el.type, el.ghost_type); + if (not new_numbering.exists(el.type, el.ghost_type)) { + auto nb_element = mesh.getNbElement(el.type, el.ghost_type); auto & numbering = new_numbering.alloc(nb_element, 1, el.type, el.ghost_type); for (auto && pair : enumerate(numbering)) { std::get<1>(pair) = std::get<0>(pair); } } - UInt & pos = last_element(el.type, el.ghost_type); - auto & numbering = new_numbering(el.type, el.ghost_type); - - numbering(el.element) = UInt(-1); - numbering(pos) = el.element; - --pos; + new_numbering(el) = UInt(-1); } + auto find_last_not_deleted = [](auto && array, Int start) -> Int { + do { + --start; + } while (start >= 0 and array[start] == UInt(-1)); + + return start; + }; + + auto find_first_deleted = [](auto && array, Int start) -> Int { + auto begin = array.begin(); + auto it = std::find_if(begin + start, array.end(), + [](auto & el) { return el == UInt(-1); }); + return Int(it - begin); + }; + + for (auto ghost_type : ghost_types) { + for (auto type : new_numbering.elementTypes(_ghost_type = ghost_type)) { + auto & numbering = new_numbering(type, ghost_type); + auto last_not_delete = + find_last_not_deleted(numbering, numbering.size()); + + if (last_not_delete < 0) { + continue; + } + + auto pos = find_first_deleted(numbering, 0); + + while (pos < last_not_delete) { + std::swap(numbering[pos], numbering[last_not_delete]); + last_not_delete = find_last_not_deleted(numbering, last_not_delete); + pos = find_first_deleted(numbering, pos + 1); + } + } + } this->sendEvent(event); } } // namespace akantu diff --git a/src/mesh_utils/mesh_partition.cc b/src/mesh_utils/mesh_partition.cc index 618dfd031..551f5e838 100644 --- a/src/mesh_utils/mesh_partition.cc +++ b/src/mesh_utils/mesh_partition.cc @@ -1,556 +1,557 @@ /** * @file mesh_partition.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Tue Aug 17 2010 * @date last modification: Wed Jan 24 2018 * * @brief implementation of common part of all partitioner * * * 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 "mesh_partition.hh" #include "aka_iterators.hh" #include "aka_types.hh" #include "mesh_accessor.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MeshPartition::MeshPartition(const Mesh & mesh, UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), mesh(mesh), spatial_dimension(spatial_dimension), partitions("partition", id, memory_id), ghost_partitions("ghost_partition", id, memory_id), ghost_partitions_offset("ghost_partition_offset", id, memory_id), saved_connectivity("saved_connectivity", id, memory_id) { AKANTU_DEBUG_IN(); UInt nb_total_element = 0; for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { linearized_offsets.push_back(std::make_pair(type, nb_total_element)); nb_total_element += mesh.getConnectivity(type).size(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MeshPartition::~MeshPartition() = default; /* -------------------------------------------------------------------------- */ UInt MeshPartition::linearized(const Element & element) { auto it = std::find_if(linearized_offsets.begin(), linearized_offsets.end(), [&element](auto & a) { return a.first == element.type; }); AKANTU_DEBUG_ASSERT(it != linearized_offsets.end(), "A bug might be crawling around this corner..."); return (it->second + element.element); } /* -------------------------------------------------------------------------- */ Element MeshPartition::unlinearized(UInt lin_element) { ElementType type{_not_defined}; UInt offset{0}; for (auto & pair : linearized_offsets) { if (lin_element < pair.second) continue; std::tie(type, offset) = pair; } return Element{type, lin_element - offset, _not_ghost}; } /* -------------------------------------------------------------------------- */ /** * conversion in c++ of the METIS_MeshToDual (mesh.c) function wrote by George * in Metis (University of Minnesota) */ void MeshPartition::buildDualGraph( Array & dxadj, Array & dadjncy, Array & edge_loads, std::function edge_load_func, Array & vertex_loads, std::function vertex_load_func) { CSR nodes_to_elements; MeshUtils::buildNode2Elements(mesh, nodes_to_elements); std::unordered_map> adjacent_elements; // for each elements look for its connected elements for_each_element( mesh, [&](auto && element) { const auto & conn = mesh.getConnectivity(element); std::map hits; // count the number of nodes shared with a given element for (auto && node : conn) { for (auto && connected_element : nodes_to_elements.getRow(node)) { ++hits[connected_element]; } } // define a minumum number of nodes to share to be considered as a // ajacent element UInt magic_number{conn.size()}; for (auto n : arange(mesh.getNbFacetTypes(element.type))) { magic_number = std::min( mesh.getNbNodesPerElement(mesh.getFacetType(element.type, n)), magic_number); } // check all neighbors to see which ones are "adjacent" for (auto && data : hits) { const auto & adjacent_element = data.first; // not adjacent to miself if (adjacent_element == element) continue; // not enough shared nodes if (data.second < magic_number) continue; /// Patch in order to prevent neighboring cohesive elements /// from detecting each other #if defined(AKANTU_COHESIVE_ELEMENT) auto element_kind = element.kind(); auto adjacent_element_kind = adjacent_element.kind(); if (element_kind == adjacent_element_kind && element_kind == _ek_cohesive) continue; #endif adjacent_elements[this->linearized(element)].push_back( this->linearized(adjacent_element)); } }, _spatial_dimension = mesh.getSpatialDimension(), _element_kind = _ek_not_defined); // prepare the arrays auto nb_elements{adjacent_elements.size()}; dxadj.resize(nb_elements + 1); vertex_loads.resize(nb_elements); for (auto && data : adjacent_elements) { const auto & element{data.first}; const auto & neighbors{data.second}; dxadj[element] = neighbors.size(); } /// convert the dxadj array of sizes in a csr one of offsets for (UInt i = 1; i < nb_elements; ++i) dxadj(i) += dxadj(i - 1); for (UInt i = nb_elements; i > 0; --i) dxadj(i) = dxadj(i - 1); dxadj(0) = 0; dadjncy.resize(dxadj(nb_elements)); edge_loads.resize(dadjncy.size()); // fill the different arrays for (auto && data : adjacent_elements) { const auto & element{data.first}; const auto & neighbors{data.second}; auto unlinearized_element = unlinearized(element); vertex_loads(element) = vertex_load_func(unlinearized_element); auto pos = dxadj(element); for (auto && neighbor : neighbors) { dadjncy(pos) = neighbor; edge_loads(pos) = edge_load_func(unlinearized_element, unlinearized(neighbor)); ++pos; } } } /* -------------------------------------------------------------------------- */ /** * TODO this function should most probably be rewritten in a more modern way * conversion in c++ of the GENDUALMETIS (mesh.c) function wrote by George in * Metis (University of Minnesota) */ // void MeshPartition::buildDualGraph( // Array & dxadj, Array & dadjncy, Array & edge_loads, // std::function edge_load_func, // Array & vertex_loads, // std::function vertex_load_func) { // AKANTU_DEBUG_IN(); // std::map *, UInt, UInt>> // connectivities; // UInt spatial_dimension = mesh.getSpatialDimension(); // UInt nb_total_element{0}; // for (auto & type : // mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { // auto type_p1 = mesh.getP1ElementType(type); // auto nb_nodes_per_element_p1 = mesh.getNbNodesPerElement(type_p1); // const auto & conn = mesh.getConnectivity(type, _not_ghost); // for (auto n : arange(mesh.getNbFacetTypes(type_p1))) { // auto magic_number = // mesh.getNbNodesPerElement(mesh.getFacetType(type_p1, n)); // connectivities[type] = // std::make_tuple(&conn, nb_nodes_per_element_p1, magic_number); // } // nb_total_element += conn.size(); // } // CSR node_to_elem; // MeshUtils::buildNode2Elements(mesh, node_to_elem); // dxadj.resize(nb_total_element + 1); // /// initialize the dxadj array // auto dxadj_it = dxadj.begin(); // for (auto & pair : connectivities) { // const auto & connectivity = *std::get<0>(pair.second); // auto nb_nodes_per_element_p1 = std::get<1>(pair.second); // std::fill_n(dxadj_it, connectivity.size(), nb_nodes_per_element_p1); // dxadj_it += connectivity.size(); // } // /// convert the dxadj_val array in a csr one // for (UInt i = 1; i < nb_total_element; ++i) // dxadj(i) += dxadj(i - 1); // for (UInt i = nb_total_element; i > 0; --i) // dxadj(i) = dxadj(i - 1); // dxadj(0) = 0; // dadjncy.resize(2 * dxadj(nb_total_element)); // /// weight map to determine adjacency // std::unordered_map weight_map; // for (auto & pair : connectivities) { // auto type = pair.first; // const auto & connectivity = *std::get<0>(pair.second); // auto nb_nodes_per_element = std::get<1>(pair.second); // auto magic_number = std::get<2>(pair.second); // Element element{type, 0, _not_ghost}; // for (const auto & conn : // make_view(connectivity, connectivity.getNbComponent())) { // auto linearized_el = linearized(element); // /// fill the weight map // for (UInt n : arange(nb_nodes_per_element)) { // auto && node = conn(n); // for (auto k = node_to_elem.rbegin(node); k != // node_to_elem.rend(node); // --k) { // auto & current_element = *k; // auto current_el = linearized(current_element); // AKANTU_DEBUG_ASSERT(current_el != UInt(-1), // "Linearized element not found"); // if (current_el <= linearized_el) // break; // auto weight_map_insert = // weight_map.insert(std::make_pair(current_el, 1)); // if (not weight_map_insert.second) // (weight_map_insert.first->second)++; // } // } // /// each element with a weight of the size of a facet are adjacent // for (auto & weight_pair : weight_map) { // auto & adjacent_el = weight_pair.first; // auto magic = weight_pair.second; // if (magic != magic_number) // continue; // #if defined(AKANTU_COHESIVE_ELEMENT) // /// Patch in order to prevent neighboring cohesive elements // /// from detecting each other // auto adjacent_element = unlinearized(adjacent_el); // auto el_kind = element.kind(); // auto adjacent_el_kind = adjacent_element.kind(); // if (el_kind == adjacent_el_kind && el_kind == _ek_cohesive) // continue; // #endif // UInt index_adj = dxadj(adjacent_el)++; // UInt index_lin = dxadj(linearized_el)++; // dadjncy(index_lin) = adjacent_el; // dadjncy(index_adj) = linearized_el; // } // element.element++; // weight_map.clear(); // } // } // Int k_start = 0, linerized_el = 0, j = 0; // for (auto & pair : connectivities) { // const auto & connectivity = *std::get<0>(pair.second); // auto nb_nodes_per_element_p1 = std::get<1>(pair.second); // auto nb_element = connectivity.size(); // for (UInt el = 0; el < nb_element; ++el, ++linerized_el) { // for (Int k = k_start; k < dxadj(linerized_el); ++k, ++j) // dadjncy(j) = dadjncy(k); // dxadj(linerized_el) = j; // k_start += nb_nodes_per_element_p1; // } // } // for (UInt i = nb_total_element; i > 0; --i) // dxadj(i) = dxadj(i - 1); // dxadj(0) = 0; // vertex_loads.resize(dxadj.size() - 1); // edge_loads.resize(dadjncy.size()); // UInt adj = 0; // for (UInt i = 0; i < nb_total_element; ++i) { // auto el = unlinearized(i); // vertex_loads(i) = vertex_load_func(el); // UInt nb_adj = dxadj(i + 1) - dxadj(i); // for (UInt j = 0; j < nb_adj; ++j, ++adj) { // auto el_adj_id = dadjncy(dxadj(i) + j); // auto el_adj = unlinearized(el_adj_id); // Int load = edge_load_func(el, el_adj); // edge_loads(adj) = load; // } // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void MeshPartition::fillPartitionInformation( const Mesh & mesh, const Int * linearized_partitions) { AKANTU_DEBUG_IN(); CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem); UInt linearized_el = 0; for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto & partition = partitions.alloc(nb_element, 1, type, _not_ghost); auto & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); auto & ghost_partition_offset = ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); auto & ghost_partition = ghost_partitions.alloc(0, 1, type, _ghost); const auto & connectivity = mesh.getConnectivity(type, _not_ghost); auto conn_it = connectivity.begin(connectivity.getNbComponent()); for (UInt el = 0; el < nb_element; ++el, ++linearized_el) { UInt part = linearized_partitions[linearized_el]; partition(el) = part; std::list list_adj_part; for (UInt n = 0; n < nb_nodes_per_element; ++n) { auto conn = Vector(*(conn_it + el)); UInt node = conn(n); for (const auto & adj_element : node_to_elem.getRow(node)) { UInt adj_el = linearized(adj_element); UInt adj_part = linearized_partitions[adj_el]; if (part != adj_part) { list_adj_part.push_back(adj_part); } } } list_adj_part.sort(); list_adj_part.unique(); for (auto & adj_part : list_adj_part) { ghost_part_csr.getRows().push_back(adj_part); ghost_part_csr.rowOffset(el)++; ghost_partition.push_back(adj_part); ghost_partition_offset(el)++; } } ghost_part_csr.countToCSR(); /// convert the ghost_partitions_offset array in an offset array auto & ghost_partitions_offset_ptr = ghost_partitions_offset(type, _ghost); for (UInt i = 1; i < nb_element; ++i) ghost_partitions_offset_ptr(i) += ghost_partitions_offset_ptr(i - 1); for (UInt i = nb_element; i > 0; --i) ghost_partitions_offset_ptr(i) = ghost_partitions_offset_ptr(i - 1); ghost_partitions_offset_ptr(0) = 0; } // All Facets for (Int sp = spatial_dimension - 1; sp >= 0; --sp) { for (auto & type : mesh.elementTypes(sp, _not_ghost, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type); auto & partition = partitions.alloc(nb_element, 1, type, _not_ghost); AKANTU_DEBUG_INFO("Allocating partitions for " << type); auto & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); auto & ghost_partition_offset = ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); auto & ghost_partition = ghost_partitions.alloc(0, 1, type, _ghost); AKANTU_DEBUG_INFO("Allocating ghost_partitions for " << type); const Array> & elem_to_subelem = mesh.getElementToSubelement(type, _not_ghost); // Facet loop for (UInt i(0); i < mesh.getNbElement(type, _not_ghost); ++i) { const auto & adjacent_elems = elem_to_subelem(i); if (adjacent_elems.empty()) { partition(i) = 0; continue; } Element min_elem{_max_element_type, std::numeric_limits::max(), - *ghost_type_t::end()}; + *(ghost_type_t{}.end())}; UInt min_part(std::numeric_limits::max()); std::set adjacent_parts; for (auto adj_elem : adjacent_elems) { - if (adj_elem == ElementNull) // case of boundary elements + if (adj_elem == ElementNull) { // case of boundary elements continue; + } auto adjacent_elem_part = partitions(adj_elem); if (adjacent_elem_part < min_part) { min_part = adjacent_elem_part; min_elem = adj_elem; } adjacent_parts.insert(adjacent_elem_part); } partition(i) = min_part; auto git = ghost_partitions_csr(min_elem.type, _not_ghost) .begin(min_elem.element); auto gend = ghost_partitions_csr(min_elem.type, _not_ghost) .end(min_elem.element); for (; git != gend; ++git) { adjacent_parts.insert(*git); } adjacent_parts.erase(min_part); for (auto & part : adjacent_parts) { ghost_part_csr.getRows().push_back(part); ghost_part_csr.rowOffset(i)++; ghost_partition.push_back(part); } ghost_partition_offset(i + 1) = ghost_partition_offset(i + 1) + adjacent_elems.size(); } ghost_part_csr.countToCSR(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::tweakConnectivity() { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(const_cast(mesh)); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { auto & connectivity = mesh_accessor.getConnectivity(type, _not_ghost); auto & saved_conn = saved_connectivity.alloc( connectivity.size(), connectivity.getNbComponent(), type, _not_ghost); saved_conn.copy(connectivity); for (auto && conn : make_view(connectivity, connectivity.getNbComponent())) { for (auto && node : conn) { if (mesh.isPeriodicSlave(node)) { node = mesh.getPeriodicMaster(node); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::restoreConnectivity() { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(const_cast(mesh)); for (auto && type : saved_connectivity.elementTypes( spatial_dimension, _not_ghost, _ek_not_defined)) { auto & conn = mesh_accessor.getConnectivity(type, _not_ghost); auto & saved_conn = saved_connectivity(type, _not_ghost); conn.copy(saved_conn); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool MeshPartition::hasPartitions(const ElementType & type, const GhostType & ghost_type) { return partitions.exists(type, ghost_type); } /* -------------------------------------------------------------------------- */ void MeshPartition::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "MeshPartition [" << "\n"; stream << space << " + id : " << id << "\n"; stream << space << " + nb partitions: " << nb_partitions << "\n"; stream << space << " + partitions [ " << "\n"; partitions.printself(stream, indent + 2); stream << space << " ]" << "\n"; stream << space << "]" << "\n"; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/communications_tmpl.hh b/src/synchronizer/communications_tmpl.hh index 49d2471ac..1d6447403 100644 --- a/src/synchronizer/communications_tmpl.hh +++ b/src/synchronizer/communications_tmpl.hh @@ -1,550 +1,550 @@ /** * @file communications_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Sep 07 2016 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of Communications * * * Copyright (©) 2016-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 "communications.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMUNICATIONS_TMPL_HH__ #define __AKANTU_COMMUNICATIONS_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template Communications::Communications(const Communicator & communicator) : communicator(communicator) {} /* -------------------------------------------------------------------------- */ template Communications::Communications(const Communications & other) : communicator(other.communicator) { for (auto sr : iterate_send_recv) { for (const auto & scheme_pair : other.iterateSchemes(sr)) { auto proc = scheme_pair.first; auto & other_scheme = scheme_pair.second; auto & scheme = this->createScheme(proc, sr); scheme.copy(other_scheme); } } this->invalidateSizes(); } /* -------------------------------------------------------------------------- */ template void Communications::swapSendRecv() { std::swap(schemes[_send], schemes[_recv]); } /* -------------------------------------------------------------------------- */ template class Communications::iterator { using communication_iterator = typename std::map::iterator; public: iterator() : communications(nullptr) {} iterator(scheme_iterator scheme_it, communication_iterator comm_it, Communications & communications, const SynchronizationTag & tag) : scheme_it(scheme_it), comm_it(comm_it), communications(&communications), tag(tag) {} iterator(const iterator & other) = default; iterator(iterator && other) noexcept = default; iterator & operator=(const iterator & other) = default; iterator & operator=(iterator && other) noexcept = default; iterator & operator++() { ++scheme_it; ++comm_it; return *this; } CommunicationDescriptor operator*() { AKANTU_DEBUG_ASSERT( scheme_it->first == comm_it->first, "The two iterators are not in phase, something wrong" << " happened, time to take out your favorite debugger (" << scheme_it->first << " != " << comm_it->first << ")"); return CommunicationDescriptor(comm_it->second, scheme_it->second, *communications, tag, scheme_it->first); } bool operator==(const iterator & other) const { return scheme_it == other.scheme_it && comm_it == other.comm_it; } bool operator!=(const iterator & other) const { return scheme_it != other.scheme_it || comm_it != other.comm_it; } private: scheme_iterator scheme_it; communication_iterator comm_it; Communications * communications; SynchronizationTag tag; }; /* -------------------------------------------------------------------------- */ template class Communications::tag_iterator { using internal_iterator = std::map::const_iterator; public: tag_iterator(const internal_iterator & it) : it(it) {} tag_iterator & operator++() { ++it; return *this; } SynchronizationTag operator*() { return it->first; } bool operator==(const tag_iterator & other) const { return it == other.it; } bool operator!=(const tag_iterator & other) const { return it != other.it; } private: internal_iterator it; }; /* -------------------------------------------------------------------------- */ template typename Communications::CommunicationPerProcs & Communications::getCommunications(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto comm_it = this->communications[sr].find(tag); if (comm_it == this->communications[sr].end()) AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "No known communications for the tag: " << tag); return comm_it->second; } /* ---------------------------------------------------------------------- */ template UInt Communications::getPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) const { const std::map & pending = pending_communications[sr]; auto it = pending.find(tag); if (it == pending.end()) return 0; return it->second; } /* -------------------------------------------------------------------------- */ template bool Communications::hasPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) const { return this->hasCommunication(tag) && (this->getPending(tag, sr) != 0); } /* ---------------------------------------------------------------------- */ template typename Communications::iterator Communications::begin(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); return iterator(this->schemes[sr].begin(), comms.begin(), *this, tag); } template typename Communications::iterator Communications::end(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); return iterator(this->schemes[sr].end(), comms.end(), *this, tag); } /* ---------------------------------------------------------------------- */ template typename Communications::iterator Communications::waitAny(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); auto it = comms.begin(); auto end = comms.end(); std::vector requests; for (; it != end; ++it) { auto & request = it->second.request(); if (!request.isFreed()) requests.push_back(request); } UInt req_id = communicator.waitAny(requests); if (req_id != UInt(-1)) { auto & request = requests[req_id]; UInt proc = sr == _recv ? request.getSource() : request.getDestination(); return iterator(this->schemes[sr].find(proc), comms.find(proc), *this, tag); } else { return this->end(tag, sr); } } /* ---------------------------------------------------------------------- */ template void Communications::waitAll(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { auto & comms = this->getCommunications(tag, sr); auto it = comms.begin(); auto end = comms.end(); std::vector requests; for (; it != end; ++it) { requests.push_back(it->second.request()); } communicator.waitAll(requests); } template void Communications::incrementPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) { ++(pending_communications[sr][tag]); } template void Communications::decrementPending( const SynchronizationTag & tag, const CommunicationSendRecv & sr) { --(pending_communications[sr][tag]); } template void Communications::freeRequests(const SynchronizationTag & tag, const CommunicationSendRecv & sr) { iterator it = this->begin(tag, sr); iterator end = this->end(tag, sr); for (; it != end; ++it) { (*it).freeRequest(); } } /* -------------------------------------------------------------------------- */ template typename Communications::Scheme & Communications::createScheme(UInt proc, const CommunicationSendRecv & sr) { // scheme_iterator it = schemes[sr].find(proc); // if (it != schemes[sr].end()) { // AKANTU_CUSTOM_EXCEPTION_INFO(debug::CommunicationException(), // "Communication scheme(" // << sr // << ") already created for proc: " << // proc); // } return schemes[sr][proc]; } template void Communications::resetSchemes(const CommunicationSendRecv & sr) { auto it = this->schemes[sr].begin(); auto end = this->schemes[sr].end(); for (; it != end; ++it) { it->second.resize(0); } } /* -------------------------------------------------------------------------- */ template void Communications::setCommunicationSize( const SynchronizationTag & tag, UInt proc, UInt size, const CommunicationSendRecv & sr) { // accessor that fails if it does not exists comm_size_computed[tag] = true; // TODO: need perhaps to be split based on sr auto & comms = this->communications[sr]; auto & comms_per_tag = comms.at(tag); comms_per_tag.at(proc).resize(size); } /* -------------------------------------------------------------------------- */ template void Communications::initializeCommunications( const SynchronizationTag & tag) { - for (auto t = send_recv_t::begin(); t != send_recv_t::end(); ++t) { - pending_communications[*t].insert(std::make_pair(tag, 0)); + for (auto t : send_recv_t{}) { + pending_communications[t].insert(std::make_pair(tag, 0)); - auto & comms = this->communications[*t]; + auto & comms = this->communications[t]; auto & comms_per_tag = comms.insert(std::make_pair(tag, CommunicationPerProcs())) .first->second; - for (auto pair : this->schemes[*t]) { + for (auto pair : this->schemes[t]) { comms_per_tag.emplace(std::piecewise_construct, std::forward_as_tuple(pair.first), - std::forward_as_tuple(*t)); + std::forward_as_tuple(t)); } } comm_counter.insert(std::make_pair(tag, 0)); } /* -------------------------------------------------------------------------- */ template typename Communications::tag_iterator Communications::begin_tag() { return tag_iterator(comm_counter.begin()); } template typename Communications::tag_iterator Communications::end_tag() { return tag_iterator(comm_counter.end()); } /* -------------------------------------------------------------------------- */ template typename Communications::scheme_iterator Communications::begin_scheme(const CommunicationSendRecv & sr) { return this->schemes[sr].begin(); } template typename Communications::scheme_iterator Communications::end_scheme(const CommunicationSendRecv & sr) { return this->schemes[sr].end(); } /* -------------------------------------------------------------------------- */ template typename Communications::const_scheme_iterator Communications::begin_scheme(const CommunicationSendRecv & sr) const { return this->schemes[sr].begin(); } template typename Communications::const_scheme_iterator Communications::end_scheme(const CommunicationSendRecv & sr) const { return this->schemes[sr].end(); } /* -------------------------------------------------------------------------- */ template typename Communications::scheme_iterator Communications::begin_send_scheme() { return this->begin_scheme(_send); } template typename Communications::scheme_iterator Communications::end_send_scheme() { return this->end_scheme(_send); } /* -------------------------------------------------------------------------- */ template typename Communications::const_scheme_iterator Communications::begin_send_scheme() const { return this->begin_scheme(_send); } template typename Communications::const_scheme_iterator Communications::end_send_scheme() const { return this->end_scheme(_send); } /* -------------------------------------------------------------------------- */ template typename Communications::scheme_iterator Communications::begin_recv_scheme() { return this->begin_scheme(_recv); } template typename Communications::scheme_iterator Communications::end_recv_scheme() { return this->end_scheme(_recv); } /* -------------------------------------------------------------------------- */ template typename Communications::const_scheme_iterator Communications::begin_recv_scheme() const { return this->begin_scheme(_recv); } template typename Communications::const_scheme_iterator Communications::end_recv_scheme() const { return this->end_scheme(_recv); } /* ------------------------------------------------------------------------ */ template bool Communications::hasCommunication( const SynchronizationTag & tag) const { return (communications[_send].find(tag) != communications[_send].end()); } template void Communications::incrementCounter(const SynchronizationTag & tag) { auto it = comm_counter.find(tag); if (it == comm_counter.end()) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "No counter initialized in communications for the tags: " << tag); } ++(it->second); } template UInt Communications::getCounter(const SynchronizationTag & tag) const { auto it = comm_counter.find(tag); if (it == comm_counter.end()) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "No counter initialized in communications for the tags: " << tag); } return it->second; } template bool Communications::hasCommunicationSize( const SynchronizationTag & tag) const { auto it = comm_size_computed.find(tag); if (it == comm_size_computed.end()) { return false; } return it->second; } template void Communications::invalidateSizes() { for (auto && pair : comm_size_computed) { pair.second = false; } } template bool Communications::hasPendingRecv( const SynchronizationTag & tag) const { return this->hasPending(tag, _recv); } template bool Communications::hasPendingSend( const SynchronizationTag & tag) const { return this->hasPending(tag, _send); } template const auto & Communications::getCommunicator() const { return communicator; } /* -------------------------------------------------------------------------- */ template typename Communications::iterator Communications::waitAnyRecv(const SynchronizationTag & tag) { return this->waitAny(tag, _recv); } template typename Communications::iterator Communications::waitAnySend(const SynchronizationTag & tag) { return this->waitAny(tag, _send); } template void Communications::waitAllRecv(const SynchronizationTag & tag) { this->waitAll(tag, _recv); } template void Communications::waitAllSend(const SynchronizationTag & tag) { this->waitAll(tag, _send); } template void Communications::freeSendRequests(const SynchronizationTag & tag) { this->freeRequests(tag, _send); } template void Communications::freeRecvRequests(const SynchronizationTag & tag) { this->freeRequests(tag, _recv); } /* -------------------------------------------------------------------------- */ template typename Communications::Scheme & Communications::createSendScheme(UInt proc) { return createScheme(proc, _send); } template typename Communications::Scheme & Communications::createRecvScheme(UInt proc) { return createScheme(proc, _recv); } /* -------------------------------------------------------------------------- */ template void Communications::resetSchemes() { resetSchemes(_send); resetSchemes(_recv); } /* -------------------------------------------------------------------------- */ template typename Communications::Scheme & Communications::getScheme(UInt proc, const CommunicationSendRecv & sr) { return this->schemes[sr].find(proc)->second; } /* -------------------------------------------------------------------------- */ template const typename Communications::Scheme & Communications::getScheme(UInt proc, const CommunicationSendRecv & sr) const { return this->schemes[sr].find(proc)->second; } /* -------------------------------------------------------------------------- */ template void Communications::setSendCommunicationSize( const SynchronizationTag & tag, UInt proc, UInt size) { this->setCommunicationSize(tag, proc, size, _send); } template void Communications::setRecvCommunicationSize( const SynchronizationTag & tag, UInt proc, UInt size) { this->setCommunicationSize(tag, proc, size, _recv); } } // namespace akantu #endif /* __AKANTU_COMMUNICATIONS_TMPL_HH__ */ diff --git a/src/synchronizer/dof_synchronizer.cc b/src/synchronizer/dof_synchronizer.cc index 15b020453..b76da0414 100644 --- a/src/synchronizer/dof_synchronizer.cc +++ b/src/synchronizer/dof_synchronizer.cc @@ -1,228 +1,230 @@ /** * @file dof_synchronizer.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Tue Feb 06 2018 * * @brief DOF synchronizing object implementation * * * 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 "dof_synchronizer.hh" #include "aka_iterators.hh" #include "dof_manager_default.hh" #include "mesh.hh" #include "node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A DOFSynchronizer needs a mesh and the number of degrees of freedom * per node to be created. In the constructor computes the local and global dof * number for each dof. The member * proc_informations (std vector) is resized with the number of mpi * processes. Each entry in the vector is a PerProcInformations object * that contains the interactions of the current mpi process (prank) with the * mpi process corresponding to the position of that entry. Every * ProcInformations object contains one array with the dofs that have * to be sent to prank and a second one with dofs that willl be received form * prank. * This information is needed for the asychronous communications. The * constructor sets up this information. */ DOFSynchronizer::DOFSynchronizer(DOFManagerDefault & dof_manager, const ID & id, MemoryID memory_id) : SynchronizerImpl(dof_manager.getCommunicator(), id, memory_id), dof_manager(dof_manager) { std::vector dof_ids = dof_manager.getDOFIDs(); // Transfers nodes to global equation numbers in new schemes for (const ID & dof_id : dof_ids) { registerDOFs(dof_id); } } /* -------------------------------------------------------------------------- */ DOFSynchronizer::~DOFSynchronizer() = default; /* -------------------------------------------------------------------------- */ void DOFSynchronizer::registerDOFs(const ID & dof_id) { - if (this->nb_proc == 1) + if (this->nb_proc == 1) { return; + } - if (dof_manager.getSupportType(dof_id) != _dst_nodal) + if (dof_manager.getSupportType(dof_id) != _dst_nodal) { return; + } const auto & equation_numbers = dof_manager.getLocalEquationsNumbers(dof_id); const auto & associated_nodes = dof_manager.getDOFsAssociatedNodes(dof_id); const auto & node_synchronizer = dof_manager.getMesh().getNodeSynchronizer(); const auto & node_communications = node_synchronizer.getCommunications(); auto transcode_node_to_global_dof_scheme = [this, &associated_nodes, &equation_numbers]( auto && it, auto && end, const CommunicationSendRecv & sr) -> void { for (; it != end; ++it) { auto & scheme = communications.createScheme(it->first, sr); const auto & node_scheme = it->second; for (auto & node : node_scheme) { auto an_begin = associated_nodes.begin(); auto an_it = an_begin; auto an_end = associated_nodes.end(); std::vector global_dofs_per_node; while ((an_it = std::find(an_it, an_end, node)) != an_end) { UInt pos = an_it - an_begin; UInt local_eq_num = equation_numbers(pos); UInt global_eq_num = dof_manager.localToGlobalEquationNumber(local_eq_num); global_dofs_per_node.push_back(global_eq_num); ++an_it; } std::sort(global_dofs_per_node.begin(), global_dofs_per_node.end()); std::transform(global_dofs_per_node.begin(), global_dofs_per_node.end(), global_dofs_per_node.begin(), [this](UInt g) -> UInt { UInt l = dof_manager.globalToLocalEquationNumber(g); return l; }); for (auto & leqnum : global_dofs_per_node) { scheme.push_back(leqnum); } } } }; - for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); - ++sr_it) { - auto ncs_it = node_communications.begin_scheme(*sr_it); - auto ncs_end = node_communications.end_scheme(*sr_it); + for (auto sr : send_recv_t{}) { + auto ncs_it = node_communications.begin_scheme(sr); + auto ncs_end = node_communications.end_scheme(sr); - transcode_node_to_global_dof_scheme(ncs_it, ncs_end, *sr_it); + transcode_node_to_global_dof_scheme(ncs_it, ncs_end, sr); } entities_changed = true; } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::fillEntityToSend(Array & dofs_to_send) { UInt nb_dofs = dof_manager.getLocalSystemSize(); this->entities_from_root.clear(); dofs_to_send.resize(0); for (UInt d : arange(nb_dofs)) { - if (not dof_manager.isLocalOrMasterDOF(d)) + if (not dof_manager.isLocalOrMasterDOF(d)) { continue; + } entities_from_root.push_back(d); } for (auto d : entities_from_root) { UInt global_dof = dof_manager.localToGlobalEquationNumber(d); dofs_to_send.push_back(global_dof); } } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::onNodesAdded(const Array & /*nodes_list*/) { auto dof_ids = dof_manager.getDOFIDs(); for (auto sr : iterate_send_recv) { for (auto && data : communications.iterateSchemes(sr)) { auto & scheme = data.second; scheme.resize(0); } } for (auto & dof_id : dof_ids) { registerDOFs(dof_id); } // const auto & node_synchronizer = // dof_manager.getMesh().getNodeSynchronizer(); const auto & // node_communications = node_synchronizer.getCommunications(); // std::map> nodes_per_proc[2]; // for (auto sr : iterate_send_recv) { // for (auto && data : node_communications.iterateSchemes(sr)) { // auto proc = data.first; // const auto & scheme = data.second; // for (auto node : scheme) { // nodes_per_proc[sr][proc].push_back(node); // } // } // } // std::map> dofs_per_proc[2]; // for (auto & dof_id : dof_ids) { // const auto & associated_nodes = // dof_manager.getDOFsAssociatedNodes(dof_id); const auto & // local_equation_numbers = // dof_manager.getEquationsNumbers(dof_id); // for (auto tuple : zip(associated_nodes, local_equation_numbers)) { // UInt assoc_node; // UInt local_eq_num; // std::tie(assoc_node, local_eq_num) = tuple; // for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); // ++sr_it) { // for (auto & pair : nodes_per_proc[*sr_it]) { // if (std::find(pair.second.end(), pair.second.end(), assoc_node) != // pair.second.end()) { // dofs_per_proc[*sr_it][pair.first].push_back(local_eq_num); // } // } // } // } // } // for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); // ++sr_it) { // for (auto & pair : dofs_per_proc[*sr_it]) { // std::sort(pair.second.begin(), pair.second.end(), // [this](UInt la, UInt lb) -> bool { // UInt ga = dof_manager.localToGlobalEquationNumber(la); // UInt gb = dof_manager.localToGlobalEquationNumber(lb); // return ga < gb; // }); // auto & scheme = communications.getScheme(pair.first, *sr_it); // scheme.resize(0); // for (auto leq : pair.second) { // scheme.push_back(leq); // } // } // } this->entities_changed = true; } } // namespace akantu diff --git a/test/test_model/test_common/test_non_local_toolbox/plate.geo b/test/test_model/test_common/test_non_local_toolbox/plate.geo index f83107ac4..970dcd272 100644 --- a/test/test_model/test_common/test_non_local_toolbox/plate.geo +++ b/test/test_model/test_common/test_non_local_toolbox/plate.geo @@ -1,15 +1,15 @@ -h = 0.5; +h = .5; length = 2; Point(1) = {-length/2, -length/2, 0, h}; Point(2) = {0, -length/2, 0, h}; line1[] = Extrude{length/2, 0, 0}{ Point{1}; }; line2[] = Extrude{length/2, 0, 0}{ Point{2}; }; surface1[] = Extrude{0, length, 0}{ Line{line1[1]}; }; surface2[] = Extrude{0, length, 0}{ Line{line2[1]}; }; Transfinite Surface "*"; Recombine Surface "*"; Physical Surface("mat_1") = {surface1[1]}; Physical Surface("mat_2") = {surface2[1]};