diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index dcd1d2b3c..7a0bdf1e0 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,519 +1,519 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ using MemoryID = UInt; // using Surface = std::string; // using SurfacePair= std::pair; // using SurfacePairList = std::list; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; #define AKANTU_PP_ENUM(s, data, i, elem) \ BOOST_PP_TUPLE_REM() \ elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data))) } // namespace akantu #if (defined(__GNUC__) || defined(__GNUG__)) #define AKA_GCC_VERSION \ (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) #if AKA_GCC_VERSION < 60000 #define AKANTU_ENUM_HASH(type_name) \ namespace std { \ template <> struct hash<::akantu::type_name> { \ using argument_type = ::akantu::type_name; \ size_t operator()(const argument_type & e) const noexcept { \ auto ue = underlying_type_t(e); \ return uh(ue); \ } \ \ private: \ const hash> uh{}; \ }; \ } #else #define AKANTU_ENUM_HASH(type_name) #endif // AKA_GCC_VERSION #endif // GNU #include "aka_element_classes_info.hh" namespace akantu { #define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) #define AKANTU_PP_TYPE_TO_STR(s, data, elem) \ ({BOOST_PP_CAT(data::_, elem), BOOST_PP_STRINGIZE(elem)}) #define AKANTU_PP_STR_TO_TYPE(s, data, elem) \ ({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data::_, elem)}) #define AKANTU_ENUM_DECLARE(type_name, list) \ enum class type_name { \ BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \ }; #define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ } \ AKANTU_ENUM_HASH(type_name) \ namespace aka { \ inline std::string to_string(const ::akantu::type_name & type) { \ static std::unordered_map<::akantu::type_name, std::string> convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, \ ::akantu::type_name, list))}; \ return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << aka::to_string(type); \ return stream; \ } #define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ inline std::istream & operator>>(std::istream & stream, type_name & type) { \ std::string str; \ stream >> str; \ static std::unordered_map convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, type_name, list))}; \ type = convert.at(str); \ return stream; \ } /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpacialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #ifndef SWIG // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { _model, _solid_mechanics_model, _solid_mechanics_model_cohesive, _heat_transfer_model, _structural_mechanics_model, _embedded_model }; #endif /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent _nls_lumped, ///< Case of lumped mass or equivalent matrix _nls_auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; /// Define the node/dof type enum NodeType : Int { _nt_pure_ghost = -3, _nt_master = -2, _nt_normal = -1 }; /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_ask_nodes, _gst_size, //--- SolidMechanicsModel tags --- _gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _gst_smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _gst_smm_uv, ///< synchronization of the nodal velocities and displacement _gst_smm_res, ///< synchronization of the nodal residual _gst_smm_init_mat, ///< synchronization of the data to initialize materials _gst_smm_stress, ///< synchronization of the stresses to compute the internal /// forces _gst_smmc_facets, ///< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, ///< synchronization of facet global connectivity _gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet /// synch _gst_smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _gst_giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _gst_ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gst_gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _gst_htm_temperature, ///< synchronization of the nodal temperature _gst_htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _gst_htm_phi, ///< synchronization of the nodal level set value phi _gst_htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _gst_mnl_for_average, ///< synchronization of data to average in non local /// material _gst_mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _gst_nh_criterion, // --- General tags --- _gst_test, ///< Test tag _gst_user_1, ///< tag for user simulations _gst_user_2, ///< tag for user simulations _gst_material_id, ///< synchronization of the material ids _gst_for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _gst_cf_nodal, ///< synchronization of disp, velo, and current position _gst_cf_incr, ///< synchronization of increment // --- Solver tags --- _gst_solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; -} +} // namespace akantu #ifndef SWIG AKANTU_ENUM_HASH(GhostType) #endif namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ struct TensorTrait {}; /* -------------------------------------------------------------------------- */ template using is_tensor = std::is_base_of; /* -------------------------------------------------------------------------- */ template using is_scalar = std::is_arithmetic; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ // #if defined(__INTEL_COMPILER) // /// remark #981: operands are evaluated in unspecified order // #pragma warning(disable : 981) // /// remark #383: value copied to temporary, reference to temporary used // #pragma warning(disable : 383) // #endif // defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ } // namespace akantu #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic number * is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/common/aka_iterators.hh b/src/common/aka_iterators.hh index b9a33e8b1..78a7496bd 100644 --- a/src/common/aka_iterators.hh +++ b/src/common/aka_iterators.hh @@ -1,354 +1,459 @@ /** * @file aka_iterators.hh * * @author Nicolas Richart * * @date creation: Fri Aug 11 2017 * @date last modification: Mon Jan 29 2018 * * @brief iterator interfaces * * @section LICENSE * * 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 "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_ITERATORS_HH__ #define __AKANTU_AKA_ITERATORS_HH__ namespace akantu { namespace tuple { /* ------------------------------------------------------------------------ */ namespace details { template struct Foreach { template static inline bool not_equal(Tuple && a, Tuple && b) { if (std::get(std::forward(a)) == std::get(std::forward(b))) return false; return Foreach::not_equal(std::forward(a), std::forward(b)); } }; /* ---------------------------------------------------------------------- */ template <> struct Foreach<0> { template static inline bool not_equal(Tuple && a, Tuple && b) { return std::get<0>(std::forward(a)) != std::get<0>(std::forward(b)); } }; template decltype(auto) make_tuple_no_decay(Ts &&... args) { return std::tuple(std::forward(args)...); } template void foreach_impl(F && func, Tuple && tuple, std::index_sequence &&) { (void)std::initializer_list{ (std::forward(func)(std::get(std::forward(tuple))), 0)...}; } template decltype(auto) transform_impl(F && func, Tuple && tuple, std::index_sequence &&) { return make_tuple_no_decay( std::forward(func)(std::get(std::forward(tuple)))...); } } // namespace details /* ------------------------------------------------------------------------ */ template bool are_not_equal(Tuple && a, Tuple && b) { return details::Foreach>::value>:: not_equal(std::forward(a), std::forward(b)); } template void foreach (F && func, Tuple && tuple) { return details::foreach_impl( std::forward(func), std::forward(tuple), std::make_index_sequence< std::tuple_size>::value>{}); } template decltype(auto) transform(F && func, Tuple && tuple) { return details::transform_impl( std::forward(func), std::forward(tuple), std::make_index_sequence< std::tuple_size>::value>{}); } } // namespace tuple /* -------------------------------------------------------------------------- */ namespace iterators { namespace { template struct is_at_least_category : std::is_same::iterator_category, category>, category> {}; - } + } // namespace template class ZipIterator { public: using value_type = std::tuple::value_type...>; using difference_type = std::common_type_t< typename std::iterator_traits::difference_type...>; using pointer = std::tuple::pointer...>; using reference = std::tuple::reference...>; using iterator_category = std::input_iterator_tag; // std::common_type_t::iterator_category...>; private: using tuple_t = std::tuple; public: explicit ZipIterator(tuple_t iterators) : iterators(std::move(iterators)) {} // input iterator ++it ZipIterator & operator++() { tuple::foreach ([](auto && it) { ++it; }, iterators); return *this; } // input iterator it++ ZipIterator operator++(int) { auto cpy = *this; tuple::foreach ([](auto && it) { ++it; }, iterators); return cpy; } // input iterator it != other_it bool operator!=(const ZipIterator & other) const { return tuple::are_not_equal(iterators, other.iterators); } // input iterator dereference *it decltype(auto) operator*() { return tuple::transform([](auto && it) -> decltype(auto) { return *it; }, iterators); } bool operator==(const ZipIterator & other) const { return not tuple::are_not_equal(iterators, other.iterators); } private: tuple_t iterators; }; } // namespace iterators /* -------------------------------------------------------------------------- */ template decltype(auto) zip_iterator(std::tuple && iterators_tuple) { auto zip = iterators::ZipIterator( std::forward(iterators_tuple)); return zip; } /* -------------------------------------------------------------------------- */ namespace containers { template class ZipContainer { using containers_t = std::tuple; public: explicit ZipContainer(Containers &&... containers) : containers(std::forward(containers)...) {} decltype(auto) begin() const { return zip_iterator( tuple::transform([](auto && c) { return c.begin(); }, std::forward(containers))); } decltype(auto) end() const { return zip_iterator( tuple::transform([](auto && c) { return c.end(); }, std::forward(containers))); } decltype(auto) begin() { return zip_iterator( tuple::transform([](auto && c) { return c.begin(); }, std::forward(containers))); } decltype(auto) end() { return zip_iterator( tuple::transform([](auto && c) { return c.end(); }, std::forward(containers))); } private: containers_t containers; }; template class Range { public: explicit Range(Iterator && it1, Iterator && it2) : iterators(std::forward(it1), std::forward(it2)) {} decltype(auto) begin() const { return std::get<0>(iterators); } decltype(auto) begin() { return std::get<0>(iterators); } decltype(auto) end() const { return std::get<1>(iterators); } decltype(auto) end() { return std::get<1>(iterators); } private: std::tuple iterators; }; } // namespace containers /* -------------------------------------------------------------------------- */ template decltype(auto) zip(Containers &&... conts) { return containers::ZipContainer( std::forward(conts)...); } template decltype(auto) range(Iterator && it1, Iterator && it2) { return containers::Range(std::forward(it1), std::forward(it2)); } /* -------------------------------------------------------------------------- */ /* Arange */ /* -------------------------------------------------------------------------- */ namespace iterators { template class ArangeIterator { public: using value_type = T; using pointer = T *; using reference = T &; using difference_type = size_t; using iterator_category = std::input_iterator_tag; constexpr ArangeIterator(T value, T step) : value(value), step(step) {} constexpr ArangeIterator(const ArangeIterator &) = default; constexpr ArangeIterator & operator++() { value += step; return *this; } constexpr T operator*() const { return value; } constexpr bool operator==(const ArangeIterator & other) const { return (value == other.value) and (step == other.step); } constexpr bool operator!=(const ArangeIterator & other) const { return not operator==(other); } private: T value{0}; const T step{1}; }; } // namespace iterators namespace containers { template class ArangeContainer { public: using iterator = iterators::ArangeIterator; + using const_iterator = iterators::ArangeIterator; constexpr ArangeContainer(T start, T stop, T step = 1) : start(start), stop((stop - start) % step == 0 ? stop : start + (1 + (stop - start) / step) * step), step(step) {} explicit constexpr ArangeContainer(T stop) : ArangeContainer(0, stop, 1) {} constexpr T operator[](size_t i) { T val = start + i * step; assert(val < stop && "i is out of range"); return val; } constexpr T size() { return (stop - start) / step; } constexpr iterator begin() { return iterator(start, step); } constexpr iterator end() { return iterator(stop, step); } private: const T start{0}, stop{0}, step{1}; }; } // namespace containers template >::value>> inline decltype(auto) arange(const T & stop) { return containers::ArangeContainer(stop); } template >::value>> inline constexpr decltype(auto) arange(const T1 & start, const T2 & stop) { return containers::ArangeContainer>(start, stop); } template >::value>> inline constexpr decltype(auto) arange(const T1 & start, const T2 & stop, const T3 & step) { return containers::ArangeContainer>( start, stop, step); } /* -------------------------------------------------------------------------- */ template inline constexpr decltype(auto) enumerate(Container && container, size_t start_ = 0) { auto stop = std::forward(container).size(); decltype(stop) start = start_; return zip(arange(start, stop), std::forward(container)); } +namespace iterators { + template + class transform_adaptor_iterator { + public: + using value_type = decltype(std::declval()( + std::declval())); + using difference_type = typename iterator_t::difference_type; + using pointer = std::decay_t *; + using reference = value_type &; + using iterator_category = typename iterator_t::iterator_category; + + transform_adaptor_iterator(iterator_t it, operator_t && op) + : it(std::move(it)), op(op) {} + transform_adaptor_iterator(const transform_adaptor_iterator &) = default; + + transform_adaptor_iterator & operator++() { + ++it; + return *this; + } + + decltype(auto) operator*() const { return op(*it); } + + bool operator==(const transform_adaptor_iterator & other) const { + return (it == other.it); + } + + bool operator!=(const transform_adaptor_iterator & other) const { + return not operator==(other); + } + + private: + iterator_t it; + operator_t op; + }; + + template + decltype(auto) make_transform_adaptor_iterator(iterator_t it, + operator_t && op) { + return transform_adaptor_iterator( + it, std::forward(op)); + } + +} // namespace iterators + +namespace containers { + template + class TransformIteratorAdaptor { + public: + using const_iterator = typename std::decay_t::const_iterator; + using iterator = typename std::decay_t::iterator; + + TransformIteratorAdaptor(container_t && cont, operator_t && op) + : cont(std::forward(cont)), + op(std::forward(op)) {} + + decltype(auto) begin() const { + return iterators::make_transform_adaptor_iterator(cont.begin(), op); + } + decltype(auto) begin() { + return iterators::make_transform_adaptor_iterator(cont.begin(), op); + } + + decltype(auto) end() const { + return iterators::make_transform_adaptor_iterator(cont.end(), op); + } + decltype(auto) end() { + return iterators::make_transform_adaptor_iterator(cont.end(), op); + } + + private: + container_t cont; + operator_t op; + }; +} // namespace containers + +template +decltype(auto) make_transform_adaptor(container_t && cont, operator_t && op) { + return containers::TransformIteratorAdaptor( + std::forward(cont), std::forward(op)); +} + +template +decltype(auto) make_keys_adaptor(container_t && cont) { + return make_transform_adaptor( + std::forward(cont), + [](auto && pair) -> decltype(pair.first) { return pair.first; }); +} + +template +decltype(auto) make_values_adaptor(container_t && cont) { + return make_transform_adaptor( + std::forward(cont), + [](auto && pair) -> decltype(pair.second) { return pair.second; }); +} + +template +decltype(auto) make_dereference_adaptor(container_t && cont) { + return make_transform_adaptor( + std::forward(cont), + [](auto && value) -> decltype(*value) { return *value; }); +} + } // namespace akantu namespace std { template struct iterator_traits<::akantu::iterators::ZipIterator> { using iterator_category = forward_iterator_tag; using value_type = typename ::akantu::iterators::ZipIterator::value_type; using difference_type = typename ::akantu::iterators::ZipIterator::difference_type; using pointer = typename ::akantu::iterators::ZipIterator::pointer; using reference = typename ::akantu::iterators::ZipIterator::reference; }; -} +} // namespace std #endif /* __AKANTU_AKA_ITERATORS_HH__ */ diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index 242dd2ef4..c5a93e998 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,316 +1,324 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Feb 07 2018 * * @brief Stores information relevent to the notion of element and nodes * groups. * * @section LICENSE * * 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_GROUP_MANAGER_HH__ #define __AKANTU_GROUP_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" +#include "aka_iterators.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementGroup; class NodeGroup; class Mesh; class Element; class ElementSynchronizer; template class CommunicationBufferTemplated; namespace dumper { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: #ifdef SWIGPYTHON public: using ElementGroups = std::map; using NodeGroups = std::map; private: #else using ElementGroups = std::map; using NodeGroups = std::map; #endif public: using GroupManagerTypeSet = std::set; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GroupManager(const Mesh & mesh, const ID & id = "group_manager", const MemoryID & memory_id = 0); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ public: using node_group_iterator = NodeGroups::iterator; using element_group_iterator = ElementGroups::iterator; using const_node_group_iterator = NodeGroups::const_iterator; using const_element_group_iterator = ElementGroups::const_iterator; #ifndef SWIG #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, \ param_in, param_out) \ inline BOOST_PP_CAT(BOOST_PP_CAT(const_, group_type), _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ }; \ \ inline BOOST_PP_CAT(group_type, _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ } #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \ AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION( \ group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY()) AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name); #endif +public: + decltype(auto) iterateNodeGroups() { + return make_dereference_adaptor(make_values_adaptor(node_groups)); + } + decltype(auto) iterateNodeGroups() const { + return make_dereference_adaptor(make_values_adaptor(node_groups)); + } /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* ------------------------------------------------------------------------ */ public: class ClusteringFilter { public: virtual bool operator()(const Element &) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create a node group from another node group but filtered template NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// destroy a node group void destroyNodeGroup(const std::string & group_name); /// create an element group and the associated node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension = _all_dimensions, bool replace_group = false); /// create an element group from another element group but filtered template ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); /// destroy all element groups and the associated node groups void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// create groups based on values stored in a given mesh data template void createGroupsFromMeshData(const std::string & dataset_name); /// create boundaries group by a clustering algorithm \todo extend to parallel UInt createBoundaryGroupFromGeometry(); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); /// register an elemental field to the given group name (overloading for /// ElementalPartionField) #ifndef SWIG template class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// ElementalField) template class ret_type, template class, bool> class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template class dump_type> dumper::Field * createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); template class ftype> dumper::Field * createNodalField(const ftype * field, const std::string & group_name, UInt padding_size = 0); template class ftype> dumper::Field * createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size); protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated & comm_buffer) const; /// take a buffer and create the missing groups localy void checkAndAddGroups(CommunicationBufferTemplated & buffer); /// register an elemental field to the given group name template inline dumper::Field * createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem); /// register an elemental field to the given group name template inline dumper::Field * createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); #endif /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; UInt getNbNodeGroups() { return node_groups.size(); }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id to create element and node groups ID id; /// memory_id to create element and node groups MemoryID memory_id; /// list of the node groups managed NodeGroups node_groups; /// list of the element groups managed ElementGroups element_groups; /// Mesh to which the element belongs const Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const GroupManager & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AKANTU_GROUP_MANAGER_HH__ */ diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index 5c16acdb5..7f77b46de 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,585 +1,586 @@ /** * @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 * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "group_manager_inline_impl.cc" #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" /* -------------------------------------------------------------------------- */ #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), lower_bounds(spatial_dimension, 0.), upper_bounds(spatial_dimension, 0.), size(spatial_dimension, 0.), local_lower_bounds(spatial_dimension, 0.), local_upper_bounds(spatial_dimension, 0.), - mesh_data("mesh_data", id, memory_id), communicator(&communicator) { + 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"); 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)); } } /* -------------------------------------------------------------------------- */ Mesh & Mesh::initMeshFacets(const ID & id) { AKANTU_DEBUG_IN(); if (!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_type = this->nodes_type; 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()); } /// transfers the the mesh physical names to the mesh facets if (not this->hasData("physical_names")) { AKANTU_DEBUG_OUT(); return *mesh_facets; } if (not mesh_facets->hasData("physical_names")) { - mesh_facets->registerData("physical_names"); + mesh_facets->registerElementalData("physical_names"); } 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) { MeshIO mesh_io; mesh_io.read(filename, *this, mesh_io_type); type_iterator it = this->firstType(spatial_dimension, _not_ghost, _ek_not_defined); type_iterator last = this->lastType(spatial_dimension, _not_ghost, _ek_not_defined); 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 << " where read."); this->nb_global_nodes = this->nodes->size(); this->computeBoundingBox(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ 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::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += 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(); for (UInt k = 0; k < spatial_dimension; ++k) { local_lower_bounds(k) = std::numeric_limits::max(); local_upper_bounds(k) = -std::numeric_limits::max(); } for (UInt i = 0; i < nodes->size(); ++i) { // if(!isPureGhostNode(i)) for (UInt k = 0; k < spatial_dimension; ++k) { local_lower_bounds(k) = std::min(local_lower_bounds[k], (*nodes)(i, k)); local_upper_bounds(k) = std::max(local_upper_bounds[k], (*nodes)(i, k)); } } if (this->is_distributed) { Matrix reduce_bounds(spatial_dimension, 2); for (UInt k = 0; k < spatial_dimension; ++k) { reduce_bounds(k, 0) = local_lower_bounds(k); reduce_bounds(k, 1) = -local_upper_bounds(k); } communicator->allReduce(reduce_bounds, SynchronizerOperation::_min); for (UInt k = 0; k < spatial_dimension; ++k) { lower_bounds(k) = reduce_bounds(k, 0); upper_bounds(k) = -reduce_bounds(k, 1); } } else { this->lower_bounds = this->local_lower_bounds; this->upper_bounds = this->local_upper_bounds; } size = upper_bounds - lower_bounds; 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, const ElementKind & element_kind) { ElementTypeMap nb_data_per_elem; for (auto type : elementTypes(spatial_dimension, _not_ghost, element_kind)) { 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, const ElementKind & element_kind); template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind) { dumper::Field * field = nullptr; ElementTypeMapArray * internal = nullptr; try { internal = &(this->getData(field_id)); } catch (...) { return nullptr; } ElementTypeMap nb_data_per_elem = this->getNbDataPerElem(*internal, element_kind); field = this->createElementalField( *internal, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); return field; } template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); template dumper::Field * Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); #endif /* -------------------------------------------------------------------------- */ void Mesh::distribute() { this->distribute(Communicator::getStaticCommunicator()); } /* -------------------------------------------------------------------------- */ void Mesh::distribute(Communicator & communicator) { 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) { // return; // } #ifdef AKANTU_USE_SCOTCH Int prank = this->communicator->whoAmI(); if (prank == 0) { MeshPartitionScotch partition(*this, spatial_dimension); partition.partitionate(psize); 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 this->computeBoundingBox(); } this->is_distributed = true; } /* -------------------------------------------------------------------------- */ 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); 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) 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); 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; } this->sendEvent(event); } } // namespace akantu diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh index a9e0b4a2a..a3418af96 100644 --- a/src/mesh/mesh.hh +++ b/src/mesh/mesh.hh @@ -1,622 +1,610 @@ /** * @file mesh.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Feb 19 2018 * * @brief the class representing the meshes * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_event_handler_manager.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element.hh" #include "element_class.hh" #include "element_type_map.hh" #include "group_manager.hh" #include "mesh_data.hh" #include "mesh_events.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class Communicator; class ElementSynchronizer; class NodeSynchronizer; class MeshGlobalDataUpdater; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Array, and the connectivity. The connectivity are stored in by element * types. * * In order to loop on all element you have to loop on all types like this : * @code{.cpp} for(auto & type : mesh.elementTypes()) { UInt nb_element = mesh.getNbElement(type); const Array & conn = mesh.getConnectivity(type); for(UInt e = 0; e < nb_element; ++e) { ... } } or for_each_element(mesh, [](Element & element) { std::cout << element << std::endl }); @endcode */ class Mesh : protected Memory, public EventHandlerManager, public GroupManager, + public MeshData, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: /// default constructor used for chaining, the last parameter is just to /// differentiate constructors Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator); public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// mesh not distributed and not using the default communicator Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id = "mesh", const MemoryID & memory_id = 0); /// constructor that use an existing nodes coordinates array, by knowing its /// ID // Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id, // const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); ~Mesh() override; /// read the mesh from a file void read(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); /// write the mesh to a file void write(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); private: /// initialize the connectivity to NULL and other stuff void init(); /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// patitionate the mesh among the processors involved in their computation virtual void distribute(Communicator & communicator); virtual void distribute(); /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Array & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; // /// extract coordinates of nodes from a reversed element // inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, // UInt * connectivity, // UInt n_nodes); /// add a Array of connectivity for the type . inline void addConnectivityType(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ template inline void sendEvent(Event & event) { // if(event.getList().size() != 0) EventHandlerManager::sendEvent(event); } /// prepare the event to remove the elements listed void eraseElements(const Array & elements); /* ------------------------------------------------------------------------ */ template inline void removeNodesFromArray(Array & vect, const Array & new_numbering); /// initialize normals void initNormals(); /// init facets' mesh Mesh & initMeshFacets(const ID & id = "mesh_facets"); /// define parent mesh void defineMeshParent(const Mesh & mesh); /// get global connectivity array void getGlobalConnectivity(ElementTypeMapArray & global_connectivity); public: void getAssociatedElements(const Array & node_list, Array & elements); private: /// fills the nodes_to_elements structure void fillNodesToElements(); /// update the global ids, nodes type, ... std::tuple updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event); void registerGlobalDataUpdater( std::unique_ptr && global_data_updater); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the id of the mesh AKANTU_GET_MACRO(ID, Memory::id, const ID &); /// get the id of the mesh AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &); /// get the spatial dimension of the mesh = number of component of the /// coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Array aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Array &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array &); /// get the normals for the elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt); /// get the Array of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array &); AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global id of a node inline UInt getNodeLocalId(UInt global_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Array AKANTU_GET_MACRO(NodesType, *nodes_type, const Array &); protected: AKANTU_GET_MACRO_NOT_CONST(NodesType, *nodes_type, Array &); public: inline NodeType getNodeType(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; AKANTU_GET_MACRO(LowerBounds, lower_bounds, const Vector &); AKANTU_GET_MACRO(UpperBounds, upper_bounds, const Vector &); AKANTU_GET_MACRO(LocalLowerBounds, local_lower_bounds, const Vector &); AKANTU_GET_MACRO(LocalUpperBounds, local_upper_bounds, const Vector &); /// get the connectivity Array for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); AKANTU_GET_MACRO(Connectivities, connectivities, const ElementTypeMapArray &); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get the number of element for a given ghost_type and a given dimension inline UInt getNbElement(const UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_not_defined) const; /// compute the barycenter of a given element inline void getBarycenter(const Element & element, Vector & barycenter) const; void getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const; #ifndef SWIG /// get the element connected to a subelement (element of lower dimension) const auto & getElementToSubelement() const; /// get the element connected to a subelement const auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the element connected to a subelement auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the elements connected to a subelement const auto & getElementToSubelement(const Element & element) const; /// get the subelement (element of lower dimension) connected to a element const auto & getSubelementToElement() const; /// get the subelement connected to an element const auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the subelement connected to an element auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the subelement (element of lower dimension) connected to a element VectorProxy getSubelementToElement(const Element & element) const; /// get connectivity of a given element inline VectorProxy getConnectivity(const Element & element) const; protected: inline auto & getElementToSubelement(const Element & element); inline VectorProxy getSubelementToElement(const Element & element); inline VectorProxy getConnectivity(const Element & element); #endif public: - /// register a new ElementalTypeMap in the MeshData - template - inline ElementTypeMapArray & registerData(const std::string & data_name); - - template - inline bool hasData(const ID & data_name, const ElementType & el_type, - const GhostType & ghost_type = _not_ghost) const; - /// get a name field associated to the mesh template inline const Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get a name field associated to the mesh template inline Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost); - /// tells if the data data_name is contained in the mesh or not - inline bool hasData(const ID & data_name) const; - /// get a name field associated to the mesh template inline const ElementTypeMapArray & getData(const ID & data_name) const; /// get a name field associated to the mesh template inline ElementTypeMapArray & getData(const ID & data_name); template ElementTypeMap getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); template dumper::Field * createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); /// templated getter returning the pointer to data in MeshData (modifiable) template inline Array & getDataPointer(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost, UInt nb_component = 1, bool size_to_nb_element = true, bool resize_with_parent = false); template inline Array & getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_); /// Facets mesh accessor inline const Mesh & getMeshFacets() const; inline Mesh & getMeshFacets(); /// Parent mesh accessor inline const Mesh & getMeshParent() const; inline bool isMeshFacets() const { return this->is_mesh_facets; } /// defines is the mesh is distributed or not inline bool isDistributed() const { return this->is_distributed; } #ifndef SWIG /// return the dumper from a group and and a dumper name DumperIOHelper & getGroupDumper(const std::string & dumper_name, const std::string & group_name); #endif /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type, UInt t); #ifndef SWIG /// get local connectivity of a facet for a given facet type static inline auto getFacetLocalConnectivity(const ElementType & type, UInt t = 0); /// get connectivity of facets for a given element inline auto getFacetConnectivity(const Element & element, UInt t = 0) const; #endif /// get the number of type of the surface element associated to a given /// element type static inline UInt getNbFacetTypes(const ElementType & type, UInt t = 0); #ifndef SWIG /// get the type of the surface element associated to a given element static inline constexpr auto getFacetType(const ElementType & type, UInt t = 0); /// get all the type of the surface element associated to a given element static inline constexpr auto getAllFacetTypes(const ElementType & type); #endif /// get the number of nodes in the given element list static inline UInt getNbNodesPerElementList(const Array & elements); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ using type_iterator = ElementTypeMapArray::type_iterator; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; template ElementTypesIteratorHelper elementTypes(pack &&... _pack) const; inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.firstType(dim, ghost_type, kind); } inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.lastType(dim, ghost_type, kind); } AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer, const ElementSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer, ElementSynchronizer &); AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer, const NodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer, NodeSynchronizer &); // AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator // &); #ifndef SWIG AKANTU_GET_MACRO(Communicator, *communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &); #endif /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshAccessor; friend class MeshUtils; AKANTU_GET_MACRO(NodesPointer, *nodes, Array &); /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline Array & getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Array and create it if necessary inline Array & getNodesTypePointer(); /// get a pointer to the connectivity Array for the given type and create it /// if necessary inline Array & getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the ghost element counter inline Array & getGhostsCounters(const ElementType & type, const GhostType & ghost_type = _ghost) { AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost, "No ghost counter for _not_ghost elements"); return ghosts_counters(type, ghost_type); } /// get a pointer to the element_to_subelement Array for the given type and /// create it if necessary inline Array> & getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Array for the given type and /// create it if necessary inline Array & getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); - AKANTU_GET_MACRO_NOT_CONST(MeshData, mesh_data, MeshData &); - /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of the nodes coordinates std::shared_ptr> nodes; /// global node ids std::shared_ptr> nodes_global_ids; /// node type, -3 pure ghost, -2 master for the node, -1 normal node, i in /// [0-N] slave node and master is proc i std::shared_ptr> nodes_type; /// global number of nodes; UInt nb_global_nodes{0}; /// all class of elements present in this mesh (for heterogenous meshes) ElementTypeMapArray connectivities; /// count the references on ghost elements ElementTypeMapArray ghosts_counters; /// map to normals for all class of elements present in this mesh ElementTypeMapArray normals; /// the spatial dimension of this mesh UInt spatial_dimension{0}; /// min of coordinates Vector lower_bounds; /// max of coordinates Vector upper_bounds; /// size covered by the mesh on each direction Vector size; /// local min of coordinates Vector local_lower_bounds; /// local max of coordinates Vector local_upper_bounds; /// Extra data loaded from the mesh file - MeshData mesh_data; + //MeshData mesh_data; /// facets' mesh std::unique_ptr mesh_facets; /// parent mesh (this is set for mesh_facets meshes) const Mesh * mesh_parent{nullptr}; /// defines if current mesh is mesh_facets or not bool is_mesh_facets{false}; /// defines if the mesh is centralized or distributed bool is_distributed{false}; /// Communicator on which mesh is distributed Communicator * communicator; /// Element synchronizer std::unique_ptr element_synchronizer; /// Node synchronizer std::unique_ptr node_synchronizer; using NodesToElements = std::vector>>; /// class to update global data using external knowledge std::unique_ptr global_data_updater; /// This info is stored to simplify the dynamic changes NodesToElements nodes_to_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ #include "element_type_map_tmpl.hh" #include "mesh_inline_impl.cc" #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/mesh/mesh_accessor.hh b/src/mesh/mesh_accessor.hh index c841b5f1b..fb2364af3 100644 --- a/src/mesh/mesh_accessor.hh +++ b/src/mesh/mesh_accessor.hh @@ -1,150 +1,148 @@ /** * @file mesh_accessor.hh * * @author Nicolas Richart * * @date creation: Tue Jun 30 2015 * @date last modification: Tue Sep 19 2017 * * @brief this class allow to access some private member of mesh it is used for * IO for examples * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_ACCESSOR_HH__ #define __AKANTU_MESH_ACCESSOR_HH__ namespace akantu { class NodeSynchronizer; class ElementSynchronizer; class MeshGlobalDataUpdater; } // namespace akantu namespace akantu { class MeshAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: explicit MeshAccessor(Mesh & mesh) : _mesh(mesh) {} virtual ~MeshAccessor() = default; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the global number of nodes inline UInt getNbGlobalNodes() const { return this->_mesh.nb_global_nodes; } /// set the global number of nodes inline void setNbGlobalNodes(UInt nb_global_nodes) { this->_mesh.nb_global_nodes = nb_global_nodes; } /// set the mesh as being distributed inline void setDistributed() { this->_mesh.is_distributed = true; } /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline auto & getNodesGlobalIds() { return this->_mesh.getNodesGlobalIdsPointer(); } /// get a pointer to the nodes_type Array and create it if necessary inline auto & getNodesType() { return this->_mesh.getNodesTypePointer(); } /// get a pointer to the coordinates Array inline auto & getNodes() { return this->_mesh.getNodesPointer(); } /// get a pointer to the coordinates Array inline auto getNodesSharedPtr() { return this->_mesh.nodes; } /// get a pointer to the connectivity Array for the given type and create it /// if necessary inline auto & getConnectivity(const ElementType & type, const GhostType & ghost_type = _not_ghost) { return this->_mesh.getConnectivityPointer(type, ghost_type); } /// get the ghost element counter inline auto & getGhostsCounters(const ElementType & type, const GhostType & ghost_type = _ghost) { return this->_mesh.getGhostsCounters(type, ghost_type); } /// get a pointer to the element_to_subelement Array for the given type and /// create it if necessary inline auto & getElementToSubelement(const ElementType & type, const GhostType & ghost_type = _not_ghost) { return this->_mesh.getElementToSubelementPointer(type, ghost_type); } /// get a pointer to the subelement_to_element Array for the given type and /// create it if necessary inline auto & getSubelementToElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) { return this->_mesh.getSubelementToElementPointer(type, ghost_type); } template inline auto & getData(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost, UInt nb_component = 1, bool size_to_nb_element = true, bool resize_with_parent = false) { return this->_mesh.getDataPointer(data_name, el_type, ghost_type, nb_component, size_to_nb_element, resize_with_parent); } - auto & getMeshData() { return this->_mesh.getMeshData(); } - /// get the node synchonizer auto & getNodeSynchronizer() { return *this->_mesh.node_synchronizer; } /// get the element synchonizer auto & getElementSynchronizer() { return *this->_mesh.element_synchronizer; } decltype(auto) updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { return this->_mesh.updateGlobalData(nodes_event, elements_event); } void registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->_mesh.registerGlobalDataUpdater( std::forward>( global_data_updater)); } private: Mesh & _mesh; }; } // namespace akantu #endif /* __AKANTU_MESH_ACCESSOR_HH__ */ diff --git a/src/mesh/mesh_data.cc b/src/mesh/mesh_data.cc index ead9bcf7c..31a9c97b2 100644 --- a/src/mesh/mesh_data.cc +++ b/src/mesh/mesh_data.cc @@ -1,47 +1,40 @@ /** * @file mesh_data.cc * * @author Dana Christen * * @date creation: Fri Apr 13 2012 * @date last modification: Mon Jun 19 2017 * * @brief Stores generic data loaded from the mesh file * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #include "mesh_data.hh" #include "mesh.hh" namespace akantu { MeshData::MeshData(const ID & _id, const ID & parent_id, const MemoryID & mem_id) - : Memory(parent_id + ":" + _id, mem_id) {} + : _id(parent_id + ":" + _id), _memory_id(mem_id) {} -MeshData::~MeshData() { - ElementalDataMap::iterator it, end; - for (it = elemental_data.begin(); it != elemental_data.end(); ++it) { - delete it->second; - } -} - -} // akantu +} // namespace akantu diff --git a/src/mesh/mesh_data.hh b/src/mesh/mesh_data.hh index 45650efae..fb5b4faf1 100644 --- a/src/mesh/mesh_data.hh +++ b/src/mesh/mesh_data.hh @@ -1,153 +1,189 @@ /** * @file mesh_data.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Mon Dec 18 2017 * * @brief Stores generic data loaded from the mesh file * * @section LICENSE * * 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_MESH_DATA_HH__ #define __AKANTU_MESH_DATA_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "element_type_map.hh" #include #include /* -------------------------------------------------------------------------- */ namespace akantu { #define AKANTU_MESH_DATA_TYPES \ ((_tc_int, Int))((_tc_uint, UInt))((_tc_real, Real))( \ (_tc_element, Element))((_tc_std_string, std::string))( \ (_tc_std_vector_element, std::vector)) #define AKANTU_MESH_DATA_TUPLE_FIRST_ELEM(s, data, elem) \ BOOST_PP_TUPLE_ELEM(2, 0, elem) enum MeshDataTypeCode : int { BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_MESH_DATA_TUPLE_FIRST_ELEM, , AKANTU_MESH_DATA_TYPES)), _tc_unknown }; -class MeshData : public Memory { +enum class MeshDataType { + _nodal, + _elemental, +}; + +class MeshData { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: using TypeCode = MeshDataTypeCode; - using ElementalDataMap = std::map; - using StringVector = std::vector; + using ElementalDataMap = + std::map>; + using NodalDataMap = std::map>; using TypeCodeMap = std::map; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshData(const ID & id = "mesh_data", const ID & parent_id = "", const MemoryID & memory_id = 0); - ~MeshData() override; /* ------------------------------------------------------------------------ */ /* Methods and accessors */ /* ------------------------------------------------------------------------ */ public: /// Register new elemental data (and alloc data) with check if the name is /// new - template void registerElementalData(const ID & name); + template ElementTypeMapArray & registerElementalData(const ID & name); inline void registerElementalData(const ID & name, TypeCode type); + /// Register new nodal data (and alloc data) with check if the name is + /// new + template + Array & registerNodalData(const ID & name, UInt nb_components = 1); + inline void registerNodalData(const ID & name, UInt nb_components, + TypeCode type); + /// tells if the given array exists template - bool hasDataArray(const ID & data_name, const ElementType & el_type, - const GhostType & ghost_type = _not_ghost) const; + bool hasData(const ID & data_name, const ElementType & el_type, + const GhostType & ghost_type = _not_ghost) const; /// tells if the given data exists - bool hasData(const ID & data_name) const; + bool hasData(const ID & data_name, + MeshDataType type = MeshDataType::_elemental) const; /// Get an existing elemental data array template const Array & getElementalDataArray(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; template Array & getElementalDataArray(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// Get an elemental data array, if it does not exist: allocate it template Array & getElementalDataArrayAlloc(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost, UInt nb_component = 1); /// get the names of the data stored in elemental_data - inline void getTagNames(StringVector & tags, const ElementType & type, + inline auto getTagNames(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; + /// get the names of the data stored in elemental_data + inline auto getTagNames() const; + /// get the type of the data stored in elemental_data template TypeCode getTypeCode() const; - inline TypeCode getTypeCode(const ID & name) const; + inline TypeCode + getTypeCode(const ID & name, + MeshDataType type = MeshDataType::_elemental) const; template inline UInt getNbComponentTemplated(const ID & name, const ElementType & el_type, const GhostType & ghost_type) const; inline UInt getNbComponent(const ID & name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; + inline UInt getNbComponent(const ID & name) const; + /// Get an existing elemental data template const ElementTypeMapArray & getElementalData(const ID & name) const; template ElementTypeMapArray & getElementalData(const ID & name); + template + Array & getNodalData(const ID & name, UInt nb_components = 1); + template const Array & getNodalData(const ID & name) const; + private: /// Register new elemental data (add alloc data) template - ElementTypeMapArray * allocElementalData(const ID & name); + ElementTypeMapArray & allocElementalData(const ID & name); + + /// Register new nodal data (add alloc data) + template + Array & allocNodalData(const ID & name, UInt nb_components); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: + ID _id; + UInt _memory_id{0}; + /// Map when elemental data is stored as ElementTypeMap ElementalDataMap elemental_data; + + /// Map when elemental data is stored as ElementTypeMap + NodalDataMap nodal_data; + /// Map when elementalType of the data stored in elemental_data - TypeCodeMap typecode_map; + std::map typecode_map{ + {MeshDataType::_elemental, {}}, {MeshDataType::_nodal, {}}}; }; -} // akantu +} // namespace akantu #include "mesh_data_tmpl.hh" #undef AKANTU_MESH_DATA_TUPLE_FIRST_ELEM #endif /* __AKANTU_MESH_DATA_HH__ */ diff --git a/src/mesh/mesh_data_tmpl.hh b/src/mesh/mesh_data_tmpl.hh index 7de088649..a627a5d90 100644 --- a/src/mesh/mesh_data_tmpl.hh +++ b/src/mesh/mesh_data_tmpl.hh @@ -1,276 +1,378 @@ /** * @file mesh_data_tmpl.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Tue Feb 20 2018 * * @brief Stores generic data loaded from the mesh file * * @section LICENSE * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "mesh_data.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_DATA_TMPL_HH__ #define __AKANTU_MESH_DATA_TMPL_HH__ namespace akantu { #define MESH_DATA_GET_TYPE(r, data, type) \ template <> \ inline MeshDataTypeCode \ MeshData::getTypeCode() const { \ return BOOST_PP_TUPLE_ELEM(2, 0, type); \ } /* -------------------------------------------------------------------------- */ // get the type of the data stored in elemental_data template inline MeshDataTypeCode MeshData::getTypeCode() const { AKANTU_ERROR("Type " << debug::demangle(typeid(T).name()) << "not implemented by MeshData."); } /* -------------------------------------------------------------------------- */ BOOST_PP_SEQ_FOR_EACH(MESH_DATA_GET_TYPE, void, AKANTU_MESH_DATA_TYPES) #undef MESH_DATA_GET_TYPE -inline MeshDataTypeCode MeshData::getTypeCode(const ID & name) const { - auto it = typecode_map.find(name); - if (it == typecode_map.end()) +inline MeshDataTypeCode MeshData::getTypeCode(const ID & name, + MeshDataType type) const { + auto it = typecode_map.at(type).find(name); + if (it == typecode_map.at(type).end()) AKANTU_EXCEPTION("No dataset named " << name << " found."); return it->second; } /* -------------------------------------------------------------------------- */ // Register new elemental data templated (and alloc data) with check if the // name is new -template void MeshData::registerElementalData(const ID & name) { +template +ElementTypeMapArray & MeshData::registerElementalData(const ID & name) { auto it = elemental_data.find(name); if (it == elemental_data.end()) { - allocElementalData(name); + return allocElementalData(name); } else { AKANTU_DEBUG_INFO("Data named " << name << " already registered."); + return getElementalData(name); } } /* -------------------------------------------------------------------------- */ // Register new elemental data of a given MeshDataTypeCode with check if the // name is new #define AKANTU_MESH_DATA_CASE_MACRO(r, name, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ registerElementalData(name); \ break; \ } inline void MeshData::registerElementalData(const ID & name, MeshDataTypeCode type) { switch (type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_MESH_DATA_CASE_MACRO, name, AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Type " << type << "not implemented by MeshData."); } } #undef AKANTU_MESH_DATA_CASE_MACRO /* -------------------------------------------------------------------------- */ /// Register new elemental data (and alloc data) template -ElementTypeMapArray * MeshData::allocElementalData(const ID & name) { - auto * dataset = new ElementTypeMapArray(name, id, memory_id); - elemental_data[name] = dataset; - typecode_map[name] = getTypeCode(); - return dataset; +ElementTypeMapArray & MeshData::allocElementalData(const ID & name) { + auto dataset = + std::make_unique>(name, _id, _memory_id); + auto * dataset_typed = dataset.get(); + elemental_data[name] = std::move(dataset); + typecode_map[MeshDataType::_elemental][name] = getTypeCode(); + return *dataset_typed; +} + +/* -------------------------------------------------------------------------- */ +// Register new nodal data templated (and alloc data) with check if the +// name is new +template +Array & MeshData::registerNodalData(const ID & name, UInt nb_components) { + auto it = nodal_data.find(name); + if (it == nodal_data.end()) { + return allocNodalData(name, nb_components); + } else { + AKANTU_DEBUG_INFO("Data named " << name << " already registered."); + return getNodalData(name); + } +} + +/* -------------------------------------------------------------------------- */ +// Register new elemental data of a given MeshDataTypeCode with check if the +// name is new +#define AKANTU_MESH_NODAL_DATA_CASE_MACRO(r, name, elem) \ + case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ + registerNodalData(name, nb_components); \ + break; \ + } + +inline void MeshData::registerNodalData(const ID & name, UInt nb_components, + MeshDataTypeCode type) { + switch (type) { + BOOST_PP_SEQ_FOR_EACH(AKANTU_MESH_NODAL_DATA_CASE_MACRO, name, + AKANTU_MESH_DATA_TYPES) + default: + AKANTU_ERROR("Type " << type << "not implemented by MeshData."); + } +} +#undef AKANTU_MESH_NODAL_DATA_CASE_MACRO + +/* -------------------------------------------------------------------------- */ +/// Register new elemental data (and alloc data) +template +Array & MeshData::allocNodalData(const ID & name, UInt nb_components) { + auto dataset = + std::make_unique>(0, nb_components, T(), _id + ":" + name); + auto * dataset_typed = dataset.get(); + nodal_data[name] = std::move(dataset); + typecode_map[MeshDataType::_nodal][name] = getTypeCode(); + return *dataset_typed; +} + +/* -------------------------------------------------------------------------- */ +template +const Array & MeshData::getNodalData(const ID & name) const { + auto it = nodal_data.find(name); + if (it == nodal_data.end()) + AKANTU_EXCEPTION("No nodal dataset named " << name << " found."); + return dynamic_cast &>(*(it->second.get())); +} + +/* -------------------------------------------------------------------------- */ +// Get an existing elemental data +template +Array & MeshData::getNodalData(const ID & name, UInt nb_components) { + auto it = nodal_data.find(name); + if (it == nodal_data.end()) + return allocNodalData(name, nb_components); + + return dynamic_cast &>(*(it->second.get())); } /* -------------------------------------------------------------------------- */ template const ElementTypeMapArray & MeshData::getElementalData(const ID & name) const { auto it = elemental_data.find(name); if (it == elemental_data.end()) AKANTU_EXCEPTION("No dataset named " << name << " found."); - return dynamic_cast &>(*(it->second)); + return dynamic_cast &>(*(it->second.get())); } /* -------------------------------------------------------------------------- */ // Get an existing elemental data template ElementTypeMapArray & MeshData::getElementalData(const ID & name) { auto it = elemental_data.find(name); if (it == elemental_data.end()) AKANTU_EXCEPTION("No dataset named " << name << " found."); - return dynamic_cast &>(*(it->second)); + return dynamic_cast &>(*(it->second.get())); } /* -------------------------------------------------------------------------- */ template -bool MeshData::hasDataArray(const ID & name, const ElementType & elem_type, - const GhostType & ghost_type) const { +bool MeshData::hasData(const ID & name, const ElementType & elem_type, + const GhostType & ghost_type) const { auto it = elemental_data.find(name); if (it == elemental_data.end()) return false; auto & elem_map = dynamic_cast &>(*(it->second)); return elem_map.exists(elem_type, ghost_type); } /* -------------------------------------------------------------------------- */ -inline bool MeshData::hasData(const ID & name) const { - auto it = elemental_data.find(name); - return (it != elemental_data.end()); +inline bool MeshData::hasData(const ID & name, MeshDataType type) const { + if (type == MeshDataType::_elemental) { + auto it = elemental_data.find(name); + return (it != elemental_data.end()); + } + + if (type == MeshDataType::_nodal) { + auto it = nodal_data.find(name); + return (it != nodal_data.end()); + } + + return false; } /* -------------------------------------------------------------------------- */ template const Array & MeshData::getElementalDataArray(const ID & name, const ElementType & elem_type, const GhostType & ghost_type) const { auto it = elemental_data.find(name); if (it == elemental_data.end()) { AKANTU_EXCEPTION("Data named " << name << " not registered for type: " << elem_type << " - ghost_type:" << ghost_type << "!"); } return dynamic_cast &>(*(it->second))( elem_type, ghost_type); } template Array & MeshData::getElementalDataArray(const ID & name, const ElementType & elem_type, const GhostType & ghost_type) { auto it = elemental_data.find(name); if (it == elemental_data.end()) { AKANTU_EXCEPTION("Data named " << name << " not registered for type: " << elem_type << " - ghost_type:" << ghost_type << "!"); } - return dynamic_cast &>(*(it->second))(elem_type, - ghost_type); + return dynamic_cast &>(*(it->second.get()))( + elem_type, ghost_type); } /* -------------------------------------------------------------------------- */ // Get an elemental data array, if it does not exist: allocate it template Array & MeshData::getElementalDataArrayAlloc(const ID & name, const ElementType & elem_type, const GhostType & ghost_type, UInt nb_component) { auto it = elemental_data.find(name); ElementTypeMapArray * dataset; if (it == elemental_data.end()) { - dataset = allocElementalData(name); + dataset = &allocElementalData(name); } else { - dataset = dynamic_cast *>(it->second); + dataset = dynamic_cast *>(it->second.get()); } AKANTU_DEBUG_ASSERT( - getTypeCode() == typecode_map.find(name)->second, + getTypeCode() == + typecode_map.at(MeshDataType::_elemental).find(name)->second, "Function getElementalDataArrayAlloc called with the wrong type!"); if (!(dataset->exists(elem_type, ghost_type))) { dataset->alloc(0, nb_component, elem_type, ghost_type); } return (*dataset)(elem_type, ghost_type); } /* -------------------------------------------------------------------------- */ #define AKANTU_MESH_DATA_CASE_MACRO(r, name, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ nb_comp = getNbComponentTemplated( \ name, el_type, ghost_type); \ break; \ } inline UInt MeshData::getNbComponent(const ID & name, const ElementType & el_type, const GhostType & ghost_type) const { - auto it = typecode_map.find(name); + auto it = typecode_map.at(MeshDataType::_elemental).find(name); UInt nb_comp(0); - if (it == typecode_map.end()) { + if (it == typecode_map.at(MeshDataType::_elemental).end()) { AKANTU_EXCEPTION("Could not determine the type held in dataset " << name << " for type: " << el_type << " - ghost_type:" << ghost_type << "."); } MeshDataTypeCode type = it->second; switch (type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_MESH_DATA_CASE_MACRO, name, AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR( "Could not call the correct instance of getNbComponentTemplated."); break; } return nb_comp; } #undef AKANTU_MESH_DATA_CASE_MACRO /* -------------------------------------------------------------------------- */ template inline UInt MeshData::getNbComponentTemplated(const ID & name, const ElementType & el_type, const GhostType & ghost_type) const { return getElementalDataArray(name, el_type, ghost_type).getNbComponent(); } +/* -------------------------------------------------------------------------- */ +inline UInt MeshData::getNbComponent(const ID & name) const { + auto it = nodal_data.find(name); + if (it == nodal_data.end()) { + AKANTU_EXCEPTION("No nodal dataset registered with the name" << name + << "."); + } + + return it->second->getNbComponent(); +} + /* -------------------------------------------------------------------------- */ // get the names of the data stored in elemental_data #define AKANTU_MESH_DATA_CASE_MACRO(r, name, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ ElementTypeMapArray * dataset; \ dataset = \ dynamic_cast *>( \ - it->second); \ + it->second.get()); \ exists = dataset->exists(el_type, ghost_type); \ break; \ } -inline void MeshData::getTagNames(StringVector & tags, - const ElementType & el_type, +inline auto MeshData::getTagNames(const ElementType & el_type, const GhostType & ghost_type) const { - tags.clear(); + std::vector tags; bool exists(false); auto it = elemental_data.begin(); auto it_end = elemental_data.end(); for (; it != it_end; ++it) { MeshDataTypeCode type = getTypeCode(it->first); switch (type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_MESH_DATA_CASE_MACRO, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not determine the proper type to (dynamic-)cast."); break; } if (exists) { tags.push_back(it->first); } } + + return tags; } #undef AKANTU_MESH_DATA_CASE_MACRO +/* -------------------------------------------------------------------------- */ +inline auto MeshData::getTagNames() const { + std::vector tags; + for (auto && data : nodal_data) { + tags.push_back(std::get<0>(data)); + } + return tags; +} + /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_MESH_DATA_TMPL_HH__ */ diff --git a/src/mesh/mesh_inline_impl.cc b/src/mesh/mesh_inline_impl.cc index e730ce663..4f4b64503 100644 --- a/src/mesh/mesh_inline_impl.cc +++ b/src/mesh/mesh_inline_impl.cc @@ -1,693 +1,675 @@ /** * @file mesh_inline_impl.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Thu Jul 15 2010 * @date last modification: Mon Dec 18 2017 * * @brief Implementation of the inline functions of the mesh class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "element_class.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_INLINE_IMPL_CC__ #define __AKANTU_MESH_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline ElementKind Element::kind() const { return Mesh::getKind(type); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template Mesh::ElementTypesIteratorHelper Mesh::elementTypes(pack &&... _pack) const { return connectivities.elementTypes(_pack...); } /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh) : new_numbering(mesh.getNbNodes(), 1, "new_numbering") {} /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh, const ID & new_numbering_id) : new_numbering(new_numbering_id, mesh.getID(), mesh.getMemoryID()) {} /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewElementsEvent & event) { this->nodes_to_elements.resize(nodes->size()); for (const auto & elem : event.getList()) { const Array & conn = connectivities(elem.type, elem.ghost_type); UInt nb_nodes_per_elem = this->getNbNodesPerElement(elem.type); for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = conn(elem.element, n); if (not nodes_to_elements[node]) nodes_to_elements[node] = std::make_unique>(); nodes_to_elements[node]->insert(elem); } } EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewNodesEvent & event) { this->computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedElementsEvent & event) { this->connectivities.onElementsRemoved(event.getNewNumbering()); this->fillNodesToElements(); this->computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedNodesEvent & event) { const auto & new_numbering = event.getNewNumbering(); this->removeNodesFromArray(*nodes, new_numbering); if (nodes_global_ids and not mesh_parent) this->removeNodesFromArray(*nodes_global_ids, new_numbering); if (nodes_type and not mesh_parent) this->removeNodesFromArray(*nodes_type, new_numbering); if (not nodes_to_elements.empty()) { std::vector>> tmp( nodes_to_elements.size()); auto it = nodes_to_elements.begin(); UInt new_nb_nodes = 0; for (auto new_i : new_numbering) { if (new_i != UInt(-1)) { tmp[new_i] = std::move(*it); ++new_nb_nodes; } ++it; } tmp.resize(new_nb_nodes); std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin()); } computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template inline void Mesh::removeNodesFromArray(Array & vect, const Array & new_numbering) { Array tmp(vect.size(), vect.getNbComponent()); UInt nb_component = vect.getNbComponent(); UInt new_nb_nodes = 0; for (UInt i = 0; i < new_numbering.size(); ++i) { UInt new_i = new_numbering(i); if (new_i != UInt(-1)) { T * to_copy = vect.storage() + i * nb_component; std::uninitialized_copy(to_copy, to_copy + nb_component, tmp.storage() + new_i * nb_component); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if (not nodes_global_ids) { nodes_global_ids = std::make_shared>( nodes->size(), 1, getID() + ":nodes_global_ids"); for (auto && global_ids : enumerate(*nodes_global_ids)) { std::get<1>(global_ids) = std::get<0>(global_ids); } } AKANTU_DEBUG_OUT(); return *nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getNodesTypePointer() { AKANTU_DEBUG_IN(); if (not nodes_type) { nodes_type = std::make_shared>(nodes->size(), 1, _nt_normal); } AKANTU_DEBUG_OUT(); return *nodes_type; } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getConnectivityPointer(const ElementType & type, const GhostType & ghost_type) { if (connectivities.exists(type, ghost_type)) return connectivities(type, ghost_type); if (ghost_type != _not_ghost) { ghosts_counters.alloc(0, 1, type, ghost_type, 1); } AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); return connectivities.alloc(0, nb_nodes_per_element, type, ghost_type); } /* -------------------------------------------------------------------------- */ inline Array> & Mesh::getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type) { return getDataPointer>("element_to_subelement", type, ghost_type, 1, true); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type) { auto & array = getDataPointer( "subelement_to_element", type, ghost_type, getNbFacetsPerElement(type), true, is_mesh_facets, ElementNull); return array; } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement() const { return getData>("element_to_subelement"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) const { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(const Element & element) const { return getData>("element_to_subelement")(element); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelement(const Element & element) { return getData>("element_to_subelement")(element); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement() const { return getData("subelement_to_element"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) const { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElement(const Element & element) const { const auto & sub_to_element = this->getSubelementToElement(element.type, element.ghost_type); auto it = sub_to_element.begin(sub_to_element.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElement(const Element & element) { auto & sub_to_element = this->getSubelementToElement(element.type, element.ghost_type); auto it = sub_to_element.begin(sub_to_element.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent) { - Array & tmp = mesh_data.getElementalDataArrayAlloc( + Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) tmp.resize(mesh_parent->getNbElement(el_type, ghost_type)); else tmp.resize(this->getNbElement(el_type, ghost_type)); } else { tmp.resize(0); } return tmp; } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_) { - Array & tmp = mesh_data.getElementalDataArrayAlloc( + Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) tmp.resize(mesh_parent->getNbElement(el_type, ghost_type), defaul_); else tmp.resize(this->getNbElement(el_type, ghost_type), defaul_); } else { tmp.resize(0); } return tmp; } -/* -------------------------------------------------------------------------- */ -template -inline bool Mesh::hasData(const ID & data_name, const ElementType & el_type, - const GhostType & ghost_type) const { - return mesh_data.hasDataArray(data_name, el_type, ghost_type); -} - /* -------------------------------------------------------------------------- */ template inline const Array & Mesh::getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type) const { - return mesh_data.getElementalDataArray(data_name, el_type, ghost_type); + return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type) { - return mesh_data.getElementalDataArray(data_name, el_type, ghost_type); -} - -/* -------------------------------------------------------------------------- */ -inline bool Mesh::hasData(const ID & data_name) const { - return mesh_data.hasData(data_name); + return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline const ElementTypeMapArray & Mesh::getData(const ID & data_name) const { - return mesh_data.getElementalData(data_name); + return this->getElementalData(data_name); } /* -------------------------------------------------------------------------- */ template inline ElementTypeMapArray & Mesh::getData(const ID & data_name) { - return mesh_data.getElementalData(data_name); + return this->getElementalData(data_name); } -/* -------------------------------------------------------------------------- */ -template -inline ElementTypeMapArray & Mesh::registerData(const ID & data_name) { - this->mesh_data.registerElementalData(data_name); - return this->getData(data_name); -} /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const ElementType & type, const GhostType & ghost_type) const { try { const Array & conn = connectivities(type, ghost_type); return conn.size(); } catch (...) { return 0; } } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & kind) const { AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == UInt(-1), "spatial_dimension is " << spatial_dimension << " and is greater than 3 !"); UInt nb_element = 0; for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) nb_element += getNbElement(type, ghost_type); return nb_element; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(const Element & element, Vector & barycenter) const { Vector conn = getConnectivity(element); Matrix local_coord(spatial_dimension, conn.size()); auto node_begin = make_view(*nodes, spatial_dimension).begin(); for(auto && node : enumerate(conn)) { local_coord(std::get<0>(node)) = Vector(node_begin[std::get<1>(node)]); } Math::barycenter(local_coord.storage(), conn.size(), spatial_dimension, barycenter.storage()); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElement(const ElementType & type) { UInt nb_nodes_per_element = 0; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass::getNbNodesPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT); #undef GET_NB_NODES_PER_ELEMENT return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getP1ElementType(const ElementType & type) { ElementType p1_type = _not_defined; #define GET_P1_TYPE(type) p1_type = ElementClass::getP1ElementType() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_P1_TYPE); #undef GET_P1_TYPE return p1_type; } /* -------------------------------------------------------------------------- */ inline ElementKind Mesh::getKind(const ElementType & type) { ElementKind kind = _ek_not_defined; #define GET_KIND(type) kind = ElementClass::getKind() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_KIND); #undef GET_KIND return kind; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getSpatialDimension(const ElementType & type) { UInt spatial_dimension = 0; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass::getSpatialDimension() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION); #undef GET_SPATIAL_DIMENSION return spatial_dimension; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetTypes(const ElementType & type, __attribute__((unused)) UInt t) { UInt nb = 0; #define GET_NB_FACET_TYPE(type) nb = ElementClass::getNbFacetTypes() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET_TYPE); #undef GET_NB_FACET_TYPE return nb; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getFacetType(const ElementType & type, UInt t) { #define GET_FACET_TYPE(type) return ElementClass::getFacetType(t); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return _not_defined; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getAllFacetTypes(const ElementType & type) { #define GET_FACET_TYPE(type) return ElementClass::getFacetTypes(); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return ElementClass<_not_defined>::getFacetTypes(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) n_facet = ElementClass::getNbFacetsPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetLocalConnectivity(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); #define GET_FACET_CON(type) \ AKANTU_DEBUG_OUT(); \ return ElementClass::getFacetLocalConnectivityPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON); #undef GET_FACET_CON AKANTU_DEBUG_OUT(); return ElementClass<_not_defined>::getFacetLocalConnectivityPerElement(0); // This avoid a compilation warning but will certainly // also cause a segfault if reached } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetConnectivity(const Element & element, UInt t) const { AKANTU_DEBUG_IN(); Matrix local_facets(getFacetLocalConnectivity(element.type, t)); Matrix facets(local_facets.rows(), local_facets.cols()); const Array & conn = connectivities(element.type, element.ghost_type); for (UInt f = 0; f < facets.rows(); ++f) { for (UInt n = 0; n < facets.cols(); ++n) { facets(f, n) = conn(element.element, local_facets(f, n)); } } AKANTU_DEBUG_OUT(); return facets; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivity(const Element & element) const { const auto & conn = connectivities(element.type, element.ghost_type); auto it = conn.begin(conn.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivity(const Element & element) { auto & conn = connectivities(element.type, element.ghost_type); auto it = conn.begin(conn.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ template inline void Mesh::extractNodalValuesFromElement( const Array & nodal_values, T * local_coord, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const { for (UInt n = 0; n < n_nodes; ++n) { memcpy(local_coord + n * nb_degree_of_freedom, nodal_values.storage() + connectivity[n] * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(T)); } } /* -------------------------------------------------------------------------- */ inline void Mesh::addConnectivityType(const ElementType & type, const GhostType & ghost_type) { getConnectivityPointer(type, ghost_type); } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPureGhostNode(UInt n) const { return nodes_type ? ((*nodes_type)(n) == _nt_pure_ghost) : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return nodes_type ? ((*nodes_type)(n) == _nt_master) || ((*nodes_type)(n) == _nt_normal) : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == _nt_normal : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == _nt_master : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return nodes_type ? (*nodes_type)(n) >= 0 : false; } /* -------------------------------------------------------------------------- */ inline NodeType Mesh::getNodeType(UInt local_id) const { return nodes_type ? (*nodes_type)(local_id) : _nt_normal; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeLocalId(UInt global_id) const { AKANTU_DEBUG_ASSERT(nodes_global_ids != nullptr, "The global ids are note set."); return nodes_global_ids->find(global_id); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->size(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElementList(const Array & elements) { UInt nb_nodes_per_element = 0; UInt nb_nodes = 0; ElementType current_element_type = _not_defined; Array::const_iterator el_it = elements.begin(); Array::const_iterator el_end = elements.end(); for (; el_it != el_end; ++el_it) { const Element & el = *el_it; if (el.type != current_element_type) { current_element_type = el.type; nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type); } nb_nodes += nb_nodes_per_element; } return nb_nodes; } /* -------------------------------------------------------------------------- */ inline Mesh & Mesh::getMeshFacets() { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshFacets() const { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshParent() const { if (!this->mesh_parent) AKANTU_SILENT_EXCEPTION( "No parent mesh is defined! This is only valid in a mesh_facets"); return *this->mesh_parent; } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_MESH_INLINE_IMPL_CC__ */ diff --git a/src/mesh_utils/mesh_utils_distribution.cc b/src/mesh_utils/mesh_utils_distribution.cc index 416eb5c45..59bd034e8 100644 --- a/src/mesh_utils/mesh_utils_distribution.cc +++ b/src/mesh_utils/mesh_utils_distribution.cc @@ -1,174 +1,176 @@ /** * @file mesh_utils_distribution.cc * * @author Nicolas Richart * * @date creation: Tue Nov 08 2016 * @date last modification: Tue Nov 07 2017 * * @brief Implementation of the methods of mesh utils distribute * * @section LICENSE * * 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 "mesh_utils_distribution.hh" #include "element_info_per_processor.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_partition.hh" #include "mesh_utils.hh" #include "node_info_per_processor.hh" #include "node_synchronizer.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace MeshUtilsDistribution { /* ------------------------------------------------------------------------ */ void distributeMeshCentralized(Mesh & mesh, UInt, const MeshPartition & partition) { MeshAccessor mesh_accessor(mesh); ElementSynchronizer & element_synchronizer = mesh_accessor.getElementSynchronizer(); NodeSynchronizer & node_synchronizer = mesh_accessor.getNodeSynchronizer(); const Communicator & comm = element_synchronizer.getCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); mesh_accessor.setNbGlobalNodes(mesh.getNbNodes()); auto & gids = mesh_accessor.getNodesGlobalIds(); if (nb_proc == 1) return; gids.resize(0); mesh.synchronizeGroupNames(); AKANTU_DEBUG_ASSERT( partition.getNbPartition() == nb_proc, "The number of partition does not match the number of processors: " << partition.getNbPartition() << " != " << nb_proc); /** * connectivity and communications scheme construction */ UInt count = 0; /* --- MAIN LOOP ON TYPES --- */ for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { /// \todo change this ugly way to avoid a problem if an element /// type is present in the mesh but not in the partitions try { partition.getPartition(type, _not_ghost); } catch (...) { continue; } MasterElementInfoPerProc proc_infos(element_synchronizer, count, my_rank, type, partition); proc_infos.synchronizeConnectivities(); proc_infos.synchronizePartitions(); proc_infos.synchronizeTags(); proc_infos.synchronizeGroups(); ++count; } { /// Ending the synchronization of elements by sending a stop message MasterElementInfoPerProc proc_infos(element_synchronizer, count, my_rank, _not_defined, partition); ++count; } /** * Nodes synchronization */ MasterNodeInfoPerProc node_proc_infos(node_synchronizer, count, my_rank); node_proc_infos.synchronizeNodes(); node_proc_infos.synchronizeTypes(); node_proc_infos.synchronizeGroups(); + node_proc_infos.synchronizeTags(); MeshUtils::fillElementToSubElementsData(mesh); mesh_accessor.setDistributed(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ void distributeMeshCentralized(Mesh & mesh, UInt root) { MeshAccessor mesh_accessor(mesh); ElementSynchronizer & element_synchronizer = mesh_accessor.getElementSynchronizer(); NodeSynchronizer & node_synchronizer = mesh_accessor.getNodeSynchronizer(); const Communicator & comm = element_synchronizer.getCommunicator(); UInt nb_proc = comm.getNbProc(); mesh_accessor.getNodesGlobalIds().resize(0); if (nb_proc == 1) return; mesh.synchronizeGroupNames(); /** * connectivity and communications scheme construction on distant * processors */ UInt count = 0; bool need_synchronize = true; do { /* --------<<<<-SIZE--------------------------------------------------- */ SlaveElementInfoPerProc proc_infos(element_synchronizer, count, root); need_synchronize = proc_infos.needSynchronize(); if (need_synchronize) { proc_infos.synchronizeConnectivities(); proc_infos.synchronizePartitions(); proc_infos.synchronizeTags(); proc_infos.synchronizeGroups(); } ++count; } while (need_synchronize); /** * Nodes synchronization */ SlaveNodeInfoPerProc node_proc_infos(node_synchronizer, count, root); node_proc_infos.synchronizeNodes(); node_proc_infos.synchronizeTypes(); node_proc_infos.synchronizeGroups(); + node_proc_infos.synchronizeTags(); MeshUtils::fillElementToSubElementsData(mesh); mesh_accessor.setDistributed(); } } // namespace MeshUtilsDistribution } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc index a422e9b0d..29bd66b81 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc @@ -1,179 +1,179 @@ /** * @file embedded_interface_intersector.cc * * @author Lucas Frerot * * @date creation: Fri May 01 2015 * @date last modification: Tue Feb 20 2018 * * @brief Class that loads the interface from mesh and computes intersections * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "embedded_interface_intersector.hh" #include "mesh_segment_intersector.hh" /// Helper macro for types in the mesh. Creates an intersector and computes /// intersection queries #define INTERFACE_INTERSECTOR_CASE(dim, type) \ do { \ MeshSegmentIntersector intersector(this->mesh, interface_mesh); \ name_to_primitives_it = name_to_primitives_map.begin(); \ for (; name_to_primitives_it != name_to_primitives_end; \ ++name_to_primitives_it) { \ intersector.setPhysicalName(name_to_primitives_it->first); \ intersector.buildResultFromQueryList(name_to_primitives_it->second); \ } \ } while (0) #define INTERFACE_INTERSECTOR_CASE_2D(type) INTERFACE_INTERSECTOR_CASE(2, type) #define INTERFACE_INTERSECTOR_CASE_3D(type) INTERFACE_INTERSECTOR_CASE(3, type) namespace akantu { EmbeddedInterfaceIntersector::EmbeddedInterfaceIntersector( Mesh & mesh, const Mesh & primitive_mesh) : MeshGeomAbstract(mesh), interface_mesh(mesh.getSpatialDimension(), "interface_mesh"), primitive_mesh(primitive_mesh) { // Initiating mesh connectivity and data interface_mesh.addConnectivityType(_segment_2, _not_ghost); interface_mesh.addConnectivityType(_segment_2, _ghost); - interface_mesh.registerData("associated_element") + interface_mesh.registerElementalData("associated_element") .alloc(0, 1, _segment_2); - interface_mesh.registerData("physical_names") + interface_mesh.registerElementalData("physical_names") .alloc(0, 1, _segment_2); } EmbeddedInterfaceIntersector::~EmbeddedInterfaceIntersector() {} void EmbeddedInterfaceIntersector::constructData(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); const UInt dim = this->mesh.getSpatialDimension(); if (dim == 1) AKANTU_ERROR( "No embedded model in 1D. Deactivate intersection initialization"); Array * physical_names = NULL; try { physical_names = &const_cast &>( this->primitive_mesh.getData("physical_names", _segment_2)); } catch (debug::Exception & e) { AKANTU_ERROR("You must define physical names to reinforcements in " "order to use the embedded model"); throw e; } const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_segment_2); Array::const_vector_iterator connectivity = primitive_mesh.getConnectivity(_segment_2).begin(nb_nodes_per_element); Array::scalar_iterator names_it = physical_names->begin(), names_end = physical_names->end(); std::map> name_to_primitives_map; // Loop over the physical names and register segment lists in // name_to_primitives_map for (; names_it != names_end; ++names_it) { UInt element_id = names_it - physical_names->begin(); const Vector el_connectivity = connectivity[element_id]; K::Segment_3 segment = this->createSegment(el_connectivity); name_to_primitives_map[*names_it].push_back(segment); } // Loop over the background types of the mesh Mesh::type_iterator type_it = this->mesh.firstType(dim, _not_ghost), type_end = this->mesh.lastType(dim, _not_ghost); std::map>::iterator name_to_primitives_it, name_to_primitives_end = name_to_primitives_map.end(); for (; type_it != type_end; ++type_it) { // Used in AKANTU_BOOST_ELEMENT_SWITCH ElementType type = *type_it; AKANTU_DEBUG_INFO("Computing intersections with background element type " << type); switch (dim) { case 1: break; case 2: // Compute intersections for supported 2D elements AKANTU_BOOST_ELEMENT_SWITCH(INTERFACE_INTERSECTOR_CASE_2D, (_triangle_3)(_triangle_6)); break; case 3: // Compute intersections for supported 3D elements AKANTU_BOOST_ELEMENT_SWITCH(INTERFACE_INTERSECTOR_CASE_3D, (_tetrahedron_4)); break; } } AKANTU_DEBUG_OUT(); } K::Segment_3 EmbeddedInterfaceIntersector::createSegment(const Vector & connectivity) { AKANTU_DEBUG_IN(); K::Point_3 *source = NULL, *target = NULL; const Array & nodes = this->primitive_mesh.getNodes(); if (this->mesh.getSpatialDimension() == 2) { source = new K::Point_3(nodes(connectivity(0), 0), nodes(connectivity(0), 1), 0.); target = new K::Point_3(nodes(connectivity(1), 0), nodes(connectivity(1), 1), 0.); } else if (this->mesh.getSpatialDimension() == 3) { source = new K::Point_3(nodes(connectivity(0), 0), nodes(connectivity(0), 1), nodes(connectivity(0), 2)); target = new K::Point_3(nodes(connectivity(1), 0), nodes(connectivity(1), 1), nodes(connectivity(1), 2)); } K::Segment_3 segment(*source, *target); delete source; delete target; AKANTU_DEBUG_OUT(); return segment; } } // namespace akantu #undef INTERFACE_INTERSECTOR_CASE #undef INTERFACE_INTERSECTOR_CASE_2D #undef INTERFACE_INTERSECTOR_CASE_3D diff --git a/src/synchronizer/communication_buffer.hh b/src/synchronizer/communication_buffer.hh index 5e0b0da8d..737ab407a 100644 --- a/src/synchronizer/communication_buffer.hh +++ b/src/synchronizer/communication_buffer.hh @@ -1,182 +1,185 @@ /** * @file communication_buffer.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 08 2017 * * @brief Buffer for packing and unpacking data * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" #include "element.hh" #ifndef __AKANTU_COMMUNICATION_BUFFER_HH__ #define __AKANTU_COMMUNICATION_BUFFER_HH__ namespace akantu { template class CommunicationBufferTemplated { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: explicit CommunicationBufferTemplated(UInt size) : buffer(size, 1) { ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); }; CommunicationBufferTemplated() : CommunicationBufferTemplated(0) {} CommunicationBufferTemplated(const CommunicationBufferTemplated & other) { buffer = other.buffer; ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); } CommunicationBufferTemplated & operator=(const CommunicationBufferTemplated & other) { if (this != &other) { buffer = other.buffer; ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); } return *this; } virtual ~CommunicationBufferTemplated() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// reset to "empty" inline void reset(); /// resize the internal buffer inline void resize(UInt size); + /// resize the internal buffer allocate always + inline void reserve(UInt size); + /// clear buffer context inline void clear(); private: inline void packResize(UInt size); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: inline char * storage() { return buffer.storage(); }; inline const char * storage() const { return buffer.storage(); }; /* ------------------------------------------------------------------------ */ /* Operators */ /* ------------------------------------------------------------------------ */ public: /// printing tool template inline std::string extractStream(UInt packet_size); /// packing data template inline CommunicationBufferTemplated & operator<<(const T & to_pack); template inline CommunicationBufferTemplated & operator<<(const Vector & to_pack); template inline CommunicationBufferTemplated & operator<<(const Matrix & to_pack); template inline CommunicationBufferTemplated & operator<<(const std::vector & to_pack); /// unpacking data template inline CommunicationBufferTemplated & operator>>(T & to_unpack); template inline CommunicationBufferTemplated & operator>>(Vector & to_unpack); template inline CommunicationBufferTemplated & operator>>(Matrix & to_unpack); template inline CommunicationBufferTemplated & operator>>(std::vector & to_unpack); inline CommunicationBufferTemplated & operator<<(const std::string & to_pack); inline CommunicationBufferTemplated & operator>>(std::string & to_unpack); private: template inline void packIterable(T & to_pack); template inline void unpackIterable(T & to_pack); /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: template static inline UInt sizeInBuffer(const T & data); template static inline UInt sizeInBuffer(const Vector & data); template static inline UInt sizeInBuffer(const Matrix & data); template static inline UInt sizeInBuffer(const std::vector & data); static inline UInt sizeInBuffer(const std::string & data); /// return the size in bytes of the stored values inline UInt getPackedSize() const { return ptr_pack - buffer.storage(); }; /// return the size in bytes of data left to be unpacked inline UInt getLeftToUnpack() const { return buffer.size() - (ptr_unpack - buffer.storage()); }; /// return the global size allocated inline UInt size() const { return buffer.size(); }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// current position for packing char * ptr_pack; /// current position for unpacking char * ptr_unpack; /// storing buffer Array buffer; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined(AKANTU_INCLUDE_INLINE_IMPL) #include "communication_buffer_inline_impl.cc" #endif using CommunicationBuffer = CommunicationBufferTemplated; using DynamicCommunicationBuffer = CommunicationBufferTemplated; } // namespace akantu #endif /* __AKANTU_COMMUNICATION_BUFFER_HH__ */ diff --git a/src/synchronizer/communication_buffer_inline_impl.cc b/src/synchronizer/communication_buffer_inline_impl.cc index d242a4b42..f434d3af0 100644 --- a/src/synchronizer/communication_buffer_inline_impl.cc +++ b/src/synchronizer/communication_buffer_inline_impl.cc @@ -1,318 +1,329 @@ /** * @file communication_buffer_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Thu Apr 14 2011 * @date last modification: Wed Nov 08 2017 * * @brief CommunicationBuffer inline implementation * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ template template inline UInt CommunicationBufferTemplated::sizeInBuffer(const T &) { return sizeof(T); } template template inline UInt CommunicationBufferTemplated::sizeInBuffer(const Vector & data) { UInt size = data.size() * sizeof(T); return size; } template template inline UInt CommunicationBufferTemplated::sizeInBuffer(const Matrix & data) { UInt size = data.size() * sizeof(T); return size; } template template inline UInt CommunicationBufferTemplated::sizeInBuffer( const std::vector & data) { UInt size = data.size() * sizeof(T) + sizeof(size_t); return size; } template inline UInt CommunicationBufferTemplated::sizeInBuffer( const std::string & data) { UInt size = data.size() * sizeof(std::string::value_type) + sizeof(size_t); return size; } /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::packResize(UInt size) { if (!is_static) { char * values = buffer.storage(); buffer.resize(buffer.size() + size); ptr_pack = buffer.storage() + (ptr_pack - values); ptr_unpack = buffer.storage() + (ptr_unpack - values); } } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator<<(const T & to_pack) { UInt size = sizeInBuffer(to_pack); packResize(size); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.size()) >= (ptr_pack + size), "Packing too much data in the CommunicationBufferTemplated"); memcpy(ptr_pack, reinterpret_cast(&to_pack), size); ptr_pack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(T & to_unpack) { UInt size = sizeInBuffer(to_unpack); alignas(alignof(T)) std::array aligned_ptr; memcpy(aligned_ptr.data(), ptr_unpack, size); auto * tmp = reinterpret_cast(aligned_ptr.data()); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.size()) >= (ptr_unpack + size), "Unpacking too much data in the CommunicationBufferTemplated"); to_unpack = *tmp; // memcpy(reinterpret_cast(&to_unpack), ptr_unpack, size); ptr_unpack += size; return *this; } /* -------------------------------------------------------------------------- */ /* Specialization */ /* -------------------------------------------------------------------------- */ /** * Vector */ /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator<<(const Vector & to_pack) { UInt size = sizeInBuffer(to_pack); packResize(size); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.size()) >= (ptr_pack + size), "Packing too much data in the CommunicationBufferTemplated"); memcpy(ptr_pack, to_pack.storage(), size); ptr_pack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(Vector & to_unpack) { UInt size = sizeInBuffer(to_unpack); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.size()) >= (ptr_unpack + size), "Unpacking too much data in the CommunicationBufferTemplated"); memcpy(to_unpack.storage(), ptr_unpack, size); ptr_unpack += size; return *this; } /** * Matrix */ /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator<<(const Matrix & to_pack) { UInt size = sizeInBuffer(to_pack); packResize(size); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.size()) >= (ptr_pack + size), "Packing too much data in the CommunicationBufferTemplated"); memcpy(ptr_pack, to_pack.storage(), size); ptr_pack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(Matrix & to_unpack) { UInt size = sizeInBuffer(to_unpack); AKANTU_DEBUG_ASSERT( (buffer.storage() + buffer.size()) >= (ptr_unpack + size), "Unpacking too much data in the CommunicationBufferTemplated"); memcpy(to_unpack.storage(), ptr_unpack, size); ptr_unpack += size; return *this; } /* -------------------------------------------------------------------------- */ template template inline void CommunicationBufferTemplated::packIterable(T & to_pack) { operator<<(size_t(to_pack.size())); auto it = to_pack.begin(); auto end = to_pack.end(); for (; it != end; ++it) operator<<(*it); } /* -------------------------------------------------------------------------- */ template template inline void CommunicationBufferTemplated::unpackIterable(T & to_unpack) { size_t size; operator>>(size); to_unpack.resize(size); auto it = to_unpack.begin(); auto end = to_unpack.end(); for (; it != end; ++it) operator>>(*it); } /** * std::vector */ /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated:: operator<<(const std::vector & to_pack) { packIterable(to_pack); return *this; } /* -------------------------------------------------------------------------- */ template template inline CommunicationBufferTemplated & CommunicationBufferTemplated:: operator>>(std::vector & to_unpack) { unpackIterable(to_unpack); return *this; } /** * std::string */ /* -------------------------------------------------------------------------- */ template inline CommunicationBufferTemplated & CommunicationBufferTemplated:: operator<<(const std::string & to_pack) { packIterable(to_pack); return *this; } /* -------------------------------------------------------------------------- */ template inline CommunicationBufferTemplated & CommunicationBufferTemplated::operator>>(std::string & to_unpack) { unpackIterable(to_unpack); return *this; } /* -------------------------------------------------------------------------- */ template template inline std::string CommunicationBufferTemplated::extractStream(UInt block_size) { std::stringstream str; auto * ptr = reinterpret_cast(buffer.storage()); UInt sz = buffer.size() / sizeof(T); UInt sz_block = block_size / sizeof(T); UInt n_block = 0; for (UInt i = 0; i < sz; ++i) { if (i % sz_block == 0) { str << std::endl << n_block << " "; ++n_block; } str << *ptr << " "; ++ptr; } return str.str(); } /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::resize(UInt size) { if (!is_static) { buffer.resize(0); } else { buffer.resize(size); } reset(); #ifndef AKANTU_NDEBUG clear(); #endif } +/* -------------------------------------------------------------------------- */ +template +inline void CommunicationBufferTemplated::reserve(UInt size) { + char * values = buffer.storage(); + auto nb_packed = ptr_pack - values; + + buffer.resize(size); + ptr_pack = buffer.storage() + nb_packed; + ptr_unpack = buffer.storage() + (ptr_unpack - values); +} + /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::clear() { buffer.clear(); } /* -------------------------------------------------------------------------- */ template inline void CommunicationBufferTemplated::reset() { ptr_pack = buffer.storage(); ptr_unpack = buffer.storage(); } /* -------------------------------------------------------------------------- */ // template // inline CommunicationBufferTemplated & // CommunicationBufferTemplated::packMeshData (const MeshData & // to_pack, const ElementType & type) { // UInt size = to_pack.size(); // operator<<(size); // typename std::vector::iterator it = to_pack.begin(); // typename std::vector::iterator end = to_pack.end(); // for(;it != end; ++it) operator<<(*it); // return *this; //} diff --git a/src/synchronizer/communicator.hh b/src/synchronizer/communicator.hh index e3ac8f462..ff6119989 100644 --- a/src/synchronizer/communicator.hh +++ b/src/synchronizer/communicator.hh @@ -1,436 +1,445 @@ /** * @file communicator.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 15 2017 * * @brief Class handling the parallel communications * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" #include "aka_event_handler_manager.hh" #include "communication_buffer.hh" #include "communication_request.hh" #include "communicator_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STATIC_COMMUNICATOR_HH__ #define __AKANTU_STATIC_COMMUNICATOR_HH__ namespace akantu { namespace debug { class CommunicationException : public Exception { public: CommunicationException() : Exception("An exception happen during a communication process.") {} }; } // namespace debug /// @enum SynchronizerOperation reduce operation that the synchronizer can /// perform enum class SynchronizerOperation { _sum, _min, _max, _prod, _land, _band, _lor, _bor, _lxor, _bxor, _min_loc, _max_loc, _null }; enum class CommunicationMode { _auto, _synchronous, _ready }; namespace { int _any_source = -1; } } // namespace akantu namespace akantu { struct CommunicatorInternalData { virtual ~CommunicatorInternalData() = default; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class Communicator : public EventHandlerManager { struct private_member {}; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Communicator(int & argc, char **& argv, const private_member &); ~Communicator() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Point to Point */ /* ------------------------------------------------------------------------ */ template inline void receive(Array & values, Int sender, Int tag) const { return this->receiveImpl( values.storage(), values.size() * values.getNbComponent(), sender, tag); } template inline void receive(Tensor & values, Int sender, Int tag, std::enable_if_t::value> * = nullptr) const { return this->receiveImpl(values.storage(), values.size(), sender, tag); } - template - inline void receive(CommunicationBufferTemplated & values, + + inline void receive(CommunicationBufferTemplated & values, Int sender, Int tag) const { return this->receiveImpl(values.storage(), values.size(), sender, tag); } + + inline void receive(CommunicationBufferTemplated & values, + Int sender, Int tag) const { + CommunicationStatus status; + this->probe(sender, tag, status); + values.reserve(status.size()); + return this->receiveImpl(values.storage(), values.size(), sender, tag); + } template inline void receive(T & values, Int sender, Int tag, std::enable_if_t::value> * = nullptr) const { return this->receiveImpl(&values, 1, sender, tag); } /* ------------------------------------------------------------------------ */ template inline void send(const Array & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto) const { return this->sendImpl(values.storage(), values.size() * values.getNbComponent(), receiver, tag, mode); } template inline void send(const Tensor & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto, std::enable_if_t::value> * = nullptr) const { return this->sendImpl(values.storage(), values.size(), receiver, tag, mode); } + template inline void send(const CommunicationBufferTemplated & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto) const { return this->sendImpl(values.storage(), values.size(), receiver, tag, mode); } template inline void send(const T & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto, std::enable_if_t::value> * = nullptr) const { return this->sendImpl(&values, 1, receiver, tag, mode); } /* ------------------------------------------------------------------------ */ template inline CommunicationRequest asyncSend(const Array & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto) const { return this->asyncSendImpl(values.storage(), values.size() * values.getNbComponent(), receiver, tag, mode); } template inline CommunicationRequest asyncSend(const std::vector & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto) const { return this->asyncSendImpl(values.data(), values.size(), receiver, tag, mode); } template inline CommunicationRequest asyncSend(const Tensor & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto, std::enable_if_t::value> * = nullptr) const { return this->asyncSendImpl(values.storage(), values.size(), receiver, tag, mode); } template inline CommunicationRequest asyncSend(const CommunicationBufferTemplated & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto) const { return this->asyncSendImpl(values.storage(), values.size(), receiver, tag, mode); } template inline CommunicationRequest asyncSend(const T & values, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto, std::enable_if_t::value> * = nullptr) const { return this->asyncSendImpl(&values, 1, receiver, tag, mode); } /* ------------------------------------------------------------------------ */ template inline CommunicationRequest asyncReceive(Array & values, Int sender, Int tag) const { return this->asyncReceiveImpl( values.storage(), values.size() * values.getNbComponent(), sender, tag); } template inline CommunicationRequest asyncReceive(std::vector & values, Int sender, Int tag) const { return this->asyncReceiveImpl(values.data(), values.size(), sender, tag); } template ::value>> inline CommunicationRequest asyncReceive(Tensor & values, Int sender, Int tag) const { return this->asyncReceiveImpl(values.storage(), values.size(), sender, tag); } template inline CommunicationRequest asyncReceive(CommunicationBufferTemplated & values, Int sender, Int tag) const { return this->asyncReceiveImpl(values.storage(), values.size(), sender, tag); } /* ------------------------------------------------------------------------ */ template void probe(Int sender, Int tag, CommunicationStatus & status) const; template bool asyncProbe(Int sender, Int tag, CommunicationStatus & status) const; /* ------------------------------------------------------------------------ */ /* Collectives */ /* ------------------------------------------------------------------------ */ template inline void allReduce(Array & values, const SynchronizerOperation & op) const { this->allReduceImpl(values.storage(), values.size() * values.getNbComponent(), op); } template inline void allReduce(Tensor & values, const SynchronizerOperation & op, std::enable_if_t::value> * = nullptr) const { this->allReduceImpl(values.storage(), values.size(), op); } template inline void allReduce(T & values, const SynchronizerOperation & op, std::enable_if_t::value> * = nullptr) const { this->allReduceImpl(&values, 1, op); } /* ------------------------------------------------------------------------ */ template inline void allGather(Array & values) const { AKANTU_DEBUG_ASSERT(UInt(psize) == values.size(), "The array size is not correct"); this->allGatherImpl(values.storage(), values.getNbComponent()); } template ::value>> inline void allGather(Tensor & values) const { AKANTU_DEBUG_ASSERT(values.size() / UInt(psize) > 0, "The vector size is not correct"); this->allGatherImpl(values.storage(), values.size() / UInt(psize)); } /* ------------------------------------------------------------------------ */ template inline void allGatherV(Array & values, const Array & sizes) const { this->allGatherVImpl(values.storage(), sizes.storage()); } /* ------------------------------------------------------------------------ */ template inline void reduce(Array & values, const SynchronizerOperation & op, int root = 0) const { this->reduceImpl(values.storage(), values.size() * values.getNbComponent(), op, root); } /* ------------------------------------------------------------------------ */ template inline void gather(Tensor & values, int root = 0, std::enable_if_t::value> * = nullptr) const { this->gatherImpl(values.storage(), values.getNbComponent(), root); } template inline void gather(T values, int root = 0, std::enable_if_t::value> * = nullptr) const { this->gatherImpl(&values, 1, root); } /* ------------------------------------------------------------------------ */ template inline void gather(Tensor & values, Array & gathered, std::enable_if_t::value> * = nullptr) const { AKANTU_DEBUG_ASSERT(values.size() == gathered.getNbComponent(), "The array size is not correct"); gathered.resize(psize); this->gatherImpl(values.data(), values.size(), gathered.storage(), gathered.getNbComponent()); } template inline void gather(T values, Array & gathered, std::enable_if_t::value> * = nullptr) const { this->gatherImpl(&values, 1, gathered.storage(), 1); } /* ------------------------------------------------------------------------ */ template inline void gatherV(Array & values, const Array & sizes, int root = 0) const { this->gatherVImpl(values.storage(), sizes.storage(), root); } /* ------------------------------------------------------------------------ */ template inline void broadcast(Array & values, int root = 0) const { this->broadcastImpl(values.storage(), values.size() * values.getNbComponent(), root); } template inline void broadcast(CommunicationBufferTemplated & values, int root = 0) const { this->broadcastImpl(values.storage(), values.size(), root); } template inline void broadcast(T & values, int root = 0) const { this->broadcastImpl(&values, 1, root); } /* ------------------------------------------------------------------------ */ void barrier() const; CommunicationRequest asyncBarrier() const; /* ------------------------------------------------------------------------ */ /* Request handling */ /* ------------------------------------------------------------------------ */ bool test(CommunicationRequest & request) const; bool testAll(std::vector & request) const; void wait(CommunicationRequest & request) const; void waitAll(std::vector & requests) const; UInt waitAny(std::vector & requests) const; inline void freeCommunicationRequest(CommunicationRequest & request) const; inline void freeCommunicationRequest(std::vector & requests) const; template inline void receiveAnyNumber(std::vector & send_requests, MsgProcessor && processor, Int tag) const; protected: template void sendImpl(const T * buffer, Int size, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto) const; template void receiveImpl(T * buffer, Int size, Int sender, Int tag) const; template CommunicationRequest asyncSendImpl( const T * buffer, Int size, Int receiver, Int tag, const CommunicationMode & mode = CommunicationMode::_auto) const; template CommunicationRequest asyncReceiveImpl(T * buffer, Int size, Int sender, Int tag) const; template void allReduceImpl(T * values, int nb_values, const SynchronizerOperation & op) const; template void allGatherImpl(T * values, int nb_values) const; template void allGatherVImpl(T * values, int * nb_values) const; template void reduceImpl(T * values, int nb_values, const SynchronizerOperation & op, int root = 0) const; template void gatherImpl(T * values, int nb_values, int root = 0) const; template void gatherImpl(T * values, int nb_values, T * gathered, int nb_gathered = 0) const; template void gatherVImpl(T * values, int * nb_values, int root = 0) const; template void broadcastImpl(T * values, int nb_values, int root = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: Int getNbProc() const { return psize; }; Int whoAmI() const { return prank; }; static Communicator & getStaticCommunicator(); static Communicator & getStaticCommunicator(int & argc, char **& argv); int getMaxTag() const; int getMinTag() const; AKANTU_GET_MACRO(CommunicatorData, (*communicator_data), decltype(auto)); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: static std::unique_ptr static_communicator; protected: Int prank{0}; Int psize{1}; std::unique_ptr communicator_data; }; inline std::ostream & operator<<(std::ostream & stream, const CommunicationRequest & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "communicator_inline_impl.hh" #endif /* __AKANTU_STATIC_COMMUNICATOR_HH__ */ diff --git a/src/synchronizer/element_info_per_processor_tmpl.hh b/src/synchronizer/element_info_per_processor_tmpl.hh index 3fb055721..87dfcb647 100644 --- a/src/synchronizer/element_info_per_processor_tmpl.hh +++ b/src/synchronizer/element_info_per_processor_tmpl.hh @@ -1,147 +1,147 @@ /** * @file element_info_per_processor_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Tue Feb 20 2018 * * @brief Helper classes to create the distributed synchronizer and distribute * a mesh * * @section LICENSE * * 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 "element_group.hh" #include "element_info_per_processor.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_INFO_PER_PROCESSOR_TMPL_HH__ #define __AKANTU_ELEMENT_INFO_PER_PROCESSOR_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template void ElementInfoPerProc::fillMeshDataTemplated(BufferType & buffer, const std::string & tag_name, UInt nb_component) { AKANTU_DEBUG_ASSERT(this->mesh.getNbElement(this->type) == nb_local_element, "Did not got enought informations for the tag " << tag_name << " and the element type " << this->type << ":" << "_not_ghost." << " Got " << nb_local_element << " values, expected " << mesh.getNbElement(this->type)); - MeshData & mesh_data = this->getMeshData(); - mesh_data.registerElementalData(tag_name); - Array & data = mesh_data.getElementalDataArrayAlloc( + + mesh.registerElementalData(tag_name); + Array & data = mesh.getElementalDataArrayAlloc( tag_name, this->type, _not_ghost, nb_component); data.resize(nb_local_element); /// unpacking the data, element by element for (UInt i(0); i < nb_local_element; ++i) { for (UInt j(0); j < nb_component; ++j) { buffer >> data(i, j); } } AKANTU_DEBUG_ASSERT(mesh.getNbElement(this->type, _ghost) == nb_ghost_element, "Did not got enought informations for the tag " << tag_name << " and the element type " << this->type << ":" << "_ghost." << " Got " << nb_ghost_element << " values, expected " << mesh.getNbElement(this->type, _ghost)); - Array & data_ghost = mesh_data.getElementalDataArrayAlloc( + Array & data_ghost = mesh.getElementalDataArrayAlloc( tag_name, this->type, _ghost, nb_component); data_ghost.resize(nb_ghost_element); /// unpacking the ghost data, element by element for (UInt j(0); j < nb_ghost_element; ++j) { for (UInt k(0); k < nb_component; ++k) { buffer >> data_ghost(j, k); } } } /* -------------------------------------------------------------------------- */ template void ElementInfoPerProc::fillMeshData(BufferType & buffer, const std::string & tag_name, const MeshDataTypeCode & type_code, UInt nb_component) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ fillMeshDataTemplated(buffer, tag_name, \ nb_component); \ break; \ } switch (type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not determine the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ template void ElementInfoPerProc::fillElementGroupsFromBuffer( CommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); Element el; el.type = type; for (auto ghost_type : ghost_types) { UInt nb_element = mesh.getNbElement(type, ghost_type); el.ghost_type = ghost_type; for (UInt e = 0; e < nb_element; ++e) { el.element = e; std::vector element_to_group; buffer >> element_to_group; AKANTU_DEBUG_ASSERT(e < mesh.getNbElement(type, ghost_type), "The mesh does not have the element " << e); for (auto && element : element_to_group) { mesh.getElementGroup(element).add(el, false, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // akantu #endif /* __AKANTU_ELEMENT_INFO_PER_PROCESSOR_TMPL_HH__ */ diff --git a/src/synchronizer/master_element_info_per_processor.cc b/src/synchronizer/master_element_info_per_processor.cc index 74053d36e..27514c1de 100644 --- a/src/synchronizer/master_element_info_per_processor.cc +++ b/src/synchronizer/master_element_info_per_processor.cc @@ -1,463 +1,457 @@ /** * @file master_element_info_per_processor.cc * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Tue Feb 20 2018 * * @brief Helper class to distribute a mesh * * @section LICENSE * * 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 "aka_iterators.hh" #include "communicator.hh" #include "element_group.hh" #include "element_info_per_processor.hh" #include "element_synchronizer.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MasterElementInfoPerProc::MasterElementInfoPerProc( ElementSynchronizer & synchronizer, UInt message_cnt, UInt root, ElementType type, const MeshPartition & partition) : ElementInfoPerProc(synchronizer, message_cnt, root, type), partition(partition), all_nb_local_element(nb_proc, 0), all_nb_ghost_element(nb_proc, 0), all_nb_element_to_send(nb_proc, 0) { Vector size(5); size(0) = (UInt)type; if (type != _not_defined) { nb_nodes_per_element = Mesh::getNbNodesPerElement(type); nb_element = mesh.getNbElement(type); const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); for (UInt el = 0; el < nb_element; ++el) { this->all_nb_local_element[partition_num(el)]++; for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { this->all_nb_ghost_element[*part]++; } this->all_nb_element_to_send[partition_num(el)] += ghost_partition.getNbCols(el) + 1; } /// tag info - std::vector tag_names; - this->getMeshData().getTagNames(tag_names, type); + auto && tag_names = this->mesh.getTagNames(type); this->nb_tags = tag_names.size(); size(4) = nb_tags; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { size(1) = this->all_nb_local_element[p]; size(2) = this->all_nb_ghost_element[p]; size(3) = this->all_nb_element_to_send[p]; AKANTU_DEBUG_INFO( "Sending connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_SIZES) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_SIZES)); } else { this->nb_local_element = this->all_nb_local_element[p]; this->nb_ghost_element = this->all_nb_ghost_element[p]; } } } else { for (UInt p = 0; p < this->nb_proc; ++p) { if (p != this->root) { AKANTU_DEBUG_INFO( "Sending empty connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_SIZES) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_SIZES)); } } } } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizeConnectivities() { const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector> buffers(this->nb_proc); const auto & connectivities = this->mesh.getConnectivity(this->type, _not_ghost); /// copying the local connectivity for (auto && part_conn : zip(partition_num, make_view(connectivities, this->nb_nodes_per_element))) { auto && part = std::get<0>(part_conn); auto && conn = std::get<1>(part_conn); for (UInt i = 0; i < conn.size(); ++i) { buffers[part].push_back(conn[i]); } } /// copying the connectivity of ghost element for (auto && tuple : enumerate(make_view(connectivities, this->nb_nodes_per_element))) { auto && el = std::get<0>(tuple); auto && conn = std::get<1>(tuple); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { UInt proc = *part; for (UInt i = 0; i < conn.size(); ++i) { buffers[proc].push_back(conn[i]); } } } #ifndef AKANTU_NDEBUG for (auto p : arange(this->nb_proc)) { UInt size = this->nb_nodes_per_element * (this->all_nb_local_element[p] + this->all_nb_ghost_element[p]); AKANTU_DEBUG_ASSERT( buffers[p].size() == size, "The connectivity data packed in the buffer are not correct"); } #endif /// send all connectivity and ghost information to all processors std::vector requests; for (auto p : arange(this->nb_proc)) { if (p == this->root) continue; auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_CONNECTIVITY); AKANTU_DEBUG_INFO("Sending connectivities to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } Array & old_nodes = this->getNodesGlobalIds(); /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, buffers[root], all_nb_local_element[root], all_nb_ghost_element[root], type, old_nodes); comm.waitAll(requests); comm.freeCommunicationRequest(requests); } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizePartitions() { const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector> buffers(this->partition.getNbPartition()); /// splitting the partition information to send them to processors Vector count_by_proc(nb_proc, 0); for (UInt el = 0; el < nb_element; ++el) { UInt proc = partition_num(el); buffers[proc].push_back(ghost_partition.getNbCols(el)); UInt i(0); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[proc].push_back(*part); } } for (UInt el = 0; el < nb_element; ++el) { UInt i(0); for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[*part].push_back(partition_num(el)); } } #ifndef AKANTU_NDEBUG for (UInt p = 0; p < this->nb_proc; ++p) { AKANTU_DEBUG_ASSERT(buffers[p].size() == (this->all_nb_ghost_element[p] + this->all_nb_element_to_send[p]), "Data stored in the buffer are most probably wrong"); } #endif std::vector requests; /// last data to compute the communication scheme for (UInt p = 0; p < this->nb_proc; ++p) { if (p == this->root) continue; auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_PARTITIONS); AKANTU_DEBUG_INFO("Sending partition informations to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } if (Mesh::getSpatialDimension(this->type) == this->mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); this->fillCommunicationScheme(buffers[this->rank]); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeTags() { AKANTU_DEBUG_IN(); if (this->nb_tags == 0) { AKANTU_DEBUG_OUT(); return; } UInt mesh_data_sizes_buffer_length; - auto & mesh_data = this->getMeshData(); /// tag info - std::vector tag_names; - mesh_data.getTagNames(tag_names, type); + auto tag_names = mesh.getTagNames(type); + // Make sure the tags are sorted (or at least not in random order), // because they come from a map !! std::sort(tag_names.begin(), tag_names.end()); // Sending information about the tags in mesh_data: name, data type and // number of components of the underlying array associated to the current // type DynamicCommunicationBuffer mesh_data_sizes_buffer; for (auto && tag_name : tag_names) { mesh_data_sizes_buffer << tag_name; - mesh_data_sizes_buffer << mesh_data.getTypeCode(tag_name); - mesh_data_sizes_buffer << mesh_data.getNbComponent(tag_name, type); + mesh_data_sizes_buffer << mesh.getTypeCode(tag_name); + mesh_data_sizes_buffer << mesh.getNbComponent(tag_name, type); } mesh_data_sizes_buffer_length = mesh_data_sizes_buffer.size(); AKANTU_DEBUG_INFO( "Broadcasting the size of the information about the mesh data tags: (" << mesh_data_sizes_buffer_length << ")."); comm.broadcast(mesh_data_sizes_buffer_length, root); AKANTU_DEBUG_INFO( "Broadcasting the information about the mesh data tags, addr " << (void *)mesh_data_sizes_buffer.storage()); if (mesh_data_sizes_buffer_length != 0) comm.broadcast(mesh_data_sizes_buffer, root); if (mesh_data_sizes_buffer_length != 0) { // Sending the actual data to each processor std::vector buffers(nb_proc); // Loop over each tag for the current type for (auto && tag_name : tag_names) { // Type code of the current tag (i.e. the tag named *names_it) this->fillTagBuffer(buffers, tag_name); } std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == root) continue; auto && tag = Tag::genTag(this->rank, this->message_count, Tag::_MESH_DATA); AKANTU_DEBUG_INFO("Sending " << buffers[p].size() << " bytes of mesh data to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } // Loop over each tag for the current type for (auto && tag_name : tag_names) { // Reinitializing the mesh data on the master this->fillMeshData(buffers[root], tag_name, - mesh_data.getTypeCode(tag_name), - mesh_data.getNbComponent(tag_name, type)); + mesh.getTypeCode(tag_name), + mesh.getNbComponent(tag_name, type)); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MasterElementInfoPerProc::fillTagBufferTemplated( std::vector & buffers, const std::string & tag_name) { - MeshData & mesh_data = this->getMeshData(); - - const auto & data = mesh_data.getElementalDataArray(tag_name, type); + const auto & data = mesh.getElementalDataArray(tag_name, type); const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); // Not possible to use the iterator because it potentially triggers the // creation of complex // type templates (such as akantu::Vector< std::vector > which don't // implement the right interface // (e.g. operator<< in that case). // typename Array::template const_iterator< Vector > data_it = // data.begin(data.getNbComponent()); // typename Array::template const_iterator< Vector > data_end = // data.end(data.getNbComponent()); const T * data_it = data.storage(); const T * data_end = data.storage() + data.size() * data.getNbComponent(); const UInt * part = partition_num.storage(); /// copying the data, element by element for (; data_it != data_end; ++part) { for (UInt j(0); j < data.getNbComponent(); ++j, ++data_it) { buffers[*part] << *data_it; } } data_it = data.storage(); /// copying the data for the ghost element for (UInt el(0); data_it != data_end; data_it += data.getNbComponent(), ++el) { auto it = ghost_partition.begin(el); auto end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; for (UInt j(0); j < data.getNbComponent(); ++j) { buffers[proc] << data_it[j]; } } } } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::fillTagBuffer( std::vector & buffers, const std::string & tag_name) { - MeshData & mesh_data = this->getMeshData(); - #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ this->fillTagBufferTemplated(buffers, \ tag_name); \ break; \ } - MeshDataTypeCode data_type_code = mesh_data.getTypeCode(tag_name); + MeshDataTypeCode data_type_code = mesh.getTypeCode(tag_name); switch (data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); std::vector buffers(nb_proc); using ElementToGroup = std::vector>; ElementToGroup element_to_group(nb_element); for (auto & eg : ElementGroupsIterable(mesh)) { const auto & name = eg.getName(); for (const auto & element : eg.getElements(type, _not_ghost)) { element_to_group[element].push_back(name); } auto eit = eg.begin(type, _not_ghost); if (eit != eg.end(type, _not_ghost)) const_cast &>(eg.getElements(type)).empty(); } const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); /// copying the data, element by element for (auto && pair : zip(partition_num, element_to_group)) { buffers[std::get<0>(pair)] << std::get<1>(pair); } /// copying the data for the ghost element for (auto && pair : enumerate(element_to_group)) { auto && el = std::get<0>(pair); auto it = ghost_partition.begin(el); auto end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; buffers[proc] << std::get<1>(pair); } } std::vector requests; for (UInt p = 0; p < this->nb_proc; ++p) { if (p == this->rank) continue; auto && tag = Tag::genTag(this->rank, p, Tag::_ELEMENT_GROUP); AKANTU_DEBUG_INFO("Sending element groups to proc " << p << " TAG(" << tag << ")"); requests.push_back(comm.asyncSend(buffers[p], p, tag)); } this->fillElementGroupsFromBuffer(buffers[this->rank]); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/node_info_per_processor.cc b/src/synchronizer/node_info_per_processor.cc index a450c7761..d5e222476 100644 --- a/src/synchronizer/node_info_per_processor.cc +++ b/src/synchronizer/node_info_per_processor.cc @@ -1,499 +1,632 @@ /** * @file node_info_per_processor.cc * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Wed Nov 08 2017 * * @brief Please type the brief for file: Helper classes to create the * distributed synchronizer and distribute a mesh * * @section LICENSE * * 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 "node_info_per_processor.hh" #include "communicator.hh" #include "node_group.hh" #include "node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NodeInfoPerProc::NodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root) : MeshAccessor(synchronizer.getMesh()), synchronizer(synchronizer), comm(synchronizer.getCommunicator()), rank(comm.whoAmI()), nb_proc(comm.getNbProc()), root(root), mesh(synchronizer.getMesh()), spatial_dimension(synchronizer.getMesh().getSpatialDimension()), message_count(message_cnt) {} /* -------------------------------------------------------------------------- */ template void NodeInfoPerProc::fillNodeGroupsFromBuffer(CommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); std::vector> node_to_group; buffer >> node_to_group; AKANTU_DEBUG_ASSERT(node_to_group.size() == mesh.getNbGlobalNodes(), "Not the good amount of nodes where transmitted"); const auto & global_nodes = mesh.getGlobalNodesIds(); - auto nbegin = global_nodes.begin(); - auto nit = global_nodes.begin(); - auto nend = global_nodes.end(); - for (; nit != nend; ++nit) { - auto it = node_to_group[*nit].begin(); - auto end = node_to_group[*nit].end(); - - for (; it != end; ++it) { - mesh.getNodeGroup(*it).add(nit - nbegin, false); + for (auto && data : enumerate(global_nodes)) { + for (const auto & node : node_to_group[std::get<1>(data)]) { + mesh.getNodeGroup(node).add(std::get<0>(data), false); } } - GroupManager::const_node_group_iterator ngi = mesh.node_group_begin(); - GroupManager::const_node_group_iterator nge = mesh.node_group_end(); - for (; ngi != nge; ++ngi) { - NodeGroup & ng = *(ngi->second); - ng.optimize(); + for (auto && ng_data : mesh.iterateNodeGroups()) { + ng_data.optimize(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::fillNodesType() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); Array & nodes_type = this->getNodesType(); Array nodes_set(nb_nodes); nodes_set.set(0); enum NodeSet { NORMAL_SET = 1, GHOST_SET = 2, }; Array already_seen(nb_nodes, 1, false); for (UInt g = _not_ghost; g <= _ghost; ++g) { auto gt = (GhostType)g; UInt set = NORMAL_SET; if (gt == _ghost) set = GHOST_SET; already_seen.set(false); Mesh::type_iterator it = mesh.firstType(_all_dimensions, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type, gt); Array::const_vector_iterator conn_it = mesh.getConnectivity(type, gt).begin(nb_nodes_per_element); for (UInt e = 0; e < nb_element; ++e, ++conn_it) { const Vector & conn = *conn_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { AKANTU_DEBUG_ASSERT(conn(n) < nb_nodes, "Node " << conn(n) << " bigger than number of nodes " << nb_nodes); if (!already_seen(conn(n))) { nodes_set(conn(n)) += set; already_seen(conn(n)) = true; } } } } } for (UInt i = 0; i < nb_nodes; ++i) { if (nodes_set(i) == NORMAL_SET) nodes_type(i) = _nt_normal; else if (nodes_set(i) == GHOST_SET) nodes_type(i) = _nt_pure_ghost; else if (nodes_set(i) == (GHOST_SET + NORMAL_SET)) nodes_type(i) = _nt_master; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::fillCommunicationScheme(const Array & master_info) { AKANTU_DEBUG_IN(); Communications & communications = this->synchronizer.getCommunications(); { // send schemes auto it = master_info.begin_reinterpret(2, master_info.size() / 2); auto end = master_info.end_reinterpret(2, master_info.size() / 2); std::map> send_array_per_proc; for (; it != end; ++it) { const Vector & send_info = *it; send_array_per_proc[send_info(0)].push_back(send_info(1)); } for (auto & send_schemes : send_array_per_proc) { auto & scheme = communications.createSendScheme(send_schemes.first); auto & sends = send_schemes.second; std::sort(sends.begin(), sends.end()); std::transform(sends.begin(), sends.end(), sends.begin(), [this](UInt g) -> UInt { return mesh.getNodeLocalId(g); }); scheme.copy(sends); } } { // receive schemes const Array & nodes_type = this->getNodesType(); std::map> recv_array_per_proc; UInt node = 0; for (auto & node_type : nodes_type) { if (Int(node_type) >= 0) { recv_array_per_proc[node_type].push_back(mesh.getNodeGlobalId(node)); } ++node; } for (auto & recv_schemes : recv_array_per_proc) { auto & scheme = communications.createRecvScheme(recv_schemes.first); auto & recvs = recv_schemes.second; std::sort(recvs.begin(), recvs.end()); std::transform(recvs.begin(), recvs.end(), recvs.begin(), [this](UInt g) -> UInt { return mesh.getNodeLocalId(g); }); scheme.copy(recvs); } } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void NodeInfoPerProc::fillNodalData(DynamicCommunicationBuffer & buffer, + std::string tag_name) { + +#define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, _, elem) \ + case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ + auto & nodal_data = \ + mesh.getNodalData(tag_name); \ + nodal_data.resize(mesh.getNbNodes()); \ + for (auto && data : make_view(nodal_data)) { \ + for (auto i[[gnu::unused]] : arange(nodal_data.getNbComponent())) { \ + buffer >> data; \ + } \ + } \ + break; \ + } + + MeshDataTypeCode data_type_code = + mesh.getTypeCode(tag_name, MeshDataType::_nodal); + switch (data_type_code) { + BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , + AKANTU_MESH_DATA_TYPES) + default: + AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!"); + break; + } +#undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA +} + /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ MasterNodeInfoPerProc::MasterNodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root) : NodeInfoPerProc(synchronizer, message_cnt, root) { UInt nb_global_nodes = this->mesh.getNbGlobalNodes(); this->comm.broadcast(nb_global_nodes, this->root); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizeNodes() { this->nodes_per_proc.resize(nb_proc); this->nb_nodes_per_proc.resize(nb_proc); Array local_nodes(0, spatial_dimension); Array & nodes = this->getNodes(); for (UInt p = 0; p < nb_proc; ++p) { UInt nb_nodes = 0; // UInt * buffer; Array * nodes_to_send; Array & nodespp = nodes_per_proc[p]; if (p != root) { nodes_to_send = new Array(0, spatial_dimension); AKANTU_DEBUG_INFO("Receiving number of nodes from proc " << p << " " << Tag::genTag(p, 0, Tag::_NB_NODES)); comm.receive(nb_nodes, p, Tag::genTag(p, 0, Tag::_NB_NODES)); nodespp.resize(nb_nodes); this->nb_nodes_per_proc(p) = nb_nodes; AKANTU_DEBUG_INFO("Receiving list of nodes from proc " << p << " " << Tag::genTag(p, 0, Tag::_NODES)); comm.receive(nodespp, p, Tag::genTag(p, 0, Tag::_NODES)); } else { Array & local_ids = this->getNodesGlobalIds(); this->nb_nodes_per_proc(p) = local_ids.size(); nodespp.copy(local_ids); nodes_to_send = &local_nodes; } Array::const_scalar_iterator it = nodespp.begin(); Array::const_scalar_iterator end = nodespp.end(); /// get the coordinates for the selected nodes for (; it != end; ++it) { Vector coord(nodes.storage() + spatial_dimension * *it, spatial_dimension); nodes_to_send->push_back(coord); } if (p != root) { /// send them for distant processors AKANTU_DEBUG_INFO("Sending coordinates to proc " << p << " " << Tag::genTag(this->rank, 0, Tag::_COORDINATES)); comm.send(*nodes_to_send, p, Tag::genTag(this->rank, 0, Tag::_COORDINATES)); delete nodes_to_send; } } /// construct the local nodes coordinates nodes.copy(local_nodes); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizeTypes() { // > std::multimap> nodes_to_proc; std::vector> nodes_type_per_proc(nb_proc); // arrays containing pairs of (proc, node) std::vector> nodes_to_send_per_proc(nb_proc); for (UInt p = 0; p < nb_proc; ++p) { nodes_type_per_proc[p].resize(nb_nodes_per_proc(p)); } this->fillNodesType(); for (UInt p = 0; p < nb_proc; ++p) { auto & nodes_types = nodes_type_per_proc[p]; if (p != root) { AKANTU_DEBUG_INFO( "Receiving first nodes types from proc " << p << " " << Tag::genTag(this->rank, this->message_count, Tag::_NODES_TYPE)); comm.receive(nodes_types, p, Tag::genTag(p, 0, Tag::_NODES_TYPE)); } else { nodes_types.copy(this->getNodesType()); } // stack all processors claiming to be master for a node for (UInt local_node = 0; local_node < nb_nodes_per_proc(p); ++local_node) { if (nodes_types(local_node) == _nt_master) { UInt global_node = nodes_per_proc[p](local_node); nodes_to_proc.insert( std::make_pair(global_node, std::make_pair(p, local_node))); } } } for (UInt i = 0; i < mesh.getNbGlobalNodes(); ++i) { auto it_range = nodes_to_proc.equal_range(i); if (it_range.first == nodes_to_proc.end() || it_range.first->first != i) continue; // pick the first processor out of the multi-map as the actual master UInt master_proc = (it_range.first)->second.first; for (auto it_node = it_range.first; it_node != it_range.second; ++it_node) { UInt proc = it_node->second.first; UInt node = it_node->second.second; if (proc != master_proc) { // store the info on all the slaves for a given master nodes_type_per_proc[proc](node) = NodeType(master_proc); nodes_to_send_per_proc[master_proc].push_back(proc); nodes_to_send_per_proc[master_proc].push_back(i); } } } std::vector requests_send_type; std::vector requests_send_master_info; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { AKANTU_DEBUG_INFO("Sending nodes types to proc " << p << " " << Tag::genTag(this->rank, 0, Tag::_NODES_TYPE)); requests_send_type.push_back( comm.asyncSend(nodes_type_per_proc[p], p, Tag::genTag(this->rank, 0, Tag::_NODES_TYPE))); auto & nodes_to_send = nodes_to_send_per_proc[p]; /// push back an element to avoid a send of size 0 nodes_to_send.push_back(-1); AKANTU_DEBUG_INFO("Sending nodes master info to proc " << p << " " << Tag::genTag(this->rank, 1, Tag::_NODES_TYPE)); requests_send_master_info.push_back(comm.asyncSend( nodes_to_send, p, Tag::genTag(this->rank, 1, Tag::_NODES_TYPE))); } else { this->getNodesType().copy(nodes_type_per_proc[p]); this->fillCommunicationScheme(nodes_to_send_per_proc[root]); } } comm.waitAll(requests_send_type); comm.freeCommunicationRequest(requests_send_type); requests_send_type.clear(); comm.waitAll(requests_send_master_info); comm.freeCommunicationRequest(requests_send_master_info); } /* -------------------------------------------------------------------------- */ void MasterNodeInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); UInt nb_total_nodes = mesh.getNbGlobalNodes(); DynamicCommunicationBuffer buffer; using NodeToGroup = std::vector>; NodeToGroup node_to_group; node_to_group.resize(nb_total_nodes); GroupManager::const_node_group_iterator ngi = mesh.node_group_begin(); GroupManager::const_node_group_iterator nge = mesh.node_group_end(); for (; ngi != nge; ++ngi) { NodeGroup & ng = *(ngi->second); std::string name = ngi->first; NodeGroup::const_node_iterator nit = ng.begin(); NodeGroup::const_node_iterator nend = ng.end(); for (; nit != nend; ++nit) { node_to_group[*nit].push_back(name); } nit = ng.begin(); if (nit != nend) ng.empty(); } buffer << node_to_group; std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == this->rank) continue; AKANTU_DEBUG_INFO("Sending node groups to proc " << p << " " << Tag::genTag(this->rank, p, Tag::_NODE_GROUP)); requests.push_back(comm.asyncSend( buffer, p, Tag::genTag(this->rank, p, Tag::_NODE_GROUP))); } this->fillNodeGroupsFromBuffer(buffer); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void MasterNodeInfoPerProc::fillTagBuffers( + std::vector & buffers, + const std::string & tag_name) { +#define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, _, elem) \ + case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ + auto & nodal_data = \ + mesh.getNodalData(tag_name); \ + for (auto && data : enumerate(nodes_per_proc)) { \ + auto proc = std::get<0>(data); \ + auto & nodes = std::get<1>(data); \ + auto & buffer = buffers[proc]; \ + for (auto & node : nodes) { \ + for (auto i : arange(nodal_data.getNbComponent())) { \ + buffer << nodal_data(node, i); \ + } \ + } \ + } \ + break; \ + } + + MeshDataTypeCode data_type_code = + mesh.getTypeCode(tag_name, MeshDataType::_nodal); + switch (data_type_code) { + BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , + AKANTU_MESH_DATA_TYPES) + default: + AKANTU_ERROR("Could not obtain the type of tag" << tag_name << "!"); + break; + } +#undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA +} // namespace akantu + +/* -------------------------------------------------------------------------- */ +void MasterNodeInfoPerProc::synchronizeTags() { + /// tag info + auto tag_names = mesh.getTagNames(); + + DynamicCommunicationBuffer tags_buffer; + for (auto && tag_name : tag_names) { + tags_buffer << tag_name; + tags_buffer << mesh.getTypeCode(tag_name, MeshDataType::_nodal); + tags_buffer << mesh.getNbComponent(tag_name); + } + + UInt tags_buffer_length = tags_buffer.size(); + AKANTU_DEBUG_INFO( + "Broadcasting the size of the information about the mesh data tags: (" + << tags_buffer_length << ")."); + comm.broadcast(tags_buffer_length, root); + + if (tags_buffer_length != 0) + comm.broadcast(tags_buffer, root); + + for (auto && tag_data : enumerate(tag_names)) { + auto tag_count = std::get<0>(tag_data); + auto & tag_name = std::get<1>(tag_data); + std::vector buffers; + std::vector requests; + buffers.resize(nb_proc); + fillTagBuffers(buffers, tag_name); + for (auto && data : enumerate(buffers)) { + auto && proc = std::get<0>(data); + auto & buffer = std::get<1>(data); + if (proc == root) { + fillNodalData(buffer, tag_name); + } else { + auto && tag = Tag::genTag(this->rank, tag_count, Tag::_MESH_DATA); + requests.push_back(comm.asyncSend(buffer, proc, tag)); + } + } + + comm.waitAll(requests); + comm.freeCommunicationRequest(requests); + } +} + /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ SlaveNodeInfoPerProc::SlaveNodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root) : NodeInfoPerProc(synchronizer, message_cnt, root) { UInt nb_global_nodes = 0; comm.broadcast(nb_global_nodes, root); this->setNbGlobalNodes(nb_global_nodes); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizeNodes() { AKANTU_DEBUG_INFO("Sending list of nodes to proc " << root << " " << Tag::genTag(this->rank, 0, Tag::_NB_NODES) << " " << Tag::genTag(this->rank, 0, Tag::_NODES)); Array & local_ids = this->getNodesGlobalIds(); Array & nodes = this->getNodes(); UInt nb_nodes = local_ids.size(); comm.send(nb_nodes, root, Tag::genTag(this->rank, 0, Tag::_NB_NODES)); comm.send(local_ids, root, Tag::genTag(this->rank, 0, Tag::_NODES)); /* --------<<<<-COORDINATES---------------------------------------------- */ nodes.resize(nb_nodes); AKANTU_DEBUG_INFO("Receiving coordinates from proc " << root << " " << Tag::genTag(root, 0, Tag::_COORDINATES)); comm.receive(nodes, root, Tag::genTag(root, 0, Tag::_COORDINATES)); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizeTypes() { this->fillNodesType(); Array & nodes_types = this->getNodesType(); AKANTU_DEBUG_INFO("Sending first nodes types to proc " << root << "" << Tag::genTag(this->rank, 0, Tag::_NODES_TYPE)); comm.send(nodes_types, root, Tag::genTag(this->rank, 0, Tag::_NODES_TYPE)); AKANTU_DEBUG_INFO("Receiving nodes types from proc " << root << " " << Tag::genTag(root, 0, Tag::_NODES_TYPE)); comm.receive(nodes_types, root, Tag::genTag(root, 0, Tag::_NODES_TYPE)); AKANTU_DEBUG_INFO("Receiving nodes master info from proc " << root << " " << Tag::genTag(root, 1, Tag::_NODES_TYPE)); CommunicationStatus status; comm.probe(root, Tag::genTag(root, 1, Tag::_NODES_TYPE), status); Array nodes_master_info(status.size()); comm.receive(nodes_master_info, root, Tag::genTag(root, 1, Tag::_NODES_TYPE)); nodes_master_info.resize(nodes_master_info.size() - 1); this->fillCommunicationScheme(nodes_master_info); } /* -------------------------------------------------------------------------- */ void SlaveNodeInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Receiving node groups from proc " << root << " " << Tag::genTag(root, this->rank, Tag::_NODE_GROUP)); CommunicationStatus status; comm.probe(root, Tag::genTag(root, this->rank, Tag::_NODE_GROUP), status); CommunicationBuffer buffer(status.size()); comm.receive(buffer, root, Tag::genTag(root, this->rank, Tag::_NODE_GROUP)); this->fillNodeGroupsFromBuffer(buffer); AKANTU_DEBUG_OUT(); } -} // akantu +/* -------------------------------------------------------------------------- */ +void SlaveNodeInfoPerProc::synchronizeTags() { + UInt tags_buffer_length{0}; + comm.broadcast(tags_buffer_length, root); + if (tags_buffer_length == 0) { + return; + } + + CommunicationBuffer tags_buffer; + tags_buffer.resize(tags_buffer_length); + comm.broadcast(tags_buffer, root); + + std::vector tag_names; + while (tags_buffer.getLeftToUnpack() > 0) { + std::string name; + MeshDataTypeCode code; + UInt nb_components; + tags_buffer >> name; + tags_buffer >> code; + tags_buffer >> nb_components; + + mesh.registerNodalData(name, nb_components, code); + tag_names.push_back(name); + } + + for (auto && tag_data : enumerate(tag_names)) { + auto tag_count = std::get<0>(tag_data); + auto & tag_name = std::get<1>(tag_data); + DynamicCommunicationBuffer buffer; + auto && tag = Tag::genTag(this->root, tag_count, Tag::_MESH_DATA); + comm.receive(buffer, this->root, tag); + + fillNodalData(buffer, tag_name); + } +} + +} // namespace akantu diff --git a/src/synchronizer/node_info_per_processor.hh b/src/synchronizer/node_info_per_processor.hh index 1bb583524..cccded8a3 100644 --- a/src/synchronizer/node_info_per_processor.hh +++ b/src/synchronizer/node_info_per_processor.hh @@ -1,108 +1,116 @@ /** * @file node_info_per_processor.hh * * @author Nicolas Richart * * @date creation: Wed Mar 16 2016 * @date last modification: Wed Nov 08 2017 * * @brief Helper classes to create the distributed synchronizer and distribute * a mesh * * @section LICENSE * * 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 "communication_buffer.hh" #include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NODE_INFO_PER_PROCESSOR_HH__ #define __AKANTU_NODE_INFO_PER_PROCESSOR_HH__ namespace akantu { class NodeSynchronizer; class Communicator; -} +} // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { class NodeInfoPerProc : protected MeshAccessor { public: NodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root); virtual void synchronizeNodes() = 0; virtual void synchronizeTypes() = 0; virtual void synchronizeGroups() = 0; + virtual void synchronizeTags() = 0; protected: template void fillNodeGroupsFromBuffer(CommunicationBuffer & buffer); void fillNodesType(); void fillCommunicationScheme(const Array &); + void fillNodalData(DynamicCommunicationBuffer & buffer, std::string tag_name); protected: NodeSynchronizer & synchronizer; const Communicator & comm; UInt rank; UInt nb_proc; UInt root; Mesh & mesh; UInt spatial_dimension; UInt message_count; }; /* -------------------------------------------------------------------------- */ class MasterNodeInfoPerProc : protected NodeInfoPerProc { public: MasterNodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root); void synchronizeNodes() override; void synchronizeTypes() override; void synchronizeGroups() override; + void synchronizeTags() override; private: + void fillTagBuffers(std::vector & buffers, + const std::string & tag_name); + /// get the list of nodes to send and send them std::vector> nodes_per_proc; Array nb_nodes_per_proc; }; /* -------------------------------------------------------------------------- */ class SlaveNodeInfoPerProc : protected NodeInfoPerProc { public: SlaveNodeInfoPerProc(NodeSynchronizer & synchronizer, UInt message_cnt, UInt root); void synchronizeNodes() override; void synchronizeTypes() override; void synchronizeGroups() override; + void synchronizeTags() override; private: }; -} // akantu +} // namespace akantu #endif /* __AKANTU_NODE_INFO_PER_PROCESSOR_HH__ */ diff --git a/test/test_common/CMakeLists.txt b/test/test_common/CMakeLists.txt index f95146d65..3d283d25d 100644 --- a/test/test_common/CMakeLists.txt +++ b/test/test_common/CMakeLists.txt @@ -1,37 +1,36 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Fri Sep 03 2010 # @date last modification: Tue Dec 05 2017 # # @brief configurations for common tests # # @section LICENSE # # Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . # #=============================================================================== add_mesh(test_grid_mesh circle.geo 2 1) register_test(test_grid test_grid.cc DEPENDS test_grid_mesh PACKAGE core) #register_test(test_types test_types.cc PACKAGE core) register_gtest_sources(SOURCES test_csr.cc PACKAGE core) -register_gtest_sources(SOURCES test_arange_iterator.cc PACKAGE core HEADER_ONLY) -register_gtest_sources(SOURCES test_zip_iterator.cc PACKAGE core HEADER_ONLY) +register_gtest_sources(SOURCES test_iterators.cc PACKAGE core HEADER_ONLY) register_gtest_sources(SOURCES test_array.cc PACKAGE core) register_gtest_sources(SOURCES test_tensors.cc PACKAGE core) register_gtest_test(test_common) diff --git a/test/test_common/test_arange_iterator.cc b/test/test_common/test_arange_iterator.cc deleted file mode 100644 index eebc171da..000000000 --- a/test/test_common/test_arange_iterator.cc +++ /dev/null @@ -1,71 +0,0 @@ -/** - * @file test_arange_iterator.cc - * - * @author Nicolas Richart - * - * @date creation: Fri Jun 18 2010 - * @date last modification: Sun Dec 03 2017 - * - * @brief A Documented file. - * - * @section LICENSE - * - * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ - -/* -------------------------------------------------------------------------- */ -#include "aka_iterators.hh" -/* -------------------------------------------------------------------------- */ -#include -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -TEST(TestArangeIterator, Stop) { - size_t ref_i = 0; - for (auto i : arange(10)) { - EXPECT_EQ(ref_i, i); - ++ref_i; - } -} - -TEST(TestArangeIterator, StartStop) { - size_t ref_i = 1; - for (auto i : arange(1, 10)) { - EXPECT_EQ(ref_i, i); - ++ref_i; - } -} - -TEST(TestArangeIterator, StartStopStep) { - size_t ref_i = 1; - for (auto i : arange(1, 22, 2)) { - EXPECT_EQ(ref_i, i); - ref_i += 2; - } -} - -TEST(TestArangeIterator, StartStopStepZipped) { - int ref_i1 = -1, ref_i2 = 1; - for (auto && i : zip(arange(-1, -10, -1), arange(1, 18, 2))) { - EXPECT_EQ(ref_i1, std::get<0>(i)); - EXPECT_EQ(ref_i2, std::get<1>(i)); - ref_i1 += -1; - ref_i2 += 2; - } -} diff --git a/test/test_common/test_zip_iterator.cc b/test/test_common/test_iterators.cc similarity index 63% rename from test/test_common/test_zip_iterator.cc rename to test/test_common/test_iterators.cc index e8a39e0a7..b07cf52df 100644 --- a/test/test_common/test_zip_iterator.cc +++ b/test/test_common/test_iterators.cc @@ -1,181 +1,300 @@ /** * @file test_zip_iterator.cc * * @author Nicolas Richart * * @date creation: Fri Jul 21 2017 * @date last modification: Fri Dec 08 2017 * * @brief test the zip container and iterator * * @section LICENSE * * 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 "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; +/* -------------------------------------------------------------------------- */ template class A { public: A() = default; A(T a) : a(a){}; A(const A & other) : a(other.a), copy_counter(other.copy_counter + 1), move_counter(other.move_counter) {} A & operator=(const A & other) { if (this != &other) { a = other.a; copy_counter = other.copy_counter + 1; } return *this; } A(A && other) : a(std::move(other.a)), copy_counter(std::move(other.copy_counter)), move_counter(std::move(other.move_counter) + 1) {} A & operator=(A && other) { if (this != &other) { a = std::move(other.a); copy_counter = std::move(other.copy_counter); move_counter = std::move(other.move_counter) + 1; } return *this; } A & operator*=(const T & b) { a *= b; return *this; } T a; size_t copy_counter{0}; size_t move_counter{0}; }; template struct C { struct iterator { using reference = A; using difference_type = void; using iterator_category = std::input_iterator_tag; using value_type = A; using pointer = A *; iterator(T pos) : pos(std::move(pos)) {} A operator*() { return A(pos); } bool operator!=(const iterator & other) const { return pos != other.pos; } bool operator==(const iterator & other) const { return pos == other.pos; } iterator & operator++() { ++pos; return *this; } T pos; }; C(T begin_, T end_) : begin_(std::move(begin_)), end_(std::move(end_)) {} iterator begin() { return iterator(begin_); } iterator end() { return iterator(end_); } T begin_, end_; }; class TestZipFixutre : public ::testing::Test { protected: void SetUp() override { a.reserve(size); b.reserve(size); for (size_t i = 0; i < size; ++i) { a.emplace_back(i); b.emplace_back(i + size); } } template void check(A && a, B && b, size_t pos, size_t nb_copy, size_t nb_move) { EXPECT_EQ(pos, a.a); EXPECT_EQ(nb_copy, a.copy_counter); EXPECT_EQ(nb_move, a.move_counter); EXPECT_FLOAT_EQ(pos + this->size, b.a); EXPECT_EQ(nb_copy, b.copy_counter); EXPECT_EQ(nb_move, b.move_counter); } protected: size_t size{20}; std::vector> a; std::vector> b; }; TEST_F(TestZipFixutre, SimpleTest) { size_t i = 0; for (auto && pair : zip(this->a, this->b)) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 0); ++i; } } TEST_F(TestZipFixutre, ConstTest) { size_t i = 0; const auto & ca = this->a; const auto & cb = this->b; for (auto && pair : zip(ca, cb)) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 0); EXPECT_EQ(true, std::is_const< std::remove_reference_t(pair))>>::value); EXPECT_EQ(true, std::is_const< std::remove_reference_t(pair))>>::value); ++i; } } TEST_F(TestZipFixutre, MixteTest) { size_t i = 0; const auto & cb = this->b; for (auto && pair : zip(a, cb)) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 0); EXPECT_EQ(false, std::is_const< std::remove_reference_t(pair))>>::value); EXPECT_EQ(true, std::is_const< std::remove_reference_t(pair))>>::value); ++i; } } TEST_F(TestZipFixutre, MoveTest) { size_t i = 0; for (auto && pair : zip(C(0, this->size), C(this->size, 2 * this->size))) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 1); ++i; } } + +/* -------------------------------------------------------------------------- */ +TEST(TestArangeIterator, Stop) { + size_t ref_i = 0; + for (auto i : arange(10)) { + EXPECT_EQ(ref_i, i); + ++ref_i; + } +} + +TEST(TestArangeIterator, StartStop) { + size_t ref_i = 1; + for (auto i : arange(1, 10)) { + EXPECT_EQ(ref_i, i); + ++ref_i; + } +} + +TEST(TestArangeIterator, StartStopStep) { + size_t ref_i = 1; + for (auto i : arange(1, 22, 2)) { + EXPECT_EQ(ref_i, i); + ref_i += 2; + } +} + +TEST(TestArangeIterator, StartStopStepZipped) { + int ref_i1 = -1, ref_i2 = 1; + for (auto && i : zip(arange(-1, -10, -1), arange(1, 18, 2))) { + EXPECT_EQ(ref_i1, std::get<0>(i)); + EXPECT_EQ(ref_i2, std::get<1>(i)); + ref_i1 += -1; + ref_i2 += 2; + } +} + +/* -------------------------------------------------------------------------- */ +TEST(TestTransformAdaptor, Keys) { + std::map map{ + {"1", 1}, {"2", 2}, {"3", 3}, {"3", 3}, {"4", 4}}; + + char counter = '1'; + for (auto && key : make_keys_adaptor(map)) { + EXPECT_EQ(counter, key[0]); + ++counter; + } +} + +TEST(TestTransformAdaptor, Values) { + std::map map{ + {"1", 1}, {"2", 2}, {"3", 3}, {"3", 3}, {"4", 4}}; + + int counter = 1; + for (auto && value : make_values_adaptor(map)) { + EXPECT_EQ(counter, value); + ++counter; + } +} + +static int plus1(int value) { return value + 1; } + +struct Plus { + Plus(int a) : a(a) {} + int operator()(int b) { return a + b; } + +private: + int a{0}; +}; + +TEST(TestTransformAdaptor, Lambda) { + auto && container = arange(10); + + for (auto && data : + zip(container, make_transform_adaptor(container, [](auto && value) { + return value + 1; + }))) { + EXPECT_EQ(std::get<0>(data) + 1, std::get<1>(data)); + } +} + +TEST(TestTransformAdaptor, LambdaLambda) { + std::map map{ + {"1", 1}, {"2", 2}, {"3", 3}, {"3", 3}, {"4", 4}}; + + int counter = 1; + for (auto && data : make_transform_adaptor( + make_values_adaptor(map), [](auto && value) { return value + 1; })) { + EXPECT_EQ(counter + 1, data); + ++counter; + } + + auto && container = arange(10); + + for (auto && data : + zip(container, make_transform_adaptor(container, [](auto && value) { + return value + 1; + }))) { + EXPECT_EQ(std::get<0>(data) + 1, std::get<1>(data)); + } +} + +TEST(TestTransformAdaptor, Function) { + auto && container = arange(10); + + for (auto && data : + zip(container, make_transform_adaptor(container, plus1))) { + EXPECT_EQ(std::get<0>(data) + 1, std::get<1>(data)); + } +} + +TEST(TestTransformAdaptor, Functor) { + auto && container = arange(10); + + for (auto && data : + zip(container, make_transform_adaptor(container, Plus(1)))) { + EXPECT_EQ(std::get<0>(data) + 1, std::get<1>(data)); + } +} diff --git a/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc index 7b618d76c..1edead6c1 100644 --- a/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc +++ b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc @@ -1,105 +1,105 @@ /** * @file test_fe_engine_precomputation_bernoulli_3.cc * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Jan 24 2018 * * @brief test of the fem class * * @section LICENSE * * 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 "fe_engine.hh" #include "integrator_gauss.hh" #include "shape_structural.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; /** * Reference: p. 285, example 5.7 - A First Course in the Finite Elements Method * Logan, 6th Edition, 2016 * ISBN-13: 978-1-305-63734-4 */ Matrix rotationReference() { return {{3. / 13, 4. / 13, 12. / 13}, {-4. / 5, 3. / 5, 0}, {-36. / 65, -48. / 65, 5. / 13}}; } int main(int argc, char * argv[]) { akantu::initialize(argc, argv); // debug::setDebugLevel(dblTest); constexpr ElementType type = _bernoulli_beam_3; UInt dim = ElementClass::getSpatialDimension(); Mesh mesh(dim); // Pushing nodes Vector node = {0, 0, 0}; mesh.getNodes().push_back(node); node = {3, 4, 12}; mesh.getNodes().push_back(node); // Pushing connectivity mesh.addConnectivityType(type); auto & connectivity = mesh.getConnectivity(type); Vector elem = {0, 1}; connectivity.push_back(elem); // Pushing normals auto & normals = - mesh.registerData("extra_normal").alloc(0, dim, type, _not_ghost); + mesh.registerElementalData("extra_normal").alloc(0, dim, type, _not_ghost); Vector normal = {-36. / 65, -48. / 65, 5. / 13}; normals.push_back(normal); normals.push_back(normal); using FE = FEEngineTemplate; using ShapeStruct = ShapeStructural<_ek_structural>; auto fem = std::make_unique(mesh, dim, "test_fem"); fem->initShapeFunctions(); auto & shape = dynamic_cast(fem->getShapeFunctions()); Matrix rot_ref = rotationReference(); Matrix solution(6, 6); solution.block(rot_ref, 0, 0); solution.block(rot_ref, 3, 3); for (auto && rot : make_view(shape.getRotations(type), 6, 6)) { if (!Math::are_vector_equal(6 * 6, solution.storage(), rot.storage())) return 1; } /// TODO check shape functions and shape derivatives finalize(); return 0; } diff --git a/test/test_fe_engine/test_fe_engine_precomputation_structural.cc b/test/test_fe_engine/test_fe_engine_precomputation_structural.cc index 440c99af5..c461d68fb 100644 --- a/test/test_fe_engine/test_fe_engine_precomputation_structural.cc +++ b/test/test_fe_engine/test_fe_engine_precomputation_structural.cc @@ -1,126 +1,126 @@ /** * @file test_fe_engine_precomputation_structural.cc * * @author Lucas Frerot * * @date creation: Fri Jan 26 2018 * @date last modification: Mon Feb 19 2018 * * @brief test of the structural precomputations * * @section LICENSE * * 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 "fe_engine.hh" #include "integrator_gauss.hh" #include "shape_structural.hh" #include "test_fe_engine_structural_fixture.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ // Need a special fixture for the extra normal class TestBernoulliB3 : public TestFEMStructuralFixture> { using parent = TestFEMStructuralFixture>; public: /// Load the mesh and provide extra normal direction void readMesh(std::string filename) override { parent::readMesh(filename); - auto & normals = this->mesh->registerData("extra_normal") + auto & normals = this->mesh->registerElementalData("extra_normal") .alloc(0, dim, type, _not_ghost); Vector normal = {-36. / 65, -48. / 65, 5. / 13}; normals.push_back(normal); } }; /* -------------------------------------------------------------------------- */ /// Type alias using TestBernoulliB2 = TestFEMStructuralFixture>; using TestDKT18 = TestFEMStructuralFixture>; /* -------------------------------------------------------------------------- */ /// Solution for 2D rotation matrices Matrix globalToLocalRotation(Real theta) { auto c = std::cos(theta); auto s = std::sin(theta); return {{c, s, 0}, {-s, c, 0}, {0, 0, 1}}; } /* -------------------------------------------------------------------------- */ TEST_F(TestBernoulliB2, PrecomputeRotations) { this->fem->initShapeFunctions(); using ShapeStruct = ShapeStructural<_ek_structural>; auto & shape = dynamic_cast(fem->getShapeFunctions()); auto & rot = shape.getRotations(type); Real a = std::atan(4. / 3); std::vector angles = {a, -a, 0}; Math::setTolerance(1e-15); for (auto && tuple : zip(make_view(rot, ndof, ndof), angles)) { auto rotation = std::get<0>(tuple); auto angle = std::get<1>(tuple); auto rotation_error = (rotation - globalToLocalRotation(angle)).norm(); EXPECT_NEAR(rotation_error, 0., Math::getTolerance()); } } /* -------------------------------------------------------------------------- */ TEST_F(TestBernoulliB3, PrecomputeRotations) { this->fem->initShapeFunctions(); using ShapeStruct = ShapeStructural<_ek_structural>; auto & shape = dynamic_cast(fem->getShapeFunctions()); auto & rot = shape.getRotations(type); Matrix ref = {{3. / 13, 4. / 13, 12. / 13}, {-4. / 5, 3. / 5, 0}, {-36. / 65, -48. / 65, 5. / 13}}; Matrix solution{ndof, ndof}; solution.block(ref, 0, 0); solution.block(ref, dim, dim); // The default tolerance is too much, really Math::setTolerance(1e-15); for (auto & rotation : make_view(rot, ndof, ndof)) { auto rotation_error = (rotation - solution).norm(); EXPECT_NEAR(rotation_error, 0., Math::getTolerance()); } } /* -------------------------------------------------------------------------- */ TEST_F(TestDKT18, DISABLED_PrecomputeRotations) { this->fem->initShapeFunctions(); using ShapeStruct = ShapeStructural<_ek_structural>; auto & shape = dynamic_cast(fem->getShapeFunctions()); auto & rot = shape.getRotations(type); for (auto & rotation : make_view(rot, ndof, ndof)) { std::cout << rotation << "\n"; } std::cout.flush(); } diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc index d01d23ac7..f918249f7 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc @@ -1,99 +1,99 @@ /** * @file test_embedded_element_matrix.cc * * @author Lucas Frerot * * @date creation: Wed Mar 25 2015 * @date last modification: Fri Feb 09 2018 * * @brief test of the class EmbeddedInterfaceModel * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "embedded_interface_model.hh" #include "sparse_matrix_aij.hh" #include "sparse_solver.hh" using namespace akantu; int main(int argc, char * argv[]) { debug::setDebugLevel(dblWarning); initialize("embedded_element.dat", argc, argv); constexpr UInt dim = 2; constexpr ElementType type = _segment_2; const Real height = 0.4; Mesh mesh(dim); mesh.read("triangle.msh"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); auto & nodes = reinforcement_mesh.getNodes(); nodes.push_back(Vector({0, height})); nodes.push_back(Vector({1, height})); reinforcement_mesh.addConnectivityType(type); auto & connectivity = reinforcement_mesh.getConnectivity(type); connectivity.push_back(Vector({0, 1})); Array names_vec(1, 1, "reinforcement", "reinforcement_names"); - reinforcement_mesh.registerData("physical_names") + reinforcement_mesh.registerElementalData("physical_names") .alloc(1, 1, type); reinforcement_mesh.getData("physical_names")(type).copy( names_vec); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); model.initFull(_analysis_method = _static); if (model.getInterfaceMesh().getNbElement(type) != 1) return EXIT_FAILURE; if (model.getInterfaceMesh().getSpatialDimension() != 2) return EXIT_FAILURE; try { // matrix should be singular model.solveStep(); } catch (debug::SingularMatrixException & e) { std::cerr << "Matrix is singular, relax, everything is fine :)" << std::endl; } catch (debug::Exception & e) { std::cerr << "Unexpceted error: " << e.what() << std::endl; throw e; } SparseMatrixAIJ & K = dynamic_cast(model.getDOFManager().getMatrix("K")); K.saveMatrix("stiffness.mtx"); Math::setTolerance(1e-8); // Testing the assembled stiffness matrix if (!Math::are_float_equal(K(0, 0), 1. - height) || !Math::are_float_equal(K(0, 2), height - 1.) || !Math::are_float_equal(K(2, 0), height - 1.) || !Math::are_float_equal(K(2, 2), 1. - height)) return EXIT_FAILURE; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc index f0c1fad1f..5c4642fdf 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc @@ -1,108 +1,108 @@ /** * @file test_embedded_interface_model.cc * * @author Lucas Frerot * * @date creation: Wed Mar 25 2015 * @date last modification: Wed Jan 31 2018 * * @brief Embedded model test based on potential energy * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #include #include "aka_common.hh" #include "embedded_interface_model.hh" #include "sparse_matrix.hh" using namespace akantu; int main(int argc, char * argv[]) { debug::setDebugLevel(dblWarning); initialize("material.dat", argc, argv); UInt dim = 2; Math::setTolerance(1e-7); // Mesh here is a 1x1 patch Mesh mesh(dim); mesh.read("embedded_mesh.msh"); Array nodes_vec(2, dim, "reinforcement_nodes"); nodes_vec.storage()[0] = 0; nodes_vec.storage()[1] = 0.5; nodes_vec.storage()[2] = 1; nodes_vec.storage()[3] = 0.5; Array conn_vec(1, 2, "reinforcement_connectivity"); conn_vec.storage()[0] = 0; conn_vec.storage()[1] = 1; Array names_vec(1, 1, "reinforcement", "reinforcement_names"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); reinforcement_mesh.getNodes().copy(nodes_vec); reinforcement_mesh.addConnectivityType(_segment_2); reinforcement_mesh.getConnectivity(_segment_2).copy(conn_vec); - reinforcement_mesh.registerData("physical_names") + reinforcement_mesh.registerElementalData("physical_names") .alloc(1, 1, _segment_2); reinforcement_mesh.getData("physical_names")(_segment_2) .copy(names_vec); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); model.initFull(_analysis_method = _static); Array & nodes = mesh.getNodes(); Array & forces = model.getForce(); Array & bound = model.getBlockedDOFs(); forces(2, 0) = -250; forces(5, 0) = -500; forces(8, 0) = -250; for (UInt i = 0; i < mesh.getNbNodes(); i++) { if (Math::are_float_equal(nodes(i, 0), 0.)) bound(i, 0) = true; if (Math::are_float_equal(nodes(i, 1), 0.)) bound(i, 1) = true; } model.addDumpFieldVector("displacement"); model.addDumpFieldTensor("stress"); model.setBaseNameToDumper("reinforcement", "reinforcement"); model.addDumpFieldTensorToDumper("reinforcement", "stress_embedded"); model.solveStep(); model.getDOFManager().getMatrix("K").saveMatrix("matrix_test"); model.dump(); Real pot_energy = model.getEnergy("potential"); if (std::abs(pot_energy - 7.37343e-06) > 1e-5) return EXIT_FAILURE; finalize(); return 0; } diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc index 1fa50a1b7..dd92bd914 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc @@ -1,102 +1,102 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_3.cc * * @author Lucas Frerot * * @date creation: Fri Jul 15 2011 * @date last modification: Fri Feb 09 2018 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "test_structural_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; class TestStructBernoulli3 : public TestStructuralFixture> { using parent = TestStructuralFixture>; public: void readMesh(std::string filename) override { parent::readMesh(filename); auto & normals = - this->mesh->registerData("extra_normal") + this->mesh->registerElementalData("extra_normal") .alloc(0, parent::spatial_dimension, parent::type, _not_ghost); Vector normal = {0, 0, 1}; normals.push_back(normal); normal = {0, 0, 1}; normals.push_back(normal); } void addMaterials() override { StructuralMaterial mat; mat.E = 1; mat.Iz = 1; mat.Iy = 1; mat.A = 1; mat.GJ = 1; this->model->addMaterial(mat); } void setDirichlets() { // Boundary conditions (blocking all DOFs of nodes 2 & 3) auto boundary = ++this->model->getBlockedDOFs().begin(parent::ndof); // clang-format off *boundary = {true, true, true, true, true, true}; ++boundary; *boundary = {true, true, true, true, true, true}; ++boundary; // clang-format on } void setNeumanns() override { // Forces Real P = 1; // N auto & forces = this->model->getExternalForce(); forces.clear(); forces(0, 2) = -P; // vertical force on first node } void assignMaterials() override { model->getElementMaterial(parent::type).set(0); } }; /* -------------------------------------------------------------------------- */ TEST_F(TestStructBernoulli3, TestDisplacements) { model->solveStep(); auto vz = model->getDisplacement()(0, 2); auto thy = model->getDisplacement()(0, 4); auto thx = model->getDisplacement()(0, 3); auto thz = model->getDisplacement()(0, 5); Real tol = Math::getTolerance(); EXPECT_NEAR(vz, -5. / 48, tol); EXPECT_NEAR(thy, -std::sqrt(2) / 8, tol); EXPECT_NEAR(thz, 0, tol); EXPECT_NEAR(thx, 0, tol); } diff --git a/test/test_synchronizer/test_data_distribution.cc b/test/test_synchronizer/test_data_distribution.cc index 26e1fe79b..efe19ae09 100644 --- a/test/test_synchronizer/test_data_distribution.cc +++ b/test/test_synchronizer/test_data_distribution.cc @@ -1,74 +1,88 @@ /** * @file test_data_distribution.cc * * @author Nicolas Richart * * @date creation: Fri Sep 05 2014 * @date last modification: Fri Jan 26 2018 * * @brief Test the mesh distribution on creation of a distributed synchonizer * * @section LICENSE * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "test_synchronizers_fixture.hh" /* -------------------------------------------------------------------------- */ TEST_F(TestSynchronizerFixture, DataDistribution) { - auto & barycenters = this->mesh->registerData("barycenters"); + auto & barycenters = this->mesh->registerElementalData("barycenters"); auto spatial_dimension = this->mesh->getSpatialDimension(); barycenters.initialize(*this->mesh, _spatial_dimension = _all_dimensions, _nb_component = spatial_dimension); this->initBarycenters(barycenters, *this->mesh); + + auto & gids = this->mesh->registerNodalData("gid"); + gids.resize(this->mesh->getNbNodes()); + + for(auto && data : enumerate(gids)) { + std::get<1>(data) = std::get<0>(data); + } + this->distribute(); for (auto && ghost_type : ghost_types) { for (const auto & type : this->mesh->elementTypes(_all_dimensions, ghost_type)) { auto & barycenters = this->mesh->getData("barycenters", type, ghost_type); for (auto && data : enumerate(make_view(barycenters, spatial_dimension))) { Element element{type, UInt(std::get<0>(data)), ghost_type}; Vector barycenter(spatial_dimension); this->mesh->getBarycenter(element, barycenter); auto dist = (std::get<1>(data) - barycenter).template norm(); EXPECT_NEAR(dist, 0, 1e-7); } } } + + if (this->mesh->isDistributed()) { + for(auto && data : zip(gids, this->mesh->getGlobalNodesIds())) { + EXPECT_EQ(std::get<0>(data), std::get<1>(data)); + } + } } TEST_F(TestSynchronizerFixture, DataDistributionTags) { this->distribute(); for (const auto & type : this->mesh->elementTypes(_all_dimensions)) { auto & tags = this->mesh->getData("tag_0", type); Array::const_vector_iterator tags_it = tags.begin(1); Array::const_vector_iterator tags_end = tags.end(1); // The number of tags should match the number of elements on rank" EXPECT_EQ(this->mesh->getNbElement(type), tags.size()); } }