diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 07db769d9..36edc7915 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,526 +1,553 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ using MemoryID = UInt; // using Surface = std::string; // using SurfacePair= std::pair; // using SurfacePairList = std::list; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; #define AKANTU_PP_ENUM(s, data, i, elem) \ BOOST_PP_TUPLE_REM() \ elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data))) } // namespace akantu #if (defined(__GNUC__) || defined(__GNUG__)) #define AKA_GCC_VERSION \ (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) #if AKA_GCC_VERSION < 60000 #define AKANTU_ENUM_HASH(type_name) \ namespace std { \ template <> struct hash<::akantu::type_name> { \ using argument_type = ::akantu::type_name; \ size_t operator()(const argument_type & e) const noexcept { \ auto ue = underlying_type_t(e); \ return uh(ue); \ } \ \ private: \ const hash> uh{}; \ }; \ } #else #define AKANTU_ENUM_HASH(type_name) #endif // AKA_GCC_VERSION #endif // GNU #include "aka_element_classes_info.hh" namespace akantu { #define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) #define AKANTU_PP_TYPE_TO_STR(s, data, elem) \ ({BOOST_PP_CAT(data::_, elem), BOOST_PP_STRINGIZE(elem)}) #define AKANTU_PP_STR_TO_TYPE(s, data, elem) \ ({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data::_, elem)}) #define AKANTU_ENUM_DECLARE(type_name, list) \ enum class type_name { \ BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \ }; #define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ } \ AKANTU_ENUM_HASH(type_name) \ namespace aka { \ inline std::string to_string(const ::akantu::type_name & type) { \ static std::unordered_map<::akantu::type_name, std::string> convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, \ ::akantu::type_name, list))}; \ return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << aka::to_string(type); \ return stream; \ } #define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ inline std::istream & operator>>(std::istream & stream, type_name & type) { \ std::string str; \ stream >> str; \ static std::unordered_map convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, type_name, list))}; \ type = convert.at(str); \ return stream; \ } /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #ifndef SWIG // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { _model, _solid_mechanics_model, _solid_mechanics_model_cohesive, _heat_transfer_model, _structural_mechanics_model, _embedded_model }; #endif /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent _nls_lumped, ///< Case of lumped mass or equivalent matrix _nls_auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; -/// Define the node/dof type -enum NodeType : Int { - _np_periodic_ghost_slave = -5, - _np_periodic_slave = -4, - _nt_pure_gost = -3, - _nt_master = -2, - _nt_normal = -1, - _nt_slave = 0, -}; - /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_ask_nodes, _gst_size, //--- SolidMechanicsModel tags --- _gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _gst_smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _gst_smm_uv, ///< synchronization of the nodal velocities and displacement _gst_smm_res, ///< synchronization of the nodal residual _gst_smm_init_mat, ///< synchronization of the data to initialize materials _gst_smm_stress, ///< synchronization of the stresses to compute the internal /// forces _gst_smmc_facets, ///< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, ///< synchronization of facet global connectivity _gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet /// synch _gst_smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _gst_giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _gst_ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gst_gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _gst_htm_temperature, ///< synchronization of the nodal temperature _gst_htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _gst_htm_phi, ///< synchronization of the nodal level set value phi _gst_htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _gst_mnl_for_average, ///< synchronization of data to average in non local /// material _gst_mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _gst_nh_criterion, // --- General tags --- _gst_test, ///< Test tag _gst_user_1, ///< tag for user simulations _gst_user_2, ///< tag for user simulations _gst_material_id, ///< synchronization of the material ids _gst_for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _gst_cf_nodal, ///< synchronization of disp, velo, and current position _gst_cf_incr, ///< synchronization of increment // --- Solver tags --- _gst_solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; + + +/// Define the flag that can be set to a node +enum class NodeFlag : std::uint8_t { + _normal = 0x00, + _master = 0x01, + _slave = 0x03, + _pure_ghost = 0x07, + _shared_mask = 0x0F, + _local_master_mask = 0x0E, + _periodic = 0x10, + _periodic_slave = 0x10, + _periodic_master = 0x30, + _periodic_mask = 0xF0, +}; + +inline NodeFlag operator& (const NodeFlag & a, const NodeFlag & b) { + using under = std::underlying_type_t; + return NodeFlag(under(a) & under(b)); +} + +inline NodeFlag operator| (const NodeFlag & a, const NodeFlag & b) { + using under = std::underlying_type_t; + return NodeFlag(under(a) | under(b)); +} + +inline NodeFlag & operator|= (NodeFlag & a, const NodeFlag & b) { + a = a | b; + return a; } +inline std::ostream & operator<< (std::ostream & stream, const NodeFlag & flag) { + using under = std::underlying_type_t; + stream << under(flag); + return stream; +} + +} // namespace akantu + #ifndef SWIG AKANTU_ENUM_HASH(GhostType) #endif namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ struct TensorTrait {}; /* -------------------------------------------------------------------------- */ template using is_tensor = std::is_base_of; /* -------------------------------------------------------------------------- */ template using is_scalar = std::is_arithmetic; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ // #if defined(__INTEL_COMPILER) // /// remark #981: operands are evaluated in unspecified order // #pragma warning(disable : 981) // /// remark #383: value copied to temporary, reference to temporary used // #pragma warning(disable : 383) // #endif // defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ } // namespace akantu #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic number * is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/io/dumper/dumper_iohelper.cc b/src/io/dumper/dumper_iohelper.cc index d4c9443c0..c963934e1 100644 --- a/src/io/dumper/dumper_iohelper.cc +++ b/src/io/dumper/dumper_iohelper.cc @@ -1,317 +1,319 @@ /** * @file dumper_iohelper.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Oct 26 2012 * @date last modification: Tue Feb 20 2018 * * @brief implementation of DumperIOHelper * * @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 #include "dumper_elemental_field.hh" #include "dumper_filtered_connectivity.hh" #include "dumper_iohelper.hh" #include "dumper_nodal_field.hh" #include "dumper_variable.hh" #include "mesh.hh" #if defined(AKANTU_IGFEM) #include "dumper_igfem_connectivity.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DumperIOHelper::DumperIOHelper() = default; /* -------------------------------------------------------------------------- */ DumperIOHelper::~DumperIOHelper() { for (auto it = fields.begin(); it != fields.end(); ++it) { delete it->second; } delete dumper; } /* -------------------------------------------------------------------------- */ void DumperIOHelper::setParallelContext(bool is_parallel) { UInt whoami = Communicator::getStaticCommunicator().whoAmI(); UInt nb_proc = Communicator::getStaticCommunicator().getNbProc(); if (is_parallel) dumper->setParallelContext(whoami, nb_proc); else dumper->setParallelContext(0, 1); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::setDirectory(const std::string & directory) { this->directory = directory; dumper->setPrefix(directory); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::setBaseName(const std::string & basename) { filename = basename; } /* -------------------------------------------------------------------------- */ void DumperIOHelper::setTimeStep(Real time_step) { if (!time_activated) this->dumper->activateTimeDescFiles(time_step); else this->dumper->setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::dump() { try { dumper->dump(filename, count); } catch (iohelper::IOHelperException & e) { AKANTU_ERROR( "I was not able to dump your data with a Dumper: " << e.what()); } ++count; } /* -------------------------------------------------------------------------- */ void DumperIOHelper::dump(UInt step) { this->count = step; this->dump(); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::dump(Real current_time, UInt step) { this->dumper->setCurrentTime(current_time); this->dump(step); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::registerMesh(const Mesh & mesh, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { #if defined(AKANTU_IGFEM) if (element_kind == _ek_igfem) { registerField("connectivities", new dumper::IGFEMConnectivityField( mesh.getConnectivities(), spatial_dimension, ghost_type)); } else #endif registerField("connectivities", new dumper::ElementalField(mesh.getConnectivities(), spatial_dimension, ghost_type, element_kind)); registerField("positions", new dumper::NodalField(mesh.getNodes())); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::registerFilteredMesh( const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { auto * f_connectivities = new ElementTypeMapArrayFilter( mesh.getConnectivities(), elements_filter); this->registerField("connectivities", new dumper::FilteredConnectivityField( *f_connectivities, nodes_filter, spatial_dimension, ghost_type, element_kind)); this->registerField("positions", new dumper::NodalField( mesh.getNodes(), 0, 0, &nodes_filter)); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::registerField(const std::string & field_id, dumper::Field * field) { auto it = fields.find(field_id); if (it != fields.end()) { AKANTU_DEBUG_WARNING( "The field " << field_id << " is already registered in this Dumper. Field ignored."); return; } fields[field_id] = field; field->registerToDumper(field_id, *dumper); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::unRegisterField(const std::string & field_id) { auto it = fields.find(field_id); if (it == fields.end()) { AKANTU_DEBUG_WARNING( "The field " << field_id << " is not registered in this Dumper. Nothing to do."); return; } delete it->second; fields.erase(it); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::registerVariable(const std::string & variable_id, dumper::VariableBase * variable) { auto it = variables.find(variable_id); if (it != variables.end()) { AKANTU_DEBUG_WARNING( "The Variable " << variable_id << " is already registered in this Dumper. Variable ignored."); return; } variables[variable_id] = variable; variable->registerToDumper(variable_id, *dumper); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::unRegisterVariable(const std::string & variable_id) { auto it = variables.find(variable_id); if (it == variables.end()) { AKANTU_DEBUG_WARNING( "The variable " << variable_id << " is not registered in this Dumper. Nothing to do."); return; } delete it->second; variables.erase(it); } /* -------------------------------------------------------------------------- */ template iohelper::ElemType getIOHelperType() { AKANTU_TO_IMPLEMENT(); return iohelper::MAX_ELEM_TYPE; } template <> iohelper::ElemType getIOHelperType<_point_1>() { return iohelper::POINT_SET; } template <> iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } template <> iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } template <> iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } template <> iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } template <> iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } template <> iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } template <> iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } template <> iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } template <> iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } template <> iohelper::ElemType getIOHelperType<_hexahedron_20>() { return iohelper::HEX2; } template <> iohelper::ElemType getIOHelperType<_pentahedron_6>() { return iohelper::PRISM1; } template <> iohelper::ElemType getIOHelperType<_pentahedron_15>() { return iohelper::PRISM2; } #if defined(AKANTU_COHESIVE_ELEMENT) template <> iohelper::ElemType getIOHelperType<_cohesive_1d_2>() { return iohelper::COH1D2; } template <> iohelper::ElemType getIOHelperType<_cohesive_2d_4>() { return iohelper::COH2D4; } template <> iohelper::ElemType getIOHelperType<_cohesive_2d_6>() { return iohelper::COH2D6; } template <> iohelper::ElemType getIOHelperType<_cohesive_3d_6>() { return iohelper::COH3D6; } template <> iohelper::ElemType getIOHelperType<_cohesive_3d_12>() { return iohelper::COH3D12; } template <> iohelper::ElemType getIOHelperType<_cohesive_3d_8>() { return iohelper::COH3D8; } // template <> // iohelper::ElemType getIOHelperType<_cohesive_3d_16>() { return // iohelper::COH3D16; } #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) template <> iohelper::ElemType getIOHelperType<_bernoulli_beam_2>() { return iohelper::BEAM2; } template <> iohelper::ElemType getIOHelperType<_bernoulli_beam_3>() { return iohelper::BEAM3; } #endif /* -------------------------------------------------------------------------- */ UInt getIOHelperType(ElementType type) { UInt ioh_type = iohelper::MAX_ELEM_TYPE; #define GET_IOHELPER_TYPE(type) ioh_type = getIOHelperType(); AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_IOHELPER_TYPE); #undef GET_IOHELPER_TYPE return ioh_type; } /* -------------------------------------------------------------------------- */ } // akantu namespace iohelper { -template <> DataType getDataType() { return _int; } +template <> DataType getDataType() { + return getDataType>(); +} } diff --git a/src/io/dumper/dumper_text.cc b/src/io/dumper/dumper_text.cc index acd6603e4..e379a6683 100644 --- a/src/io/dumper/dumper_text.cc +++ b/src/io/dumper/dumper_text.cc @@ -1,112 +1,112 @@ /** * @file dumper_text.cc * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Nov 07 2017 * * @brief implementation of text dumper * * @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 "dumper_text.hh" #include "communicator.hh" #include "dumper_nodal_field.hh" #include "mesh.hh" #include namespace akantu { /* -------------------------------------------------------------------------- */ DumperText::DumperText(const std::string & basename, iohelper::TextDumpMode mode, bool parallel) : DumperIOHelper() { AKANTU_DEBUG_IN(); iohelper::DumperText * dumper_text = new iohelper::DumperText(mode); this->dumper = dumper_text; this->setBaseName(basename); this->setParallelContext(parallel); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DumperText::registerMesh(const Mesh & mesh, __attribute__((unused)) UInt spatial_dimension, __attribute__((unused)) const GhostType & ghost_type, __attribute__((unused)) const ElementKind & element_kind) { registerField("position", new dumper::NodalField(mesh.getNodes())); // in parallel we need node type UInt nb_proc = mesh.getCommunicator().getNbProc(); if (nb_proc > 1) { registerField("nodes_type", - new dumper::NodalField(mesh.getNodesType())); + new dumper::NodalField(mesh.getNodesFlags())); } } /* -------------------------------------------------------------------------- */ void DumperText::registerFilteredMesh( const Mesh & mesh, __attribute__((unused)) const ElementTypeMapArray & elements_filter, const Array & nodes_filter, __attribute__((unused)) UInt spatial_dimension, __attribute__((unused)) const GhostType & ghost_type, __attribute__((unused)) const ElementKind & element_kind) { registerField("position", new dumper::NodalField( mesh.getNodes(), 0, 0, &nodes_filter)); // in parallel we need node type UInt nb_proc = mesh.getCommunicator().getNbProc(); if (nb_proc > 1) { - registerField("nodes_type", new dumper::NodalField( - mesh.getNodesType(), 0, 0, &nodes_filter)); + registerField("nodes_type", new dumper::NodalField( + mesh.getNodesFlags(), 0, 0, &nodes_filter)); } } /* -------------------------------------------------------------------------- */ void DumperText::setBaseName(const std::string & basename) { AKANTU_DEBUG_IN(); DumperIOHelper::setBaseName(basename); static_cast(this->dumper) ->setDataSubDirectory(this->filename + "-DataFiles"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DumperText::setPrecision(UInt prec) { AKANTU_DEBUG_IN(); static_cast(this->dumper)->setPrecision(prec); AKANTU_DEBUG_OUT(); } } // akantu diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index 47e2090a8..fbec19b3d 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,567 +1,574 @@ /** * @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), connectivities("connectivities", id, memory_id), ghosts_counters("ghosts_counters", id, memory_id), normals("normals", id, memory_id), spatial_dimension(spatial_dimension), - size(spatial_dimension, 0.), - bbox(spatial_dimension), bbox_local(spatial_dimension), - mesh_data("mesh_data", id, memory_id), communicator(&communicator) { + size(spatial_dimension, 0.), bbox(spatial_dimension), + bbox_local(spatial_dimension), mesh_data("mesh_data", id, memory_id), + communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); + this->nodes_flags = std::make_shared>(0, 1, NodeFlag::_normal, + id + ":nodes_flags"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id, memory_id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const { barycenter.resize(getNbElement(type, ghost_type)); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type}, std::get<1>(data)); } } /* -------------------------------------------------------------------------- */ 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_flags = this->nodes_flags; mesh_facets->nodes_global_ids = this->nodes_global_ids; MeshUtils::buildAllFacets(*this, *mesh_facets, 0); if (mesh.isDistributed()) { mesh_facets->is_distributed = true; mesh_facets->element_synchronizer = std::make_unique( *mesh_facets, mesh.getElementSynchronizer()); } /// 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"); } auto & mesh_phys_data = this->getData("physical_names"); auto & phys_data = mesh_facets->getData("physical_names"); phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); ElementTypeMapArray barycenters(getID(), "temporary_barycenters"); barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); for (auto && ghost_type : ghost_types) { for (auto && type : barycenters.elementTypes(spatial_dimension - 1, ghost_type)) { mesh_facets->getBarycenters(barycenters(type, ghost_type), type, ghost_type); } } for_each_element( mesh, [&](auto && element) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); auto norm_barycenter = barycenter.norm(); auto tolerance = Math::getTolerance(); if (norm_barycenter > tolerance) tolerance *= norm_barycenter; const auto & element_to_facet = mesh_facets->getElementToSubelement( element.type, element.ghost_type); Vector barycenter_facet(spatial_dimension); auto range = enumerate(make_view(barycenters(element.type, element.ghost_type), spatial_dimension)); #ifndef AKANTU_NDEBUG auto min_dist = std::numeric_limits::max(); #endif // this is a spacial search coded the most inefficient way. auto facet = std::find_if(range.begin(), range.end(), [&](auto && data) { auto facet = std::get<0>(data); if (element_to_facet(facet)[1] == ElementNull) return false; auto norm_distance = barycenter.distance(std::get<1>(data)); #ifndef AKANTU_NDEBUG min_dist = std::min(min_dist, norm_distance); #endif return (norm_distance < tolerance); }); if (facet == range.end()) { AKANTU_DEBUG_INFO("The element " << element << " did not find its associated facet in the " "mesh_facets! Try to decrease math tolerance. " "The closest element was at a distance of " << min_dist); return; } // set physical name phys_data(Element{element.type, UInt(std::get<0>(*facet)), element.ghost_type}) = mesh_phys_data(element); }, _spatial_dimension = spatial_dimension - 1); mesh_facets->createGroupsFromMeshData("physical_names"); } AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { + + AKANTU_DEBUG_ASSERT(not is_distributed, + "You cannot read a mesh that is already distributed"); + MeshIO mesh_io; mesh_io.read(filename, *this, mesh_io_type); - type_iterator it = - this->firstType(spatial_dimension, _not_ghost, _ek_not_defined); - type_iterator last = - this->lastType(spatial_dimension, _not_ghost, _ek_not_defined); + auto it = this->firstType(spatial_dimension, _not_ghost, _ek_not_defined); + auto last = this->lastType(spatial_dimension, _not_ghost, _ek_not_defined); if (it == last) AKANTU_EXCEPTION( "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_flags->resize(nodes->size(), NodeFlag::_normal); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::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(); bbox_local.reset(); - for(auto & pos : make_view(*nodes, spatial_dimension)) { + for (auto & pos : make_view(*nodes, spatial_dimension)) { // if(!isPureGhostNode(i)) bbox_local += pos; } if (this->is_distributed) { bbox = bbox_local.allSum(*communicator); } else { bbox = bbox_local; } size = bbox.size(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : global_connectivity.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_not_defined, _ghost_type = ghost_type)) { if (not connectivities.exists(type, ghost_type)) continue; auto & local_conn = connectivities(type, ghost_type); auto & g_connectivity = global_connectivity(type, ghost_type); UInt nb_nodes = local_conn.size() * local_conn.getNbComponent(); std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), - [&](UInt l) -> UInt { - return this->getNodeGlobalId(l); - }); + [&](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"); + if (psize > 1) { + AKANTU_DEBUG_ERROR( + "Cannot distribute a mesh without a partitioning tool"); } #endif this->computeBoundingBox(); } - this->is_distributed = true; + if (psize > 1) + this->is_distributed = true; + + this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const Array & node_list, Array & elements) { for (const auto & node : node_list) for (const auto & element : *nodes_to_elements[node]) elements.push_back(element); } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements() { Element e; UInt nb_nodes = nodes->size(); for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) this->nodes_to_elements[n]->clear(); else this->nodes_to_elements[n] = std::make_unique>(); } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); Array::const_iterator> conn_it = connectivities(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) nodes_to_elements[conn(n)]->insert(e); } } } } /* -------------------------------------------------------------------------- */ std::tuple Mesh::updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { if (global_data_updater) return this->global_data_updater->updateData(nodes_event, elements_event); else { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } } /* -------------------------------------------------------------------------- */ void Mesh::registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->global_data_updater = std::move(global_data_updater); } /* -------------------------------------------------------------------------- */ void Mesh::eraseElements(const Array & elements) { ElementTypeMap last_element; RemovedElementsEvent event(*this); 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 2b10d7ae0..7bbed509d 100644 --- a/src/mesh/mesh.hh +++ b/src/mesh/mesh.hh @@ -1,651 +1,654 @@ /** * @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_bbox.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 #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 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 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(); /// set the periodicity in a given direction void makePeriodic(const SpatialDirection & direction); protected: void makePeriodic(const SpatialDirection & direction, const Array & list_1, const Array & list_2); public: /// 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 &); + AKANTU_GET_MACRO(NodesFlags, *nodes_flags, const Array &); protected: - AKANTU_GET_MACRO_NOT_CONST(NodesType, *nodes_type, Array &); + AKANTU_GET_MACRO_NOT_CONST(NodesFlags, *nodes_flags, Array &); public: - inline NodeType getNodeType(UInt local_id) const; + inline NodeFlag getNodeFlag(UInt local_id) const; + inline Int getNodePrank(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; + inline bool isPeriodicSlave(UInt n) const; + inline bool isPeriodicMaster(UInt n) const; + const Vector & getLowerBounds() const { return bbox.getLowerBounds(); } const Vector & getUpperBounds() const { return bbox.getUpperBounds(); } AKANTU_GET_MACRO(BBox, bbox, const BBox &); const Vector & getLocalLowerBounds() const { return bbox_local.getLowerBounds(); } const Vector & getLocalUpperBounds() const { return bbox_local.getUpperBounds(); } AKANTU_GET_MACRO(LocalBBox, bbox_local, const BBox &); /// 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; } /// defines if the mesh is periodic or not inline bool isPeriodic() const { return (this->is_periodic != 0); } inline bool isPeriodic(const SpatialDirection & direction) const { return ((this->is_periodic & (1 << direction)) != 0); } #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 const std::unordered_multimap & getPeriodicMasterSlaves() const { return periodic_master_slave; } /* ------------------------------------------------------------------------ */ /* 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(); + inline Array & getNodesFlagsPointer(); /// 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; + /// node flags (shared/periodic/...) + std::shared_ptr> nodes_flags; /// processor handling the node when not local or master std::unordered_map nodes_prank; /// 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}; /// size covered by the mesh on each direction Vector size; /// global bounding box BBox bbox; /// local bounding box BBox bbox_local; /// Extra data loaded from the mesh file 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}; /// defines if the mesh is periodic (3bits, 1 per direction) char is_periodic; /// 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; /// periodicity local info std::unordered_map periodic_slave_master; std::unordered_multimap periodic_master_slave; }; /// 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..2e00ed348 100644 --- a/src/mesh/mesh_accessor.hh +++ b/src/mesh/mesh_accessor.hh @@ -1,150 +1,153 @@ /** * @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(); } + inline auto & getNodesFlags() { return this->_mesh.getNodesFlags(); } + + /// get a pointer to the nodes_type Array and create it if necessary + inline void setNodePrank(UInt node, Int prank) { this->_mesh.nodes_prank[node] = prank; } /// 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_inline_impl.cc b/src/mesh/mesh_inline_impl.cc index 9c7a5531c..129907cae 100644 --- a/src/mesh/mesh_inline_impl.cc +++ b/src/mesh/mesh_inline_impl.cc @@ -1,688 +1,691 @@ /** * @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(); - + this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal); 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) + if (nodes_global_ids and not is_mesh_facets) this->removeNodesFromArray(*nodes_global_ids, new_numbering); - if (nodes_type and not mesh_parent) - this->removeNodesFromArray(*nodes_type, new_numbering); + if (not is_mesh_facets) + this->removeNodesFromArray(*nodes_flags, 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( 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( 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); } /* -------------------------------------------------------------------------- */ 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); } /* -------------------------------------------------------------------------- */ template inline const ElementTypeMapArray & Mesh::getData(const ID & data_name) const { return mesh_data.getElementalData(data_name); } /* -------------------------------------------------------------------------- */ template inline ElementTypeMapArray & Mesh::getData(const ID & data_name) { return mesh_data.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; + return ((*nodes_flags)(n) & NodeFlag::_shared_mask) == NodeFlag::_pure_ghost; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { - return nodes_type - ? ((*nodes_type)(n) == _nt_master) || - ((*nodes_type)(n) == _nt_normal) - : true; + return ((*nodes_flags)(n) & NodeFlag::_local_master_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { - return nodes_type ? (*nodes_type)(n) == _nt_normal : true; + return ((*nodes_flags)(n) & NodeFlag::_shared_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { - return nodes_type ? (*nodes_type)(n) == _nt_master : false; + return ((*nodes_flags)(n) & NodeFlag::_shared_mask) == NodeFlag::_master; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { - return nodes_type ? (*nodes_type)(n) >= 0 : false; + return ((*nodes_flags)(n) & NodeFlag::_shared_mask) == NodeFlag::_slave; +} + +/* -------------------------------------------------------------------------- */ +inline bool Mesh::isPeriodicSlave(UInt n) const { + return ((*nodes_flags)(n) & NodeFlag::_periodic_mask) == + NodeFlag::_periodic_slave; +} + +/* -------------------------------------------------------------------------- */ +inline bool Mesh::isPeriodicMaster(UInt n) const { + return ((*nodes_flags)(n) & NodeFlag::_periodic_mask) == + NodeFlag::_periodic_master; +} + +/* -------------------------------------------------------------------------- */ +inline NodeFlag Mesh::getNodeFlag(UInt local_id) const { + return (*nodes_flags)(local_id); } /* -------------------------------------------------------------------------- */ -inline NodeType Mesh::getNodeType(UInt local_id) const { - return nodes_type ? (*nodes_type)(local_id) : _nt_normal; +inline Int Mesh::getNodePrank(UInt local_id) const { + auto it = nodes_prank.find(local_id); + return it == nodes_prank.end() ? -1 : it->second; } /* -------------------------------------------------------------------------- */ 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; - for(const auto & el : elements) { + for (const auto & el : elements) { 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/mesh_periodic.cc b/src/mesh/mesh_periodic.cc index da8f4b4ee..5afccf2c8 100644 --- a/src/mesh/mesh_periodic.cc +++ b/src/mesh/mesh_periodic.cc @@ -1,335 +1,436 @@ /** * @file mesh_pbc.cc * * @author Nicolas Richart * * @date creation Sat Feb 10 2018 * * @brief Implementation of the perdiodicity capabilities in the mesh * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "communication_tag.hh" #include "communicator.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void Mesh::makePeriodic(const SpatialDirection & direction) { Array list_1; Array list_2; Real tolerance = 1e-10; auto lower_bound = this->getLowerBounds(); auto upper_bound = this->getUpperBounds(); auto length = upper_bound(direction) - lower_bound(direction); const auto & positions = *nodes; std::cout << bbox << std::endl; for (auto && data : enumerate(make_view(positions, spatial_dimension))) { UInt node = std::get<0>(data); const auto & pos = std::get<1>(data); if (std::abs((pos(direction) - lower_bound(direction)) / length) < tolerance) { list_1.push_back(node); } if (std::abs((pos(direction) - upper_bound(direction)) / length) < tolerance) { list_2.push_back(node); } } this->makePeriodic(direction, list_1, list_2); } namespace { struct NodeInfo { NodeInfo() {} NodeInfo(UInt spatial_dimension) : position(spatial_dimension) {} NodeInfo(UInt node, const Vector & position, const SpatialDirection & direction) : node(node), position(position) { this->direction_position = position(direction); this->position(direction) = 0.; } NodeInfo(const NodeInfo & other) : node(other.node), position(other.position), direction_position(other.direction_position) {} UInt node{0}; Vector position; Real direction_position{0.}; }; // std::ostream & operator<<(std::ostream & stream, const NodeInfo & info) { // stream << info.node << " " << info.position << " " // << info.direction_position; // return stream; // } } /* -------------------------------------------------------------------------- */ // left is for lower values on direction and right for highest values void Mesh::makePeriodic(const SpatialDirection & direction, const Array & list_left, const Array & list_right) { Real tolerance = 1e-10; const auto & positions = *nodes; auto lower_bound = this->getLowerBounds(); auto upper_bound = this->getUpperBounds(); auto length = upper_bound(direction) - lower_bound(direction); lower_bound(direction) = 0; upper_bound(direction) = 0; auto prank = communicator->whoAmI(); std::cout << prank << " - left:" << list_left.size() << " - right:" << list_right.size() << std::endl; std::vector nodes_left(list_left.size()); std::vector nodes_right(list_right.size()); BBox bbox(spatial_dimension); auto to_position = [&](UInt node) { Vector pos(spatial_dimension); for (UInt s : arange(spatial_dimension)) { pos(s) = positions(node, s); } auto && info = NodeInfo(node, pos, direction); bbox += info.position; return info; }; std::transform(list_left.begin(), list_left.end(), nodes_left.begin(), to_position); BBox bbox_left = bbox; bbox.reset(); std::transform(list_right.begin(), list_right.end(), nodes_right.begin(), to_position); BBox bbox_right = bbox; std::vector new_nodes; if (is_distributed) { + /* ---------------------------------------------------------------------- */ auto extract_and_send_nodes = [&](const auto & bbox, const auto & node_list, auto & send_buffers, auto proc, auto cnt) { send_buffers.emplace_back(); auto & buffer = send_buffers.back(); for (auto & info : node_list) { if (bbox.contains(info.position) and isLocalOrMasterNode(info.node)) { Vector pos = info.position; pos(direction) = info.direction_position; buffer << getNodeGlobalId(info.node); buffer << pos; + + buffer << ((*nodes_flags)(info.node) & NodeFlag::_periodic_mask); + if (isPeriodicSlave(info.node)) { + buffer << getNodeGlobalId(periodic_slave_master[info.node]); + } + if (isPeriodicMaster(info.node)) { + buffer << periodic_master_slave.count(info.node); + auto slaves = periodic_master_slave.equal_range(info.node); + for (auto it = slaves.first; it != slaves.second; ++it) { + buffer << getNodeGlobalId(it->second); + } + } } } auto tag = Tag::genTag(prank, cnt, Tag::_PERIODIC_NODES); return communicator->asyncSend(buffer, proc, tag); }; + /* ---------------------------------------------------------------------- */ auto recv_and_extract_nodes = [&](auto & bbox, auto & node_list, auto & buffer, const auto proc, auto cnt) { if (not bbox) return; buffer.reset(); auto tag = Tag::genTag(proc, cnt, Tag::_PERIODIC_NODES); communicator->receive(buffer, proc, tag); while (not buffer.empty()) { - NodeInfo info(spatial_dimension); Vector pos(spatial_dimension); UInt global_node; + NodeFlag flag; buffer >> global_node; - buffer >> pos; - info.position = pos; + buffer >> flag; + + auto local_node = getNodeLocalId(global_node); + + if (flag == NodeFlag::_periodic_slave) { + UInt master_node; + buffer >> master_node; + auto local_master_node = getNodeLocalId(master_node); + AKANTU_DEBUG_ASSERT(local_master_node != UInt(-1), + "Should I know the master node " << master_node); + } - info.direction_position = pos(direction); - info.position(direction) = 0; - // info.distance = lower_bound.distance(info.position); + if (flag == NodeFlag::_periodic_master) { + UInt nb_slaves; + buffer >> nb_slaves; + for (auto ns[[gnu::unused]] : arange(nb_slaves)) { + UInt gslave_node; + buffer >> gslave_node; + auto lslave_node = getNodeLocalId(gslave_node); + AKANTU_DEBUG_ASSERT(lslave_node != UInt(-1), + "Should I know the slave node " << gslave_node); + } + } - info.node = getNodeLocalId(global_node); - if (info.node != UInt(-1)) + if (local_node != UInt(-1)) continue; - info.node = nodes->size(); + local_node = nodes->size(); + + NodeInfo info(local_node, pos, direction); nodes->push_back(pos); nodes_global_ids->push_back(global_node); - nodes_type.push_back(NodeType(-100)); + nodes_flags->push_back(flag | NodeFlag::_pure_ghost); new_nodes.push_back(info.node); node_list.push_back(info); nodes_prank[info.node] = proc; } }; + /* ---------------------------------------------------------------------- */ + auto && intersection_right = bbox_left.intersection(bbox_right, *communicator); auto && intersection_left = bbox_right.intersection(bbox_left, *communicator); auto prank = communicator->whoAmI(); auto nb_proc = communicator->getNbProc(); using buffers_t = std::vector; std::vector send_requests; buffers_t send_buffers; for (auto && data : zip(arange(nb_proc), intersection_left, intersection_right)) { auto proc = std::get<0>(data); if (proc == prank) continue; buffers_t::iterator it; const auto & bbox_p_send_left = std::get<1>(data); if (bbox_p_send_left) { send_requests.push_back(extract_and_send_nodes( bbox_p_send_left, nodes_left, send_buffers, proc, 0)); } const auto & bbox_p_send_right = std::get<2>(data); if (bbox_p_send_right) { send_requests.push_back(extract_and_send_nodes( bbox_p_send_right, nodes_right, send_buffers, proc, 1)); } } DynamicCommunicationBuffer buffer; for (auto && data : zip(arange(nb_proc), intersection_left, intersection_right)) { auto proc = std::get<0>(data); if (proc == prank) continue; const auto & bbox_p_send_left = std::get<1>(data); recv_and_extract_nodes(bbox_p_send_left, nodes_left, buffer, proc, 0); const auto & bbox_p_send_right = std::get<2>(data); recv_and_extract_nodes(bbox_p_send_right, nodes_right, buffer, proc, 1); } communicator->waitAll(send_requests); communicator->freeCommunicationRequest(send_requests); } auto to_sort = [&](auto && info1, auto && info2) -> bool { return info1.position < info2.position; }; std::sort(nodes_left.begin(), nodes_left.end(), to_sort); std::sort(nodes_right.begin(), nodes_right.end(), to_sort); + auto register_pair = [&](auto & master, auto & slave) { + if (master == slave) + return; + + // if pair already registered + auto master_slaves = periodic_master_slave.equal_range(master); + auto slave_it = + std::find_if(master_slaves.first, master_slaves.second, + [&](auto & pair) { return pair.second == slave; }); + if (slave_it == master_slaves.second) { + // no duplicates + periodic_master_slave.insert(std::make_pair(master, slave)); + std::cout << "adding: " << getNodeGlobalId(master) << " - " + << getNodeGlobalId(slave) << std::endl; + } + + periodic_slave_master[slave] = master; + + (*nodes_flags)[slave] |= NodeFlag::_periodic_slave; + (*nodes_flags)[master] |= NodeFlag::_periodic_master; + }; + + auto updating_master = [&](auto & old_master, auto & new_master) { + if (old_master == new_master) + return; + + auto slaves = periodic_master_slave.equal_range(old_master); + AKANTU_DEBUG_ASSERT(slaves.first != periodic_master_slave.end(), + "Cannot update master " << old_master + << ", its not a master node!"); + std::cout << "updating: " << getNodeGlobalId(old_master) << " -> " + << getNodeGlobalId(new_master) << std::endl; + auto it = slaves.first; + for (; it != slaves.second; ++it) { + auto slave = it->second; + periodic_master_slave.insert(std::make_pair(new_master, slave)); + std::cout << " - " << getNodeGlobalId(new_master) << " - " + << getNodeGlobalId(slave) << std::endl; + /* might need the register_pair here */ + // \TODO tell processor prank[slave] that master is now node1 + // instead of node2 + periodic_slave_master[slave] = new_master; + // \TODO tell processor prank[new_master] that it also has a slave on + // prank[slave] + } + periodic_master_slave.erase(old_master); + (*nodes_flags)[old_master] |= NodeFlag::_periodic_slave; + }; + auto match_found = [&](auto & info1, auto & info2) { auto node1 = info1.node; auto node2 = info2.node; - if (info1.direction_position < info2.direction_position) { - std::swap(node1, node2); + + auto node2_master = node2; + if (isPeriodicSlave(node2)) { + node2_master = periodic_slave_master[node2]; } - auto mits = periodic_master_slave.equal_range(node2); - if (mits.first != periodic_master_slave.end()) { - for (auto mit = mits.first; mit != mits.second; ++mit) { - periodic_master_slave.insert(std::make_pair(node1, mit->second)); - // \TODO tell processor prank[mit->second] that master is now node1 - // instead of node2 - periodic_slave_master[mit->second] = node1; - // \TODO tell processor prank[node1] that it also has a slave on - // prank[mit->second] - } - periodic_master_slave.erase(node2); + auto master = node1; + bool node1_side_master = false; + if (isPeriodicMaster(node1)) { + node1_side_master = true; } - auto node1_slaves = periodic_master_slave.equal_range(node1); - auto slave_it = - std::find_if(node1_slaves.first, node1_slaves.second, - [&](auto & pair) { return pair.second == node2; }); - if (slave_it == node1_slaves.second) - periodic_master_slave.insert(std::make_pair(node1, node2)); + if (isPeriodicSlave(node1)) { + node1_side_master = true; + master = periodic_slave_master[node1]; + // register_pair(master, node1); + } + + if (node1_side_master) { + register_pair(master, node2); + + if (isPeriodicSlave(node2)) { + updating_master(node2_master, master); + return; + } - periodic_slave_master[node2] = node1; + if (isPeriodicMaster(node2)) { + updating_master(node2, master); + return; + } + } else { + if (isPeriodicSlave(node2)) { + register_pair(node2_master, node1); + return; + } + + if (isPeriodicMaster(node2)) { + register_pair(node2, node1); + return; + } + + register_pair(node1, node2); + } }; auto match_pairs = [&](auto & nodes_1, auto & nodes_2) { auto it = nodes_2.begin(); for (auto && info1 : nodes_1) { if (not isLocalOrMasterNode(info1.node)) continue; auto & pos1 = info1.position; auto it_cur = it; for (; it_cur != nodes_2.end(); ++it_cur) { auto & info2 = *it_cur; auto & pos2 = info2.position; auto dist = pos1.distance(pos2) / length; if (dist < tolerance) { match_found(info1, *it_cur); it = it_cur; break; } } } }; match_pairs(nodes_left, nodes_right); match_pairs(nodes_right, nodes_left); std::cout << periodic_slave_master.size() << std::endl; - for (auto & pair : periodic_master_slave) { - std::cout << prank << " - " << pair.first << " " << pair.second - << std::endl; - } + // for (auto & pair : periodic_master_slave) { + // std::cout << prank << " - " << pair.first << " " << pair.second + // << std::endl; + // } // std::cout << prank << " - Left" << std::endl; // for (auto && data : nodes_left) { // std::cout << prank << " - " << data << " -- " << getNodeType(data.node) // << std::endl; // } // std::cout << prank << " - Right" << std::endl; // for (auto && data : nodes_right) { // std::cout << prank << " - " << data << " -- " << getNodeType(data.node) // << std::endl; //} this->is_periodic |= 1 < direction; } } // akantu diff --git a/src/mesh_utils/global_ids_updater.cc b/src/mesh_utils/global_ids_updater.cc index 56eca0db9..ee0c72503 100644 --- a/src/mesh_utils/global_ids_updater.cc +++ b/src/mesh_utils/global_ids_updater.cc @@ -1,139 +1,139 @@ /** * @file global_ids_updater.cc * * @author Marco Vocialta * * @date creation: Fri Apr 13 2012 * @date last modification: Fri Dec 08 2017 * * @brief Functions of the GlobalIdsUpdater * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "global_ids_updater.hh" #include "element_synchronizer.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { UInt GlobalIdsUpdater::updateGlobalIDs(UInt local_nb_new_nodes) { if (mesh.getCommunicator().getNbProc() == 1) return local_nb_new_nodes; UInt total_nb_new_nodes = this->updateGlobalIDsLocally(local_nb_new_nodes); if (mesh.isDistributed()) { this->synchronizeGlobalIDs(); } return total_nb_new_nodes; } UInt GlobalIdsUpdater::updateGlobalIDsLocally(UInt local_nb_new_nodes) { const auto & comm = mesh.getCommunicator(); Int rank = comm.whoAmI(); Int nb_proc = comm.getNbProc(); if (nb_proc == 1) return local_nb_new_nodes; /// resize global ids array Array & nodes_global_ids = mesh.getGlobalNodesIds(); UInt old_nb_nodes = mesh.getNbNodes() - local_nb_new_nodes; nodes_global_ids.resize(mesh.getNbNodes(), -1); /// compute the number of global nodes based on the number of old nodes Matrix local_master_nodes(2, nb_proc, 0); for (UInt n = 0; n < old_nb_nodes; ++n) if (mesh.isLocalOrMasterNode(n)) ++local_master_nodes(0, rank); /// compute amount of local or master doubled nodes for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) if (mesh.isLocalOrMasterNode(n)) ++local_master_nodes(1, rank); comm.allGather(local_master_nodes); local_master_nodes = local_master_nodes.transpose(); UInt old_global_nodes = std::accumulate(local_master_nodes(0).storage(), local_master_nodes(0).storage() + nb_proc, 0); /// update global number of nodes UInt total_nb_new_nodes = std::accumulate(local_master_nodes(1).storage(), local_master_nodes(1).storage() + nb_proc, 0); if (total_nb_new_nodes == 0) return 0; /// set global ids of local and master nodes UInt starting_index = std::accumulate(local_master_nodes(1).storage(), local_master_nodes(1).storage() + rank, old_global_nodes); for (UInt n = old_nb_nodes; n < mesh.getNbNodes(); ++n) { if (mesh.isLocalOrMasterNode(n)) { nodes_global_ids(n) = starting_index; ++starting_index; } else { nodes_global_ids(n) = -1; } } MeshAccessor mesh_accessor(mesh); mesh_accessor.setNbGlobalNodes(old_global_nodes + total_nb_new_nodes); return total_nb_new_nodes; } void GlobalIdsUpdater::synchronizeGlobalIDs() { this->reduce = true; this->synchronizer.slaveReductionOnce(*this, _gst_giu_global_conn); #ifndef AKANTU_NDEBUG - for (auto node : nodes_types) { - auto node_type = mesh.getNodeType(node.first); - if (node_type != _nt_pure_ghost) + for (auto node : nodes_flags) { + auto node_flag = mesh.getNodeFlag(node.first); + if (node_flag != NodeFlag::_pure_ghost) continue; auto n = 0u; for (auto & pair : node.second) { - if (std::get<1>(pair) == _nt_pure_ghost) + if (std::get<1>(pair) == NodeFlag::_pure_ghost) ++n; } if (n == node.second.size()) { AKANTU_DEBUG_WARNING( "The node " << n << "is ghost on all the neighboring processors"); } } #endif this->reduce = false; this->synchronizer.synchronizeOnce(*this, _gst_giu_global_conn); } } // akantu diff --git a/src/mesh_utils/global_ids_updater.hh b/src/mesh_utils/global_ids_updater.hh index e22c4bd4d..229397424 100644 --- a/src/mesh_utils/global_ids_updater.hh +++ b/src/mesh_utils/global_ids_updater.hh @@ -1,107 +1,107 @@ /** * @file global_ids_updater.hh * * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Oct 02 2015 * @date last modification: Fri Dec 08 2017 * * @brief Class that updates the global ids of new nodes that are * inserted in the mesh. The functions in this class must be called * after updating the node types * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GLOBAL_IDS_UPDATER_HH__ #define __AKANTU_GLOBAL_IDS_UPDATER_HH__ /* -------------------------------------------------------------------------- */ #include "data_accessor.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class ElementSynchronizer; } // akantu namespace akantu { class GlobalIdsUpdater : public DataAccessor { public: GlobalIdsUpdater(Mesh & mesh, ElementSynchronizer & synchronizer) : mesh(mesh), synchronizer(synchronizer) {} /// function to update and synchronize the global connectivity of /// new inserted nodes. It must be called after updating the node /// types. (It calls in sequence the functions /// updateGlobalIDsLocally and synchronizeGlobalIDs) UInt updateGlobalIDs(UInt local_nb_new_nodes); /// function to update the global connectivity (only locally) of new /// inserted nodes. It must be called after updating the node types. UInt updateGlobalIDsLocally(UInt local_nb_new_nodes); /// function to synchronize the global connectivity of new inserted /// nodes among the processors. It must be called after updating the /// node types. void synchronizeGlobalIDs(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ template inline void packUnpackGlobalConnectivity(CommunicationBuffer & buffer, const Array & elements) const; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// Reference to the mesh to update Mesh & mesh; /// distributed synchronizer to communicate the connectivity ElementSynchronizer & synchronizer; /// Tells if a reduction is taking place or not bool reduce{false}; - std::unordered_map>> nodes_types; + std::unordered_map>> nodes_flags; }; } // akantu #include "global_ids_updater_inline_impl.cc" #endif /* __AKANTU_GLOBAL_IDS_UPDATER_HH__ */ diff --git a/src/mesh_utils/global_ids_updater_inline_impl.cc b/src/mesh_utils/global_ids_updater_inline_impl.cc index 14f427431..4727529ca 100644 --- a/src/mesh_utils/global_ids_updater_inline_impl.cc +++ b/src/mesh_utils/global_ids_updater_inline_impl.cc @@ -1,128 +1,128 @@ /** * @file global_ids_updater_inline_impl.cc * * @author Marco Vocialta * * @date creation: Fri Oct 02 2015 * @date last modification: Sun Aug 13 2017 * * @brief Implementation of the inline functions of GlobalIdsUpdater * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "global_ids_updater.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ #define __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt GlobalIdsUpdater::getNbData(const Array & elements, const SynchronizationTag & tag) const { UInt size = 0; if (tag == _gst_giu_global_conn){ size += Mesh::getNbNodesPerElementList(elements) * sizeof(UInt); #ifndef AKANTU_NDEBUG - size += sizeof(NodeType) * Mesh::getNbNodesPerElementList(elements) + sizeof(int); + size += sizeof(NodeFlag) * Mesh::getNbNodesPerElementList(elements) + sizeof(int); #endif } return size; } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag != _gst_giu_global_conn) return; auto & global_nodes_ids = mesh.getGlobalNodesIds(); #ifndef AKANTU_NDEBUG buffer << int(mesh.getCommunicator().whoAmI()); #endif for (auto & element : elements) { /// get element connectivity const Vector current_conn = const_cast(mesh).getConnectivity(element); /// loop on all connectivity nodes for (auto node : current_conn) { UInt index = -1; if ((this->reduce and mesh.isLocalOrMasterNode(node)) or (not this->reduce and not mesh.isPureGhostNode(node))) { index = global_nodes_ids(node); } buffer << index; #ifndef AKANTU_NDEBUG - buffer << mesh.getNodeType(node); + buffer << mesh.getNodeFlag(node); #endif } } } /* -------------------------------------------------------------------------- */ inline void GlobalIdsUpdater::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag != _gst_giu_global_conn) return; auto & global_nodes_ids = mesh.getGlobalNodesIds(); #ifndef AKANTU_NDEBUG int proc; buffer >> proc; #endif for (auto & element : elements) { /// get element connectivity Vector current_conn = const_cast(mesh).getConnectivity(element); /// loop on all connectivity nodes for (auto node : current_conn) { UInt index; buffer >> index; #ifndef AKANTU_NDEBUG - NodeType node_type; - buffer >> node_type; + NodeFlag node_flag; + buffer >> node_flag; if (reduce) - nodes_types[node].push_back(std::make_pair(proc, node_type)); + nodes_flags[node].push_back(std::make_pair(proc, node_flag)); #endif if (index == UInt(-1)) continue; if (mesh.isSlaveNode(node)) global_nodes_ids(node) = index; } } } } // akantu #endif /* __AKANTU_GLOBAL_IDS_UPDATER_INLINE_IMPL_CC__ */ diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index e617058dc..64f3090b3 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,886 +1,884 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Thu Feb 08 2018 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager_default.hh" #include "communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver_default.hh" #include "sparse_matrix_aij.hh" #include "terms_to_assemble.hh" #include "time_step_solver_default.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = i; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { if (c_jcn >= c_irn) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addElementalMatrixToUnsymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), - dofs_type(0, 1, std::string(id + ":dofs_type")), + dofs_flag(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), - global_equation_number(0, 1, "global_equation_number"), synchronizer(nullptr) {} /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), - dofs_type(0, 1, std::string(id + ":dofs_type")), + dofs_flag(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), - jacobian_release(0), global_equation_number(0, 1, "global_equation_number"), first_global_dof_id(0), synchronizer(nullptr) { if (this->mesh->isDistributed()) this->synchronizer = std::make_unique( *this, this->id + ":dof_synchronizer", this->memory_id); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() = default; /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = array_to_assemble.size() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.size() == nb_degree_of_freedoms, "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); typename Array::const_scalar_iterator arr_it = array_to_assemble.begin_reinterpret(nb_degree_of_freedoms); Array::const_scalar_iterator equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++arr_it, ++equ_it) { global_array(*equ_it) += scale_factor * (*arr_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id) : DOFData(dof_id), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManagerDefault::getNewDOFData(const ID & dof_id) { this->dofs[dof_id] = std::make_unique(dof_id); return *this->dofs[dof_id]; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs) { // auto prank = this->communicator.whoAmI(); // auto psize = this->communicator.getNbProc(); // access the relevant data to update auto & dof_data = this->getDOFDataTyped(dof_id); const auto & support_type = dof_data.support_type; const auto & group = dof_data.group_support; if (support_type == _dst_nodal and group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_dofs, nb_pure_local_dofs, [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs, [](UInt node) -> UInt { return node; }); } // update the synchronizer if needed if (this->synchronizer) this->synchronizer->registerDOFs(dof_id); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, support_type); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, group_support); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; std::unique_ptr sm = std::make_unique(*this, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrixAIJ & sm_to_copy = this->getMatrix(matrix_to_copy_id); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) { SparseMatrix & matrix = DOFManager::getMatrix(id); return dynamic_cast(matrix); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerDefault::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; std::unique_ptr nls; switch (type) { #if defined(AKANTU_IMPLICIT) case _nls_newton_raphson: case _nls_newton_raphson_modified: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } case _nls_linear: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } #endif case _nls_lumped: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } default: AKANTU_EXCEPTION("The asked type of non linear solver is not supported by " "this dof manager"); } return this->registerNonLinearSolver(non_linear_solver_id, nls); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique( *this, type, non_linear_solver, time_step_solver_id, this->memory_id); return this->registerTimeStepSolver(time_step_solver_id, tss); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearResidual() { this->residual.resize(this->local_system_size); this->residual.clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearMatrix(const ID & mtx) { this->getMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateGlobalBlockedDofs() { auto it = this->dofs.begin(); auto end = this->dofs.end(); this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.resize(this->local_system_size); this->global_blocked_dofs.clear(); for (; it != end; ++it) { if (!this->hasBlockedDOFs(it->first)) continue; DOFData & dof_data = *it->second; this->assembleToGlobalArray(it->first, *dof_data.blocked_dofs, this->global_blocked_dofs, true); } } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = equation_number.size(); local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent()); auto loc_it = local_array.begin_reinterpret(nb_degree_of_freedoms); auto equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++loc_it, ++equ_it) { (*loc_it) = global_array(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->global_solution, solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToResidual( const ID & dof_id, const Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); this->assembleToGlobalArray(dof_id, array_to_assemble, this->residual, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToLumpedMatrix( const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); Array & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { SparseMatrixAIJ & A = this->getMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, 1.); Array tmp_residual(this->residual.size(), 1, 0.); A.matVecMul(data_cache, tmp_residual, scale_factor, 1.); this->residual += tmp_residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleLumpedMatMulVectToResidual( const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { const Array & A = this->getLumpedMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, scale_factor); auto A_it = A.begin(); auto A_end = A.end(); auto x_it = data_cache.begin(); auto r_it = this->residual.begin(); for (; A_it != A_end; ++A_it, ++x_it, ++r_it) { *r_it += *A_it * *x_it; } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); DOFData & dof_data = this->getDOFData(dof_id); AKANTU_DEBUG_ASSERT(dof_data.support_type == _dst_nodal, "This function applies only on Nodal dofs"); this->addToProfile(matrix_id, dof_id, type, ghost_type); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt nb_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); nb_element = group_elements.size(); filter_it = group_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } } AKANTU_DEBUG_ASSERT(elementary_mat.size() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = dof_data.dof->getNbComponent(); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); auto conn_begin = connectivity.begin(nb_nodes_per_element); auto conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); std::transform(element_eq_nb.begin(), element_eq_nb.end(), element_eq_nb.begin(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (filter_it) ++filter_it; else ++conn_it; if (A.getMatrixType() == _symmetric) if (elemental_matrix_type == _symmetric) this->addSymmetricElementalMatrixToSymmetric(A, *el_mat_it, element_eq_nb, A.size()); else this->addUnsymmetricElementalMatrixToSymmetric(A, *el_mat_it, element_eq_nb, A.size()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb, A.size()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { const Array & equation_number_m = this->getLocalEquationNumbers(dof_id_m); const Array & equation_number_n = this->getLocalEquationNumbers(dof_id_n); SparseMatrixAIJ & A = this->getMatrix(matrix_id); for (const auto & term : terms) { UInt gi = this->localToGlobalEquationNumber(equation_number_m(term.i())); UInt gj = this->localToGlobalEquationNumber(equation_number_n(term.j())); A.add(gi, gj, term); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const DOFData & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) return; auto mat_dof = std::make_pair(matrix_id, dof_id); auto type_pair = std::make_pair(type, ghost_type); auto prof_it = this->matrix_profiled_dofs.find(mat_dof); if (prof_it != this->matrix_profiled_dofs.end() && std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) != prof_it->second.end()) return; UInt nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt size = A.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto & connectivity = this->mesh->getConnectivity(type, ghost_type); auto cbegin = connectivity.begin(nb_nodes_per_element); auto cit = cbegin; UInt nb_elements = connectivity.size(); UInt * ge_it = nullptr; if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); ge_it = group_elements.storage(); nb_elements = group_elements.size(); } UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; Vector element_eq_nb(size_mat); for (UInt e = 0; e < nb_elements; ++e) { if (ge_it) cit = cbegin + *ge_it; this->extractElementEquationNumber( equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (ge_it) ++ge_it; else ++cit; for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = element_eq_nb(i); if (c_irn < size) { for (UInt j = 0; j < size_mat; ++j) { UInt c_jcn = element_eq_nb(j); if (c_jcn < size) { A.add(c_irn, c_jcn); } } } } } this->matrix_profiled_dofs[mat_dof].push_back(type_pair); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::applyBoundary(const ID & matrix_id) { SparseMatrixAIJ & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getValueRelease()) { auto are_equal = std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (not are_equal) J.applyBoundary(); previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getValueRelease(); } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getGlobalResidual() { if (this->synchronizer) { if (not this->global_residual) { this->global_residual = std::make_unique>(0, 1, "global_residual"); } if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->global_residual->resize(this->system_size); this->synchronizer->gather(this->residual, *this->global_residual); } else { this->synchronizer->gather(this->residual); } return *this->global_residual; } else { return this->residual; } } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getResidual() const { return this->residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::setGlobalSolution(const Array & solution) { if (this->synchronizer) { if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->synchronizer->scatter(this->global_solution, solution); } else { this->synchronizer->scatter(this->global_solution); } } else { AKANTU_DEBUG_ASSERT(solution.size() == this->global_solution.size(), "Sequential call to this function needs the solution " "to be the same size as the global_solution"); this->global_solution.copy(solution); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { DOFManager::onNodesAdded(nodes_list, event); if (this->synchronizer) this->synchronizer->onNodesAdded(nodes_list); } /* -------------------------------------------------------------------------- */ std::pair DOFManagerDefault::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { UInt nb_new_local_dofs, nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = DOFManager::updateNodalDOFs(dof_id, nodes_list); auto & dof_data = this->getDOFDataTyped(dof_id); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor(DOFManagerDefault::DOFDataDefault & dof_data, DOFManagerDefault & dof_manager) : dof_data(dof_data), dof_manager(dof_manager) { for (auto && pair : zip(dof_data.local_equation_number, dof_data.associated_nodes)) { UInt node, dof; std::tie(dof, node) = pair; dofs_per_node[node].push_back(dof); } } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_ask_nodes) { return nodes.size() * dof_data.dof->getNbComponent() * sizeof(UInt); } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_ask_nodes) { for (auto & node : nodes) { auto & dofs = dofs_per_node.at(node); for (auto & dof : dofs) { buffer << dof_manager.global_equation_number(dof); } } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { if (tag == _gst_ask_nodes) { for (auto & node : nodes) { auto & dofs = dofs_per_node[node]; for (auto dof : dofs) { UInt global_dof; buffer >> global_dof; dof_manager.global_equation_number(dof) = global_dof; dof_manager.global_to_local_mapping[global_dof] = dof; } } } } protected: std::unordered_map> dofs_per_node; DOFManagerDefault::DOFDataDefault & dof_data; DOFManagerDefault & dof_manager; }; /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateDOFsData( DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, const std::function & getNode) { auto prank = this->communicator.whoAmI(); auto psize = this->communicator.getNbProc(); // access the relevant data to update const auto & support_type = dof_data.support_type; // resize all relevant arrays this->residual.resize(this->local_system_size, 0.); - this->dofs_type.resize(this->local_system_size); + this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_solution.resize(this->local_system_size, 0.); this->global_blocked_dofs.resize(this->local_system_size, true); this->previous_global_blocked_dofs.resize(this->local_system_size, true); this->global_equation_number.resize(this->local_system_size, -1); for (auto & lumped_matrix : lumped_matrices) lumped_matrix.second->resize(this->local_system_size); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); // determine the first local/global dof id to use Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = nb_new_pure_local; this->communicator.allGather(nb_dofs_per_proc); this->first_global_dof_id += std::accumulate( nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0); UInt first_dof_id = this->local_system_size - nb_new_local_dofs; if (support_type == _dst_nodal) { dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_new_local_dofs); } // update per dof info for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number.push_back(local_eq_num); bool is_local_dof = true; // determine the dof type switch (support_type) { case _dst_nodal: { UInt node = getNode(d / dof_data.dof->getNbComponent()); - this->dofs_type(local_eq_num) = this->mesh->getNodeType(node); + this->dofs_flag(local_eq_num) = this->mesh->getNodeFlag(node); dof_data.associated_nodes.push_back(node); is_local_dof = this->mesh->isLocalOrMasterNode(node); break; } case _dst_generic: { - this->dofs_type(local_eq_num) = _nt_normal; + this->dofs_flag(local_eq_num) = NodeFlag::_normal; break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } // update global id for local dofs if (is_local_dof) { this->global_equation_number(local_eq_num) = this->first_global_dof_id; this->global_to_local_mapping[this->first_global_dof_id] = local_eq_num; ++this->first_global_dof_id; } else { this->global_equation_number(local_eq_num) = UInt(-1); } } // synchronize the global numbering for slaves if (support_type == _dst_nodal && this->synchronizer) { GlobalDOFInfoDataAccessor data_accessor(dof_data, *this); auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, _gst_ask_nodes); } this->first_global_dof_id += std::accumulate( nb_dofs_per_proc.begin() + prank + 1, nb_dofs_per_proc.end(), 0); matrix_profiled_dofs.clear(); for(auto & matrix : matrices) { matrix.second->clearProfile(); } } /* -------------------------------------------------------------------------- */ // register in factory static bool default_dof_manager_is_registered[[gnu::unused]] = DefaultDOFManagerFactory::getInstance().registerAllocator( "default", [](const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(id, mem_id); }); static bool dof_manager_is_registered[[gnu::unused]] = DOFManagerFactory::getInstance().registerAllocator( "default", [](Mesh & mesh, const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(mesh, id, mem_id); }); } // namespace akantu diff --git a/src/model/dof_manager_default.hh b/src/model/dof_manager_default.hh index 4fbf01e52..f6ca50501 100644 --- a/src/model/dof_manager_default.hh +++ b/src/model/dof_manager_default.hh @@ -1,356 +1,356 @@ /** * @file dof_manager_default.hh * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Jan 31 2018 * * @brief Default implementation of the dof manager * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_DEFAULT_HH__ #define __AKANTU_DOF_MANAGER_DEFAULT_HH__ namespace akantu { class SparseMatrixAIJ; class NonLinearSolverDefault; class TimeStepSolverDefault; class DOFSynchronizer; } // namespace akantu namespace akantu { class DOFManagerDefault : public DOFManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFManagerDefault(const ID & id = "dof_manager_default", const MemoryID & memory_id = 0); DOFManagerDefault(Mesh & mesh, const ID & id = "dof_manager_default", const MemoryID & memory_id = 0); ~DOFManagerDefault() override; protected: struct DOFDataDefault : public DOFData { explicit DOFDataDefault(const ID & dof_id); /// associated node for _dst_nodal dofs only Array associated_nodes; }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ private: void registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs); public: /// register an array of degree of freedom void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) override; /// the dof as an implied type of _dst_nodal and is defined only on a subset /// of nodes void registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support) override; /// Assemble an array to the global residual array void assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor = 1.) override; /// Assemble an array to the global lumped matrix array void assembleToLumpedMatrix(const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor = 1.) override; /** * Assemble elementary values to the global matrix. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) override; /// multiply a vector by a matrix and assemble the result to the residual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) override; /// multiply a vector by a lumped matrix and assemble the result to the /// residual void assembleLumpedMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) override; /// assemble coupling terms between to dofs void assemblePreassembledMatrix(const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) override; protected: /// Assemble an array to the global residual array template void assembleToGlobalArray(const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor); public: /// clear the residual void clearResidual() override; /// sets the matrix to 0 void clearMatrix(const ID & mtx) override; /// sets the lumped matrix to 0 void clearLumpedMatrix(const ID & mtx) override; /// update the global dofs vector virtual void updateGlobalBlockedDofs(); /// apply boundary conditions to jacobian matrix virtual void applyBoundary(const ID & matrix_id = "J"); // void getEquationsNumbers(const ID & dof_id, // Array & equation_numbers) override; protected: /// Get local part of an array corresponding to a given dofdata template void getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const; /// Get the part of the solution corresponding to the dof_id void getSolutionPerDOFs(const ID & dof_id, Array & solution_array) override; /// fill a Vector with the equation numbers corresponding to the given /// connectivity inline void extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & local_equation_number); public: /// extract a lumped matrix part corresponding to a given dof void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) override; private: /// Add a symmetric matrices to a symmetric sparse matrix void addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a unsymmetric matrices to a symmetric sparse matrix (i.e. cohesive /// elements) void addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a matrices to a unsymmetric sparse matrix void addElementalMatrixToUnsymmetric(SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); void addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ protected: std::pair updateNodalDOFs(const ID & dof_id, const Array & nodes_list) override; private: void updateDOFsData(DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, const std::function & getNode); public: /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get an instance of a new SparseMatrix SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) override; /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) override; /// Get the reference of an existing matrix SparseMatrixAIJ & getMatrix(const ID & matrix_id); /* ------------------------------------------------------------------------ */ /* Non Linear Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver NonLinearSolver & getNewNonLinearSolver( const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type) override; /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver TimeStepSolver & getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) override; /* ------------------------------------------------------------------------ */ /// Get the solution array AKANTU_GET_MACRO_NOT_CONST(GlobalSolution, global_solution, Array &); /// Set the global solution array void setGlobalSolution(const Array & solution); /// Get the global residual array across processors (SMP call) const Array & getGlobalResidual(); /// Get the residual array const Array & getResidual() const; /// Get the blocked dofs array AKANTU_GET_MACRO(GlobalBlockedDOFs, global_blocked_dofs, const Array &); /// Get the blocked dofs array AKANTU_GET_MACRO(PreviousGlobalBlockedDOFs, previous_global_blocked_dofs, const Array &); /// Get the location type of a given dof inline bool isLocalOrMasterDOF(UInt local_dof_num); /// Answer to the question is a dof a slave dof ? inline bool isSlaveDOF(UInt local_dof_num); /// get the equation numbers (in local numbering) corresponding to a dof ID inline const Array & getLocalEquationNumbers(const ID & dof_id) const; /// tells if the dof manager knows about a global dof bool hasGlobalEquationNumber(UInt global) const; /// return the local index of the global equation number inline UInt globalToLocalEquationNumber(UInt global) const; /// converts local equation numbers to global equation numbers; inline UInt localToGlobalEquationNumber(UInt local) const; /// get the array of dof types (use only if you know what you do...) - inline Int getDOFType(UInt local_id) const; + inline NodeFlag getDOFFlag(UInt local_id) const; /// get the array of dof types (use only if you know what you do...) inline const Array & getDOFsAssociatedNodes(const ID & dof_id) const; /// access the internal dof_synchronizer AKANTU_GET_MACRO_NOT_CONST(Synchronizer, *synchronizer, DOFSynchronizer &); /// access the internal dof_synchronizer bool hasSynchronizer() const { return synchronizer != nullptr; } protected: DOFData & getNewDOFData(const ID & dof_id) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend class GlobalDOFInfoDataAccessor; // using AIJMatrixMap = std::map>; // using DefaultNonLinearSolversMap = // std::map>; // using DefaultTimeStepSolversMap = // std::map>; using DOFToMatrixProfile = std::map, std::vector>>; /// contains the the dofs that where added to the profile of a given matrix. DOFToMatrixProfile matrix_profiled_dofs; /// rhs to the system of equation corresponding to the residual linked to the /// different dofs Array residual; /// rhs used only on root proc in case of parallel computing, this is the full /// gathered rhs array std::unique_ptr> global_residual; /// solution of the system of equation corresponding to the different dofs Array global_solution; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array global_blocked_dofs; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array previous_global_blocked_dofs; /// define the dofs type, local, shared, ghost - Array dofs_type; + Array dofs_flag; /// Memory cache, this is an array to keep the temporary memory needed for /// some operations, it is meant to be resized or cleared when needed Array data_cache; /// Release at last apply boundary on jacobian UInt jacobian_release{0}; /// equation number in global numbering Array global_equation_number; using equation_numbers_map = std::unordered_map; /// dual information of global_equation_number equation_numbers_map global_to_local_mapping; /// accumulator to know what would be the next global id to use UInt first_global_dof_id{0}; /// synchronizer to maintain coherency in dof fields std::unique_ptr synchronizer; }; } // namespace akantu #include "dof_manager_default_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_DEFAULT_HH__ */ diff --git a/src/model/dof_manager_default_inline_impl.cc b/src/model/dof_manager_default_inline_impl.cc index 26b764a7c..da2cb3302 100644 --- a/src/model/dof_manager_default_inline_impl.cc +++ b/src/model/dof_manager_default_inline_impl.cc @@ -1,105 +1,105 @@ /** * @file dof_manager_default_inline_impl.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Jan 31 2018 * * @brief Implementation of the DOFManagerDefault inline functions * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager_default.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_DEFAULT_INLINE_IMPL_CC__ #define __AKANTU_DOF_MANAGER_DEFAULT_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline bool DOFManagerDefault::isLocalOrMasterDOF(UInt dof_num) { - Int dof_type = this->dofs_type(dof_num); - return (dof_type == Int(_nt_normal)) || (dof_type == Int(_nt_master)); + auto dof_flag = this->dofs_flag(dof_num); + return (dof_flag & NodeFlag::_shared_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool DOFManagerDefault::isSlaveDOF(UInt dof_num) { - Int dof_type = this->dofs_type(dof_num); - return (dof_type >= 0); + auto dof_flag = this->dofs_flag(dof_num); + return (dof_flag & NodeFlag::_shared_mask) == NodeFlag::_slave; } /* -------------------------------------------------------------------------- */ inline const Array & DOFManagerDefault::getLocalEquationNumbers(const ID & dof_id) const { return this->getDOFData(dof_id).local_equation_number; } inline const Array & DOFManagerDefault::getDOFsAssociatedNodes(const ID & dof_id) const { const auto & dof_data = this->getDOFDataTyped(dof_id); return dof_data.associated_nodes; } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & element_equation_number) { for (UInt i = 0, ld = 0; i < connectivity.size(); ++i) { UInt n = connectivity(i); for (UInt d = 0; d < nb_degree_of_freedom; ++d, ++ld) { element_equation_number(ld) = equation_numbers(n * nb_degree_of_freedom + d); } } } /* -------------------------------------------------------------------------- */ inline UInt DOFManagerDefault::localToGlobalEquationNumber(UInt local) const { return this->global_equation_number(local); } /* -------------------------------------------------------------------------- */ inline bool DOFManagerDefault::hasGlobalEquationNumber(UInt global) const { auto it = this->global_to_local_mapping.find(global); return (it != this->global_to_local_mapping.end()); } /* -------------------------------------------------------------------------- */ inline UInt DOFManagerDefault::globalToLocalEquationNumber(UInt global) const { auto it = this->global_to_local_mapping.find(global); AKANTU_DEBUG_ASSERT(it != this->global_to_local_mapping.end(), "This global equation number " << global << " does not exists in " << this->id); return it->second; } /* -------------------------------------------------------------------------- */ -inline Int DOFManagerDefault::getDOFType(UInt local_id) const { - return this->dofs_type(local_id); +inline NodeFlag DOFManagerDefault::getDOFFlag(UInt local_id) const { + return this->dofs_flag(local_id); } /* -------------------------------------------------------------------------- */ } // akantu #endif /* __AKANTU_DOF_MANAGER_DEFAULT_INLINE_IMPL_CC_ */ diff --git a/src/synchronizer/communicator.cc b/src/synchronizer/communicator.cc index fc925697c..6256c57c2 100644 --- a/src/synchronizer/communicator.cc +++ b/src/synchronizer/communicator.cc @@ -1,181 +1,181 @@ /** * @file communicator.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Feb 05 2018 * * @brief implementation of the common part of the static communicator * * @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 "communicator.hh" #if defined(AKANTU_USE_MPI) #include "mpi_communicator_data.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { #if defined(AKANTU_USE_MPI) int MPICommunicatorData::is_externaly_initialized = 0; #endif UInt InternalCommunicationRequest::counter = 0; /* -------------------------------------------------------------------------- */ InternalCommunicationRequest::InternalCommunicationRequest(UInt source, UInt dest) : source(source), destination(dest) { this->id = counter++; } /* -------------------------------------------------------------------------- */ InternalCommunicationRequest::~InternalCommunicationRequest() = default; /* -------------------------------------------------------------------------- */ void InternalCommunicationRequest::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "CommunicationRequest [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + source : " << source << std::endl; stream << space << " + destination : " << destination << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ Communicator::~Communicator() { auto * event = new FinalizeCommunicatorEvent(*this); this->sendEvent(*event); delete event; } /* -------------------------------------------------------------------------- */ Communicator & Communicator::getStaticCommunicator() { AKANTU_DEBUG_IN(); if (!static_communicator) { int nb_args = 0; char ** null = nullptr; static_communicator = std::make_unique(nb_args, null, private_member{}); } AKANTU_DEBUG_OUT(); return *static_communicator; } /* -------------------------------------------------------------------------- */ Communicator & Communicator::getStaticCommunicator(int & argc, char **& argv) { if (!static_communicator) static_communicator = std::make_unique(argc, argv, private_member{}); return getStaticCommunicator(); } } // namespace akantu #ifdef AKANTU_USE_MPI #include "communicator_mpi_inline_impl.cc" #else #include "communicator_dummy_inline_impl.cc" #endif namespace akantu { /* -------------------------------------------------------------------------- */ /* Template instantiation */ /* -------------------------------------------------------------------------- */ #define AKANTU_COMM_INSTANTIATE(T) \ template void Communicator::probe(Int sender, Int tag, \ CommunicationStatus & status) const; \ template bool Communicator::asyncProbe( \ Int sender, Int tag, CommunicationStatus & status) const; \ template void Communicator::sendImpl( \ const T * buffer, Int size, Int receiver, Int tag, \ const CommunicationMode & mode) const; \ template void Communicator::receiveImpl(T * buffer, Int size, Int sender, \ Int tag) const; \ template CommunicationRequest Communicator::asyncSendImpl( \ const T * buffer, Int size, Int receiver, Int tag, \ const CommunicationMode & mode) const; \ template CommunicationRequest Communicator::asyncReceiveImpl( \ T * buffer, Int size, Int sender, Int tag) const; \ template void Communicator::allGatherImpl(T * values, int nb_values) \ const; \ template void Communicator::allGatherVImpl(T * values, int * nb_values) \ const; \ template void Communicator::gatherImpl(T * values, int nb_values, \ int root) const; \ template void Communicator::gatherImpl( \ T * values, int nb_values, T * gathered, int nb_gathered) const; \ template void Communicator::gatherVImpl(T * values, int * nb_values, \ int root) const; \ template void Communicator::broadcastImpl(T * values, int nb_values, \ int root) const; \ template void Communicator::allReduceImpl( \ T * values, int nb_values, const SynchronizerOperation & op) const #define MIN_MAX_REAL SCMinMaxLoc AKANTU_COMM_INSTANTIATE(bool); AKANTU_COMM_INSTANTIATE(Real); AKANTU_COMM_INSTANTIATE(UInt); AKANTU_COMM_INSTANTIATE(Int); AKANTU_COMM_INSTANTIATE(char); -AKANTU_COMM_INSTANTIATE(NodeType); +AKANTU_COMM_INSTANTIATE(NodeFlag); AKANTU_COMM_INSTANTIATE(MIN_MAX_REAL); #if AKANTU_INTEGER_SIZE > 4 AKANTU_COMM_INSTANTIATE(int); #endif // template void Communicator::send>( // SCMinMaxLoc * buffer, Int size, Int receiver, Int tag); // template void Communicator::receive>( // SCMinMaxLoc * buffer, Int size, Int sender, Int tag); // template CommunicationRequest // Communicator::asyncSend>( // SCMinMaxLoc * buffer, Int size, Int receiver, Int tag); // template CommunicationRequest // Communicator::asyncReceive>( // SCMinMaxLoc * buffer, Int size, Int sender, Int tag); // template void Communicator::probe>( // Int sender, Int tag, CommunicationStatus & status); // template void Communicator::allGather>( // SCMinMaxLoc * values, int nb_values); // template void Communicator::allGatherV>( // SCMinMaxLoc * values, int * nb_values); // template void Communicator::gather>( // SCMinMaxLoc * values, int nb_values, int root); // template void Communicator::gatherV>( // SCMinMaxLoc * values, int * nb_values, int root); // template void Communicator::broadcast>( // SCMinMaxLoc * values, int nb_values, int root); // template void Communicator::allReduce>( // SCMinMaxLoc * values, int nb_values, // const SynchronizerOperation & op); } // akantu diff --git a/src/synchronizer/communicator_mpi_inline_impl.cc b/src/synchronizer/communicator_mpi_inline_impl.cc index fda048047..809108db4 100644 --- a/src/synchronizer/communicator_mpi_inline_impl.cc +++ b/src/synchronizer/communicator_mpi_inline_impl.cc @@ -1,459 +1,460 @@ /** * @file communicator_mpi_inline_impl.cc * * @author Nicolas Richart * * @date creation: Tue Nov 07 2017 * @date last modification: Mon Dec 18 2017 * * @brief StaticCommunicatorMPI implementation * * @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 "mpi_communicator_data.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #if (defined(__GNUC__) || defined(__GNUG__)) #if AKA_GCC_VERSION < 60000 namespace std { template <> struct hash { using argument_type = akantu::SynchronizerOperation; size_t operator()(const argument_type & e) const noexcept { auto ue = underlying_type_t(e); return uh(ue); } private: const hash> uh{}; }; } // namespace std #endif #endif namespace akantu { class CommunicationRequestMPI : public InternalCommunicationRequest { public: CommunicationRequestMPI(UInt source, UInt dest) : InternalCommunicationRequest(source, dest), request(std::make_unique()) {} MPI_Request & getMPIRequest() { return *request; }; private: std::unique_ptr request; }; namespace { template inline MPI_Datatype getMPIDatatype(); MPI_Op getMPISynchronizerOperation(const SynchronizerOperation & op) { std::unordered_map _operations{ {SynchronizerOperation::_sum, MPI_SUM}, {SynchronizerOperation::_min, MPI_MIN}, {SynchronizerOperation::_max, MPI_MAX}, {SynchronizerOperation::_prod, MPI_PROD}, {SynchronizerOperation::_land, MPI_LAND}, {SynchronizerOperation::_band, MPI_BAND}, {SynchronizerOperation::_lor, MPI_LOR}, {SynchronizerOperation::_bor, MPI_BOR}, {SynchronizerOperation::_lxor, MPI_LXOR}, {SynchronizerOperation::_bxor, MPI_BXOR}, {SynchronizerOperation::_min_loc, MPI_MINLOC}, {SynchronizerOperation::_max_loc, MPI_MAXLOC}, {SynchronizerOperation::_null, MPI_OP_NULL}}; return _operations[op]; } template MPI_Datatype inline getMPIDatatype() { return MPI_DATATYPE_NULL; } #define SPECIALIZE_MPI_DATATYPE(type, mpi_type) \ template <> MPI_Datatype inline getMPIDatatype() { return mpi_type; } #define COMMA , SPECIALIZE_MPI_DATATYPE(char, MPI_CHAR) SPECIALIZE_MPI_DATATYPE(float, MPI_FLOAT) SPECIALIZE_MPI_DATATYPE(double, MPI_DOUBLE) SPECIALIZE_MPI_DATATYPE(long double, MPI_LONG_DOUBLE) SPECIALIZE_MPI_DATATYPE(signed int, MPI_INT) - SPECIALIZE_MPI_DATATYPE(NodeType, getMPIDatatype()) SPECIALIZE_MPI_DATATYPE(unsigned int, MPI_UNSIGNED) SPECIALIZE_MPI_DATATYPE(signed long int, MPI_LONG) SPECIALIZE_MPI_DATATYPE(unsigned long int, MPI_UNSIGNED_LONG) SPECIALIZE_MPI_DATATYPE(signed long long int, MPI_LONG_LONG) SPECIALIZE_MPI_DATATYPE(unsigned long long int, MPI_UNSIGNED_LONG_LONG) SPECIALIZE_MPI_DATATYPE(SCMinMaxLoc, MPI_DOUBLE_INT) SPECIALIZE_MPI_DATATYPE(SCMinMaxLoc, MPI_FLOAT_INT) SPECIALIZE_MPI_DATATYPE(bool, MPI_CXX_BOOL) + template <> MPI_Datatype inline getMPIDatatype() { return getMPIDatatype>(); } + inline int getMPISource(int src) { if (src == _any_source) return MPI_ANY_SOURCE; return src; } decltype(auto) convertRequests(std::vector & requests) { std::vector mpi_requests(requests.size()); for (auto && request_pair : zip(requests, mpi_requests)) { auto && req = std::get<0>(request_pair); auto && mpi_req = std::get<1>(request_pair); mpi_req = dynamic_cast(req.getInternal()) .getMPIRequest(); } return mpi_requests; } } // namespace // this is ugly but shorten the code a lot #define MPIDATA \ (*reinterpret_cast(communicator_data.get())) /* -------------------------------------------------------------------------- */ /* Implementation */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ Communicator::Communicator(int & /*argc*/, char **& /*argv*/, const private_member & /*unused*/) : communicator_data(std::make_unique()) { prank = MPIDATA.rank(); psize = MPIDATA.size(); } /* -------------------------------------------------------------------------- */ template void Communicator::sendImpl(const T * buffer, Int size, Int receiver, Int tag, const CommunicationMode & mode) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Datatype type = getMPIDatatype(); switch (mode) { case CommunicationMode::_auto: MPI_Send(buffer, size, type, receiver, tag, communicator); break; case CommunicationMode::_synchronous: MPI_Ssend(buffer, size, type, receiver, tag, communicator); break; case CommunicationMode::_ready: MPI_Rsend(buffer, size, type, receiver, tag, communicator); break; } } /* -------------------------------------------------------------------------- */ template void Communicator::receiveImpl(T * buffer, Int size, Int sender, Int tag) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Status status; MPI_Datatype type = getMPIDatatype(); MPI_Recv(buffer, size, type, getMPISource(sender), tag, communicator, &status); } /* -------------------------------------------------------------------------- */ template CommunicationRequest Communicator::asyncSendImpl(const T * buffer, Int size, Int receiver, Int tag, const CommunicationMode & mode) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); auto * request = new CommunicationRequestMPI(prank, receiver); MPI_Request & req = request->getMPIRequest(); MPI_Datatype type = getMPIDatatype(); switch (mode) { case CommunicationMode::_auto: MPI_Isend(buffer, size, type, receiver, tag, communicator, &req); break; case CommunicationMode::_synchronous: MPI_Issend(buffer, size, type, receiver, tag, communicator, &req); break; case CommunicationMode::_ready: MPI_Irsend(buffer, size, type, receiver, tag, communicator, &req); break; } return std::shared_ptr(request); } /* -------------------------------------------------------------------------- */ template CommunicationRequest Communicator::asyncReceiveImpl(T * buffer, Int size, Int sender, Int tag) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); auto * request = new CommunicationRequestMPI(sender, prank); MPI_Datatype type = getMPIDatatype(); MPI_Request & req = request->getMPIRequest(); MPI_Irecv(buffer, size, type, getMPISource(sender), tag, communicator, &req); return std::shared_ptr(request); } /* -------------------------------------------------------------------------- */ template void Communicator::probe(Int sender, Int tag, CommunicationStatus & status) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Status mpi_status; MPI_Probe(getMPISource(sender), tag, communicator, &mpi_status); MPI_Datatype type = getMPIDatatype(); int count; MPI_Get_count(&mpi_status, type, &count); status.setSource(mpi_status.MPI_SOURCE); status.setTag(mpi_status.MPI_TAG); status.setSize(count); } /* -------------------------------------------------------------------------- */ template bool Communicator::asyncProbe(Int sender, Int tag, CommunicationStatus & status) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Status mpi_status; int test; MPI_Iprobe(getMPISource(sender), tag, communicator, &test, &mpi_status); if (not test) return false; MPI_Datatype type = getMPIDatatype(); int count; MPI_Get_count(&mpi_status, type, &count); status.setSource(mpi_status.MPI_SOURCE); status.setTag(mpi_status.MPI_TAG); status.setSize(count); return true; } /* -------------------------------------------------------------------------- */ bool Communicator::test(CommunicationRequest & request) const { MPI_Status status; int flag; auto & req_mpi = dynamic_cast(request.getInternal()); MPI_Request & req = req_mpi.getMPIRequest(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Test(&req, &flag, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Test."); return (flag != 0); } /* -------------------------------------------------------------------------- */ bool Communicator::testAll(std::vector & requests) const { int are_finished; auto && mpi_requests = convertRequests(requests); MPI_Testall(mpi_requests.size(), mpi_requests.data(), &are_finished, MPI_STATUSES_IGNORE); return are_finished != 0 ? true : false; } /* -------------------------------------------------------------------------- */ void Communicator::wait(CommunicationRequest & request) const { MPI_Status status; auto & req_mpi = dynamic_cast(request.getInternal()); MPI_Request & req = req_mpi.getMPIRequest(); MPI_Wait(&req, &status); } /* -------------------------------------------------------------------------- */ void Communicator::waitAll(std::vector & requests) const { auto && mpi_requests = convertRequests(requests); MPI_Waitall(mpi_requests.size(), mpi_requests.data(), MPI_STATUSES_IGNORE); } /* -------------------------------------------------------------------------- */ UInt Communicator::waitAny(std::vector & requests) const { auto && mpi_requests = convertRequests(requests); int pos; MPI_Waitany(mpi_requests.size(), mpi_requests.data(), &pos, MPI_STATUSES_IGNORE); if (pos != MPI_UNDEFINED) { return pos; } else { return UInt(-1); } } /* -------------------------------------------------------------------------- */ void Communicator::barrier() const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Barrier(communicator); } /* -------------------------------------------------------------------------- */ CommunicationRequest Communicator::asyncBarrier() const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); auto * request = new CommunicationRequestMPI(0, 0); MPI_Request & req = request->getMPIRequest(); MPI_Ibarrier(communicator, &req); return std::shared_ptr(request); } /* -------------------------------------------------------------------------- */ template void Communicator::reduceImpl(T * values, int nb_values, const SynchronizerOperation & op, int root) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Datatype type = getMPIDatatype(); MPI_Reduce(MPI_IN_PLACE, values, nb_values, type, getMPISynchronizerOperation(op), root, communicator); } /* -------------------------------------------------------------------------- */ template void Communicator::allReduceImpl(T * values, int nb_values, const SynchronizerOperation & op) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Datatype type = getMPIDatatype(); MPI_Allreduce(MPI_IN_PLACE, values, nb_values, type, getMPISynchronizerOperation(op), communicator); } /* -------------------------------------------------------------------------- */ template void Communicator::allGatherImpl(T * values, int nb_values) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Datatype type = getMPIDatatype(); MPI_Allgather(MPI_IN_PLACE, nb_values, type, values, nb_values, type, communicator); } /* -------------------------------------------------------------------------- */ template void Communicator::allGatherVImpl(T * values, int * nb_values) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); std::vector displs(psize); displs[0] = 0; for (int i = 1; i < psize; ++i) { displs[i] = displs[i - 1] + nb_values[i - 1]; } MPI_Datatype type = getMPIDatatype(); MPI_Allgatherv(MPI_IN_PLACE, *nb_values, type, values, nb_values, displs.data(), type, communicator); } /* -------------------------------------------------------------------------- */ template void Communicator::gatherImpl(T * values, int nb_values, int root) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); T *send_buf = nullptr, *recv_buf = nullptr; if (prank == root) { send_buf = (T *)MPI_IN_PLACE; recv_buf = values; } else { send_buf = values; } MPI_Datatype type = getMPIDatatype(); MPI_Gather(send_buf, nb_values, type, recv_buf, nb_values, type, root, communicator); } /* -------------------------------------------------------------------------- */ template void Communicator::gatherImpl(T * values, int nb_values, T * gathered, int nb_gathered) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); T * send_buf = values; T * recv_buf = gathered; if (nb_gathered == 0) nb_gathered = nb_values; MPI_Datatype type = getMPIDatatype(); MPI_Gather(send_buf, nb_values, type, recv_buf, nb_gathered, type, this->prank, communicator); } /* -------------------------------------------------------------------------- */ template void Communicator::gatherVImpl(T * values, int * nb_values, int root) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); int * displs = nullptr; if (prank == root) { displs = new int[psize]; displs[0] = 0; for (int i = 1; i < psize; ++i) { displs[i] = displs[i - 1] + nb_values[i - 1]; } } T *send_buf = nullptr, *recv_buf = nullptr; if (prank == root) { send_buf = (T *)MPI_IN_PLACE; recv_buf = values; } else send_buf = values; MPI_Datatype type = getMPIDatatype(); MPI_Gatherv(send_buf, *nb_values, type, recv_buf, nb_values, displs, type, root, communicator); if (prank == root) { delete[] displs; } } /* -------------------------------------------------------------------------- */ template void Communicator::broadcastImpl(T * values, int nb_values, int root) const { MPI_Comm communicator = MPIDATA.getMPICommunicator(); MPI_Datatype type = getMPIDatatype(); MPI_Bcast(values, nb_values, type, root, communicator); } /* -------------------------------------------------------------------------- */ int Communicator::getMaxTag() const { return MPIDATA.getMaxTag(); } int Communicator::getMinTag() const { return 0; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/grid_synchronizer.cc b/src/synchronizer/grid_synchronizer.cc index 1c14f23a8..a9245572a 100644 --- a/src/synchronizer/grid_synchronizer.cc +++ b/src/synchronizer/grid_synchronizer.cc @@ -1,475 +1,475 @@ /** * @file grid_synchronizer.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Oct 03 2011 * @date last modification: Tue Nov 07 2017 * * @brief implementation of the grid synchronizer * * @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 "grid_synchronizer.hh" #include "aka_grid_dynamic.hh" #include "communicator.hh" #include "fe_engine.hh" #include "integration_point.hh" #include "mesh.hh" #include "mesh_io.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template void GridSynchronizer::createGridSynchronizer(const SpatialGrid & grid) { AKANTU_DEBUG_IN(); const Communicator & comm = this->mesh.getCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); if (nb_proc == 1) return; UInt spatial_dimension = this->mesh.getSpatialDimension(); BBox my_bounding_box(spatial_dimension); const auto & lower = grid.getLowerBounds(); const auto & upper = grid.getUpperBounds(); const auto & spacing = grid.getSpacing(); my_bounding_box.getLowerBounds() = lower - spacing; my_bounding_box.getUpperBounds() = upper + spacing; AKANTU_DEBUG_INFO( "Exchange of bounding box to detect the overlapping regions."); auto && bboxes = my_bounding_box.allGather(comm); std::vector intersects_proc(nb_proc); std::fill(intersects_proc.begin(), intersects_proc.end(), true); Matrix first_cells(spatial_dimension, nb_proc); Matrix last_cells(spatial_dimension, nb_proc); std::map> element_per_proc; // check the overlapping between my box and the one from other processors for (UInt p = 0; p < nb_proc; ++p) { if (p == my_rank) continue; const auto & proc_bounding_box = bboxes[p]; auto intersection = my_bounding_box.intersection(proc_bounding_box); Vector first_cell_p = first_cells(p); Vector last_cell_p = last_cells(p); intersects_proc[p] = intersection; if(intersects_proc[p]) { for (UInt s = 0; s < spatial_dimension; ++s) { first_cell_p(s) = grid.getCellID(intersection.getLowerBounds()(s), s); last_cell_p(s) = grid.getCellID(intersection.getUpperBounds()(s), s); } } // create the list of cells in the overlapping using CellID = typename SpatialGrid::CellID; std::vector cell_ids; if (intersects_proc[p]) { AKANTU_DEBUG_INFO("I intersects with processor " << p); CellID cell_id(spatial_dimension); // for (UInt i = 0; i < spatial_dimension; ++i) { // if(first_cell_p[i] != 0) --first_cell_p[i]; // if(last_cell_p[i] != 0) ++last_cell_p[i]; // } for (Int fd = first_cell_p(0); fd <= last_cell_p(0); ++fd) { cell_id.setID(0, fd); if (spatial_dimension == 1) { cell_ids.push_back(cell_id); } else { for (Int sd = first_cell_p(1); sd <= last_cell_p(1); ++sd) { cell_id.setID(1, sd); if (spatial_dimension == 2) { cell_ids.push_back(cell_id); } else { for (Int ld = first_cell_p(2); ld <= last_cell_p(2); ++ld) { cell_id.setID(2, ld); cell_ids.push_back(cell_id); } } } } } // get the list of elements in the cells of the overlapping std::set to_send; for (auto & cur_cell_id : cell_ids) { auto & cell = grid.getCell(cur_cell_id); for (auto & element : cell) { to_send.insert(element); } } AKANTU_DEBUG_INFO("I have prepared " << to_send.size() << " elements to send to processor " << p); auto & scheme = this->getCommunications().createSendScheme(p); std::stringstream sstr; sstr << "element_per_proc_" << p; element_per_proc.emplace( std::piecewise_construct, std::forward_as_tuple(p), std::forward_as_tuple(sstr.str(), id, memory_id)); ElementTypeMapArray & elempproc = element_per_proc[p]; for (auto elem : to_send) { ElementType type = elem.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); // /!\ this part must be slow due to the access in the // ElementTypeMapArray if (!elempproc.exists(type, _not_ghost)) elempproc.alloc(0, nb_nodes_per_element, type, _not_ghost); Vector global_connect(nb_nodes_per_element); Vector local_connect = mesh.getConnectivity(type).begin( nb_nodes_per_element)[elem.element]; for (UInt i = 0; i < nb_nodes_per_element; ++i) { global_connect(i) = mesh.getNodeGlobalId(local_connect(i)); AKANTU_DEBUG_ASSERT( global_connect(i) < mesh.getNbGlobalNodes(), "This global node send in the connectivity does not seem correct " << global_connect(i) << " corresponding to " << local_connect(i) << " from element " << elem.element); } elempproc(type).push_back(global_connect); scheme.push_back(elem); } } } AKANTU_DEBUG_INFO("I have finished to compute intersection," << " no it's time to communicate with my neighbors"); /** * Sending loop, sends the connectivity asynchronously to all concerned proc */ std::vector isend_requests; Tensor3 space(2, _max_element_type, nb_proc); for (UInt p = 0; p < nb_proc; ++p) { if (p == my_rank) continue; if (not intersects_proc[p]) continue; Matrix info_proc = space(p); auto & elempproc = element_per_proc[p]; UInt count = 0; for (auto type : elempproc.elementTypes(_all_dimensions, _not_ghost)) { Array & conn = elempproc(type, _not_ghost); Vector info = info_proc((UInt)type); info[0] = (UInt)type; info[1] = conn.size() * conn.getNbComponent(); AKANTU_DEBUG_INFO( "I have " << conn.size() << " elements of type " << type << " to send to processor " << p << " (communication tag : " << Tag::genTag(my_rank, count, DATA_TAG) << ")"); isend_requests.push_back( comm.asyncSend(info, p, Tag::genTag(my_rank, count, SIZE_TAG))); if (info[1] != 0) isend_requests.push_back(comm.asyncSend( conn, p, Tag::genTag(my_rank, count, DATA_TAG))); ++count; } Vector info = info_proc((UInt)_not_defined); info[0] = (UInt)_not_defined; info[1] = 0; isend_requests.push_back( comm.asyncSend(info, p, Tag::genTag(my_rank, count, SIZE_TAG))); } /** * Receives the connectivity and store them in the ghosts elements */ MeshAccessor mesh_accessor(mesh); auto & global_nodes_ids = mesh_accessor.getNodesGlobalIds(); - auto & nodes_type = mesh_accessor.getNodesType(); + auto & nodes_type = mesh_accessor.getNodesFlags(); std::vector isend_nodes_requests; Vector nb_nodes_to_recv(nb_proc); UInt nb_total_nodes_to_recv = 0; UInt nb_current_nodes = global_nodes_ids.size(); NewNodesEvent new_nodes; NewElementsEvent new_elements; std::map> ask_nodes_per_proc; for (UInt p = 0; p < nb_proc; ++p) { nb_nodes_to_recv(p) = 0; if (p == my_rank) continue; if (!intersects_proc[p]) continue; auto & scheme = this->getCommunications().createRecvScheme(p); ask_nodes_per_proc.emplace(std::piecewise_construct, std::forward_as_tuple(p), std::forward_as_tuple(0)); auto & ask_nodes = ask_nodes_per_proc[p]; UInt count = 0; ElementType type = _not_defined; do { Vector info(2); comm.receive(info, p, Tag::genTag(p, count, SIZE_TAG)); type = (ElementType)info[0]; if (type == _not_defined) break; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = info[1] / nb_nodes_per_element; Array tmp_conn(nb_element, nb_nodes_per_element); tmp_conn.clear(); if (info[1] != 0) comm.receive(tmp_conn, p, Tag::genTag(p, count, DATA_TAG)); AKANTU_DEBUG_INFO("I will receive " << nb_element << " elements of type " << ElementType(info[0]) << " from processor " << p << " (communication tag : " << Tag::genTag(p, count, DATA_TAG) << ")"); auto & ghost_connectivity = mesh_accessor.getConnectivity(type, _ghost); auto & ghost_counter = mesh_accessor.getGhostsCounters(type, _ghost); UInt nb_ghost_element = ghost_connectivity.size(); Element element{type, 0, _ghost}; Vector conn(nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el) { UInt nb_node_to_ask_for_elem = 0; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt gn = tmp_conn(el, n); UInt ln = global_nodes_ids.find(gn); AKANTU_DEBUG_ASSERT(gn < mesh.getNbGlobalNodes(), "This global node seems not correct " << gn << " from element " << el << " node " << n); if (ln == UInt(-1)) { global_nodes_ids.push_back(gn); - nodes_type.push_back(_nt_pure_ghost); // pure ghost node + nodes_type.push_back(NodeFlag::_pure_ghost); // pure ghost node ln = nb_current_nodes; new_nodes.getList().push_back(ln); ++nb_current_nodes; ask_nodes.push_back(gn); ++nb_node_to_ask_for_elem; } conn[n] = ln; } // all the nodes are already known locally, the element should // already exists auto c = UInt(-1); if (nb_node_to_ask_for_elem == 0) { c = ghost_connectivity.find(conn); element.element = c; } if (c == UInt(-1)) { element.element = nb_ghost_element; ++nb_ghost_element; ghost_connectivity.push_back(conn); ghost_counter.push_back(1); new_elements.getList().push_back(element); } else { ++ghost_counter(c); } scheme.push_back(element); } count++; } while (type != _not_defined); AKANTU_DEBUG_INFO("I have " << ask_nodes.size() << " missing nodes for elements coming from processor " << p << " (communication tag : " << Tag::genTag(my_rank, 0, ASK_NODES_TAG) << ")"); ask_nodes.push_back(UInt(-1)); isend_nodes_requests.push_back( comm.asyncSend(ask_nodes, p, Tag::genTag(my_rank, 0, ASK_NODES_TAG))); nb_nodes_to_recv(p) = ask_nodes.size() - 1; nb_total_nodes_to_recv += nb_nodes_to_recv(p); } comm.waitAll(isend_requests); comm.freeCommunicationRequest(isend_requests); /** * Sends requested nodes to proc */ auto & nodes = const_cast &>(mesh.getNodes()); UInt nb_nodes = nodes.size(); std::vector isend_coordinates_requests; std::map> nodes_to_send_per_proc; for (UInt p = 0; p < nb_proc; ++p) { if (p == my_rank || !intersects_proc[p]) continue; Array asked_nodes; CommunicationStatus status; AKANTU_DEBUG_INFO("Waiting list of nodes to send to processor " << p << "(communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")"); comm.probe(p, Tag::genTag(p, 0, ASK_NODES_TAG), status); UInt nb_nodes_to_send = status.size(); asked_nodes.resize(nb_nodes_to_send); AKANTU_DEBUG_INFO("I have " << nb_nodes_to_send - 1 << " nodes to send to processor " << p << " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")"); AKANTU_DEBUG_INFO("Getting list of nodes to send to processor " << p << " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")"); comm.receive(asked_nodes, p, Tag::genTag(p, 0, ASK_NODES_TAG)); nb_nodes_to_send--; asked_nodes.resize(nb_nodes_to_send); nodes_to_send_per_proc.emplace(std::piecewise_construct, std::forward_as_tuple(p), std::forward_as_tuple(0, spatial_dimension)); auto & nodes_to_send = nodes_to_send_per_proc[p]; auto node_it = nodes.begin(spatial_dimension); for (UInt n = 0; n < nb_nodes_to_send; ++n) { UInt ln = global_nodes_ids.find(asked_nodes(n)); AKANTU_DEBUG_ASSERT(ln != UInt(-1), "The node [" << asked_nodes(n) << "] requested by proc " << p << " was not found locally!"); nodes_to_send.push_back(node_it + ln); } if (nb_nodes_to_send != 0) { AKANTU_DEBUG_INFO("Sending the " << nb_nodes_to_send << " nodes to processor " << p << " (communication tag : " << Tag::genTag(p, 0, SEND_NODES_TAG) << ")"); isend_coordinates_requests.push_back(comm.asyncSend( nodes_to_send, p, Tag::genTag(my_rank, 0, SEND_NODES_TAG))); } #if not defined(AKANTU_NDEBUG) else { AKANTU_DEBUG_INFO("No nodes to send to processor " << p); } #endif } comm.waitAll(isend_nodes_requests); comm.freeCommunicationRequest(isend_nodes_requests); nodes.resize(nb_total_nodes_to_recv + nb_nodes); for (UInt p = 0; p < nb_proc; ++p) { if ((p != my_rank) && (nb_nodes_to_recv(p) > 0)) { AKANTU_DEBUG_INFO("Receiving the " << nb_nodes_to_recv(p) << " nodes from processor " << p << " (communication tag : " << Tag::genTag(p, 0, SEND_NODES_TAG) << ")"); Vector nodes_to_recv(nodes.storage() + nb_nodes * spatial_dimension, nb_nodes_to_recv(p) * spatial_dimension); comm.receive(nodes_to_recv, p, Tag::genTag(p, 0, SEND_NODES_TAG)); nb_nodes += nb_nodes_to_recv(p); } #if not defined(AKANTU_NDEBUG) else { if (p != my_rank) { AKANTU_DEBUG_INFO("No nodes to receive from processor " << p); } } #endif } comm.waitAll(isend_coordinates_requests); comm.freeCommunicationRequest(isend_coordinates_requests); mesh.sendEvent(new_nodes); mesh.sendEvent(new_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void GridSynchronizer::createGridSynchronizer( const SpatialGrid & grid); template void GridSynchronizer::createGridSynchronizer( const SpatialGrid & grid); } // namespace akantu diff --git a/src/synchronizer/node_info_per_processor.cc b/src/synchronizer/node_info_per_processor.cc index a450c7761..70759d691 100644 --- a/src/synchronizer/node_info_per_processor.cc +++ b/src/synchronizer/node_info_per_processor.cc @@ -1,499 +1,516 @@ /** * @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); } } 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(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NodeInfoPerProc::fillNodesType() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); - Array & nodes_type = this->getNodesType(); + auto & nodes_flags = this->getNodesFlags(); 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; + for (auto gt : ghost_types) { 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 (auto && type : + mesh.elementTypes(_all_dimensions, gt, _ek_not_defined)) { + const auto & connectivity = mesh.getConnectivity(type, gt); - 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) { + for (auto & conn : + make_view(connectivity, connectivity.getNbComponent())) { + for (UInt n = 0; n < conn.size(); ++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; + nodes_flags(i) = NodeFlag::_normal; else if (nodes_set(i) == GHOST_SET) - nodes_type(i) = _nt_pure_ghost; + nodes_flags(i) = NodeFlag::_pure_ghost; else if (nodes_set(i) == (GHOST_SET + NORMAL_SET)) - nodes_type(i) = _nt_master; + nodes_flags(i) = NodeFlag::_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)); + for (auto node : arange(mesh.getNbNodes())) { + if (mesh.isSlaveNode(node)) { + recv_array_per_proc[mesh.getNodePrank(node)].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(); } -/* -------------------------------------------------------------------------- */ -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ +/* -------------------------------------------------------------------------- + */ -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ 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); + std::vector> nodes_flags_per_proc(nb_proc); + std::vector> nodes_prank_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)); + nodes_flags_per_proc[p].resize(nb_nodes_per_proc(p)); + nodes_prank_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]; + auto & nodes_flags = nodes_flags_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)); + comm.receive(nodes_flags, p, Tag::genTag(p, 0, Tag::_NODES_TYPE)); } else { - nodes_types.copy(this->getNodesType()); + nodes_flags.copy(this->getNodesFlags()); } // 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) { + if ((nodes_flags(local_node) & NodeFlag::_shared_mask) == NodeFlag::_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_prank_per_proc[proc](node) = 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, + comm.asyncSend(nodes_flags_per_proc[p], p, Tag::genTag(this->rank, 0, Tag::_NODES_TYPE))); - auto & nodes_to_send = nodes_to_send_per_proc[p]; + requests_send_type.push_back( + comm.asyncSend(nodes_prank_per_proc[p], p, + Tag::genTag(this->rank, 2, Tag::_NODES_TYPE))); - /// push back an element to avoid a send of size 0 - nodes_to_send.push_back(-1); + auto & nodes_to_send = nodes_to_send_per_proc[p]; 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->getNodesFlags().copy(nodes_flags_per_proc[p]); + for(auto && data : enumerate(nodes_prank_per_proc[p])) { + auto node = std::get<0>(data); + if(mesh.isSlaveNode(node)) { + this->setNodePrank(node, std::get<1>(data)); + } + } 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(); } -/* -------------------------------------------------------------------------- */ -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ +/* -------------------------------------------------------------------------- + */ -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ 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(); + auto & nodes_types = this->getNodesFlags(); 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)); + Array nodes_prank(nodes_types.size()); + comm.receive(nodes_prank, root, Tag::genTag(root, 2, Tag::_NODES_TYPE)); + for(auto && data : enumerate(nodes_prank)) { + auto node = std::get<0>(data); + if(mesh.isSlaveNode(node)) { + this->setNodePrank(node, std::get<1>(data)); + } + } + 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); + if(nodes_master_info.size()>0) + comm.receive(nodes_master_info, root, Tag::genTag(root, 1, Tag::_NODES_TYPE)); 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 diff --git a/src/synchronizer/node_synchronizer.cc b/src/synchronizer/node_synchronizer.cc index edb4c7a28..18f461357 100644 --- a/src/synchronizer/node_synchronizer.cc +++ b/src/synchronizer/node_synchronizer.cc @@ -1,177 +1,176 @@ /** * @file node_synchronizer.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 15 2017 * * @brief Implementation of the node synchronizer * * @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 "node_synchronizer.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NodeSynchronizer::NodeSynchronizer(Mesh & mesh, const ID & id, MemoryID memory_id, const bool register_to_event_manager, EventHandlerPriority event_priority) : SynchronizerImpl(mesh.getCommunicator(), id, memory_id), mesh(mesh) { AKANTU_DEBUG_IN(); if (register_to_event_manager) { this->mesh.registerEventHandler(*this, event_priority); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NodeSynchronizer::~NodeSynchronizer() = default; /* -------------------------------------------------------------------------- */ void NodeSynchronizer::onNodesAdded(const Array & /*nodes_list*/, const NewNodesEvent &) { std::map> nodes_per_proc; // recreates fully the schemes due to changes of global ids // \TODO add an event to handle global id changes for (auto && data : communications.iterateSchemes(_recv)) { auto & scheme = data.second; scheme.resize(0); } for (auto && local_id : arange(mesh.getNbNodes())) { - auto type = mesh.getNodeType(local_id); - if (type < 0) + if (not mesh.isSlaveNode(local_id)) continue; // local, master or pure ghost auto global_id = mesh.getNodeGlobalId(local_id); - auto proc = UInt(type); + auto proc = mesh.getNodePrank(local_id); nodes_per_proc[proc].push_back(global_id); auto & scheme = communications.createScheme(proc, _recv); scheme.push_back(local_id); } std::vector send_requests; for (auto && pair : communications.iterateSchemes(_recv)) { auto proc = pair.first; // if proc not in nodes_per_proc this should insert an empty array to send send_requests.push_back(communicator.asyncSend( nodes_per_proc[proc], proc, Tag::genTag(rank, proc, 0xcafe))); } for (auto && data : communications.iterateSchemes(_send)) { auto proc = data.first; auto & scheme = data.second; CommunicationStatus status; auto tag = Tag::genTag(proc, rank, 0xcafe); communicator.probe(proc, tag, status); scheme.resize(status.size()); communicator.receive(scheme, proc, tag); std::transform(scheme.begin(), scheme.end(), scheme.begin(), [&](auto & gnode) { return mesh.getNodeLocalId(gnode); }); } // communicator.receiveAnyNumber( // send_requests, // [&](auto && proc, auto && nodes) { // auto & scheme = communications.createScheme(proc, _send); // scheme.resize(nodes.size()); // for (auto && data : enumerate(nodes)) { // auto global_id = std::get<1>(data); // auto local_id = mesh.getNodeLocalId(global_id); // AKANTU_DEBUG_ASSERT(local_id != UInt(-1), // "The global node " << global_id // << "is not known on rank " // << rank); // scheme[std::get<0>(data)] = local_id; // } // }, // Tag::genTag(rank, count, 0xcafe)); // ++count; communicator.waitAll(send_requests); communicator.freeCommunicationRequest(send_requests); } /* -------------------------------------------------------------------------- */ UInt NodeSynchronizer::sanityCheckDataSize(const Array & nodes, const SynchronizationTag & tag, bool from_comm_desc) const { UInt size = SynchronizerImpl::sanityCheckDataSize(nodes, tag, from_comm_desc); // positions size += mesh.getSpatialDimension() * sizeof(Real) * nodes.size(); return size; } /* -------------------------------------------------------------------------- */ void NodeSynchronizer::packSanityCheckData( CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & /*tag*/) const { auto dim = mesh.getSpatialDimension(); for (auto && node : nodes) { buffer << Vector(mesh.getNodes().begin(dim)[node]); } } /* -------------------------------------------------------------------------- */ void NodeSynchronizer::unpackSanityCheckData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag, UInt proc, UInt rank) const { auto dim = mesh.getSpatialDimension(); // std::set skip_conn_tags{_gst_smmc_facets_conn, // _gst_giu_global_conn}; // bool is_skip_tag_conn = skip_conn_tags.find(tag) != skip_conn_tags.end(); for (auto && node : nodes) { Vector pos_remote(dim); buffer >> pos_remote; Vector pos(mesh.getNodes().begin(dim)[node]); auto dist = pos_remote.distance(pos); if (not Math::are_float_equal(dist, 0.)) { AKANTU_EXCEPTION("Unpacking an unknown value for the node " << node << "(position " << pos << " != buffer " << pos_remote << ") [" << dist << "] - tag: " << tag << " comm from " << proc << " to " << rank); } } } /* -------------------------------------------------------------------------- */ } // namespace akantu