diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 0ae09052d..bf9cdb26a 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,477 +1,477 @@ /** * @file aka_common.hh * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Thu Jan 21 2016 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ } // namespace akantu #include "aka_element_classes_info.hh" namespace akantu { #define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) #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))) #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) \ } \ namespace std { \ inline string to_string(const ::akantu::type_name & type) { \ static unordered_map<::akantu::type_name, 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[type]; \ + return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << std::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[str]; \ + type = convert.at(str); \ return stream; \ } /// small help to use names for directions enum SpacialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_tranfser_model) \ (structural_mechanics_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) /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent _nls_lumped, ///< Case of lumped mass or equivalent matrix _nls_auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; /// Define the node/dof type enum NodeType : Int { _nt_pure_gost = -3, _nt_master = -2, _nt_normal = -1 }; /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_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_capacity, ///< synchronization of the nodal heat capacity _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, _ghost, _casper // not used but a real cute ghost }; /* -------------------------------------------------------------------------- */ 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/parser/input_file_parser.hh b/src/io/parser/input_file_parser.hh index dab04c652..201f9a9b9 100644 --- a/src/io/parser/input_file_parser.hh +++ b/src/io/parser/input_file_parser.hh @@ -1,263 +1,270 @@ /** * @file input_file_parser.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Tue May 19 2015 * * @brief Grammar definition for the input files * * @section LICENSE * * Copyright (©) 2014, 2015 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 . * */ /* -------------------------------------------------------------------------- */ // Boost /* -------------------------------------------------------------------------- */ #include -#include +#include #include #include #include -#include +#include #include #ifndef __AKANTU_INPUT_FILE_PARSER_HH__ #define __AKANTU_INPUT_FILE_PARSER_HH__ namespace spirit = boost::spirit; namespace qi = boost::spirit::qi; namespace lbs = boost::spirit::qi::labels; namespace ascii = boost::spirit::ascii; namespace phx = boost::phoenix; namespace akantu { namespace parser { -struct error_handler_ { - template struct result { - using type = void; + struct error_handler_ { + template struct result { + using type = void; + }; + + template + void operator()(qi::info const & what, Iterator err_pos, Iterator /*first*/, + Iterator /*last*/) const { + spirit::classic::file_position pos = err_pos.get_position(); + + AKANTU_EXCEPTION("Parse error [ " + << "Expecting " << what << " instead of \"" << *err_pos + << "\" ]" + << " in file " << pos.file << " line " << pos.line + << " column " << pos.column << std::endl + << "'" << err_pos.get_currentline() << "'" << std::endl + << std::setw(pos.column) << " " + << "^- here"); + } + + private: }; - template - void operator()(qi::info const & what, Iterator err_pos, - __attribute__((unused)) Iterator first, - __attribute__((unused)) Iterator last) const { - spirit::classic::file_position pos = err_pos.get_position(); - - AKANTU_EXCEPTION("Parse error [ " - << "Expecting " << what << " instead of \"" << *err_pos - << "\" ]" - << " in file " << pos.file << " line " << pos.line - << " column " << pos.column << std::endl - << "'" << err_pos.get_currentline() << "'" << std::endl - << std::setw(pos.column) << " " - << "^- here"); - } + static ParserSection & + create_subsection(const ParserType & type, const boost::optional & opt_name, + const boost::optional & opt_option, + ParserSection & sect) { + std::string option = ""; + if (opt_option) + option = *opt_option; -private: -}; + static size_t id = 12; + std::string name = "anonymous_" + std::to_string(id++); + if (opt_name) + name = *opt_name; -static ParserSection & -create_subsection(const ParserType & type, const std::string & name, - const boost::optional & option, - ParserSection & sect) { - std::string opt = ""; - if (option) - opt = *option; - ParserSection sect_tmp(name, type, opt, sect); - return sect.addSubSection(sect_tmp); -} + ParserSection sect_tmp(name, type, option, sect); + return sect.addSubSection(sect_tmp); + } -template -static bool create_parameter(boost::iterator_range & rng, - std::string & value, ParserSection & sect) { - try { - std::string name(rng.begin(), rng.end()); - name = trim(name); - spirit::classic::file_position pos = rng.begin().get_position(); - - ParserParameter param_tmp(name, value, sect); - param_tmp.setDebugInfo(pos.file, pos.line, pos.column); - sect.addParameter(param_tmp); - } catch (debug::Exception & e) { - return false; + template + static bool create_parameter(boost::iterator_range & rng, + std::string & value, ParserSection & sect) { + try { + std::string name(rng.begin(), rng.end()); + name = trim(name); + spirit::classic::file_position pos = rng.begin().get_position(); + + ParserParameter param_tmp(name, value, sect); + param_tmp.setDebugInfo(pos.file, pos.line, pos.column); + sect.addParameter(param_tmp); + } catch (debug::Exception & e) { + return false; + } + return true; } - return true; -} -static std::string concatenate(const std::string & t1, const std::string & t2) { - return (t1 + t2); -} + static std::string concatenate(const std::string & t1, + const std::string & t2) { + return (t1 + t2); + } -/* ---------------------------------------------------------------------- */ -/* Grammars definitions */ -/* ---------------------------------------------------------------------- */ -template -struct InputFileGrammar - : qi::grammar::type> { - InputFileGrammar(ParserSection * sect) - : InputFileGrammar::base_type(start, "input_file_grammar"), - parent_section(sect) { - - /* clang-format off */ + /* ---------------------------------------------------------------------- */ + /* Grammars definitions */ + /* ---------------------------------------------------------------------- */ + template + struct InputFileGrammar + : qi::grammar::type> { + InputFileGrammar(ParserSection * sect) + : InputFileGrammar::base_type(start, "input_file_grammar"), + parent_section(sect) { + + /* clang-format off */ start = mini_section(parent_section) ; mini_section = *( entry (lbs::_r1) | section(lbs::_r1) ) ; entry = ( - qi::raw[key] + qi::raw[key] >> '=' - > value - ) [ lbs::_pass = phx::bind(&create_parameter, - lbs::_1, - lbs::_2, - *lbs::_r1) ] + > value + ) [ lbs::_pass = phx::bind(&create_parameter, + lbs::_1, + lbs::_2, + *lbs::_r1) ] ; section = ( - qi::no_case[section_type] + qi::no_case[section_type] > qi::lexeme [ - section_name + -section_name > -section_option ] ) [ lbs::_a = &phx::bind(&create_subsection, lbs::_1, phx::at_c<0>(lbs::_2), phx::at_c<1>(lbs::_2), *lbs::_r1) ] > '[' > mini_section(lbs::_a) > ']' ; section_name = qi::char_("a-zA-Z_") >> *qi::char_("a-zA-Z_0-9") ; section_option = (+ascii::space >> section_name) [ lbs::_val = lbs::_2 ] ; key = qi::char_("a-zA-Z_") >> *qi::char_("a-zA-Z_0-9") ; value = ( mono_line_value [ lbs::_a = phx::bind(&concatenate, lbs::_a, lbs::_1) ] > *( '\\' > mono_line_value [ lbs::_a = phx::bind(&concatenate, lbs::_a, lbs::_1) ] ) ) [ lbs::_val = lbs::_a ] ; mono_line_value = qi::lexeme [ +(qi::char_ - (qi::char_('=') | spirit::eol | '#' | ';' | '\\')) ] ; skipper = ascii::space | "#" >> *(qi::char_ - spirit::eol) ; /* clang-format on */ #define AKANTU_SECTION_TYPE_ADD(r, data, elem) \ - (BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(ParserType::_,elem)) + (BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(ParserType::_, elem)) - section_type.add BOOST_PP_SEQ_FOR_EACH(AKANTU_SECTION_TYPE_ADD, _, - AKANTU_SECTION_TYPES); + section_type.add BOOST_PP_SEQ_FOR_EACH(AKANTU_SECTION_TYPE_ADD, _, + AKANTU_SECTION_TYPES); #undef AKANTU_SECTION_TYPE_ADD #if !defined(AKANTU_NDEBUG) && defined(AKANTU_CORE_CXX_11) - phx::function const error_handler = error_handler_(); - qi::on_error(start, - error_handler(lbs::_4, lbs::_3, lbs::_1, lbs::_2)); + phx::function const error_handler = error_handler_(); + qi::on_error(start, + error_handler(lbs::_4, lbs::_3, lbs::_1, lbs::_2)); #endif - section.name("section"); - section_name.name("section-name"); - section_option.name("section-option"); - mini_section.name("section-content"); - entry.name("parameter"); - key.name("parameter-name"); - value.name("parameter-value"); - section_type.name("section-types-list"); - mono_line_value.name("mono-line-value"); + section.name("section"); + section_name.name("section-name"); + section_option.name("section-option"); + mini_section.name("section-content"); + entry.name("parameter"); + key.name("parameter-name"); + value.name("parameter-value"); + section_type.name("section-types-list"); + mono_line_value.name("mono-line-value"); #if !defined AKANTU_NDEBUG - if (AKANTU_DEBUG_TEST(dblDebug)) { - // qi::debug(section); - qi::debug(section_name); - qi::debug(section_option); - // qi::debug(mini_section); - // qi::debug(entry); - qi::debug(key); - qi::debug(value); - qi::debug(mono_line_value); - } + if (AKANTU_DEBUG_TEST(dblDebug)) { + // qi::debug(section); + qi::debug(section_name); + qi::debug(section_option); + // qi::debug(mini_section); + // qi::debug(entry); + qi::debug(key); + qi::debug(value); + qi::debug(mono_line_value); + } #endif - } + } - const std::string & getErrorMessage() const { return error_message; }; + const std::string & getErrorMessage() const { return error_message; }; - using skipper_type = typename Skipper::type; - skipper_type skipper; + using skipper_type = typename Skipper::type; + skipper_type skipper; -private: - std::string error_message; + private: + std::string error_message; - qi::rule mini_section; - qi::rule, - skipper_type> section; - qi::rule start; - qi::rule section_name; - qi::rule section_option; - qi::rule entry; - qi::rule key; - qi::rule, skipper_type> - value; - qi::rule mono_line_value; + qi::rule mini_section; + qi::rule, + skipper_type> + section; + qi::rule start; + qi::rule section_name; + qi::rule section_option; + qi::rule entry; + qi::rule key; + qi::rule, skipper_type> + value; + qi::rule mono_line_value; - qi::symbols section_type; + qi::symbols section_type; - ParserSection * parent_section; -}; + ParserSection * parent_section; + }; } } // akantu #endif /* __AKANTU_INPUT_FILE_PARSER_HH__ */ diff --git a/src/io/parser/parameter_registry_tmpl.hh b/src/io/parser/parameter_registry_tmpl.hh index f015d5db9..3a96f7ad6 100644 --- a/src/io/parser/parameter_registry_tmpl.hh +++ b/src/io/parser/parameter_registry_tmpl.hh @@ -1,382 +1,391 @@ /** * @file parameter_registry_tmpl.hh * * @author Nicolas Richart * * @date creation: Thu Aug 09 2012 * @date last modification: Fri Mar 27 2015 * * @brief implementation of the templated part of ParameterRegistry class and * the derivated ones * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_error.hh" +#include "aka_iterators.hh" #include "parameter_registry.hh" #include "parser.hh" -#include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARAMETER_REGISTRY_TMPL_HH__ #define __AKANTU_PARAMETER_REGISTRY_TMPL_HH__ namespace akantu { namespace debug { class ParameterException : public Exception { public: ParameterException(const std::string & name, const std::string & message) : Exception(message), name(name) {} const std::string & name; }; class ParameterUnexistingException : public ParameterException { public: ParameterUnexistingException(const std::string & name, const ParameterRegistry & registery) : ParameterException( name, "Parameter " + name + " does not exists in this scope") { auto && params = registery.listParameters(); this->_info = std::accumulate(params.begin(), params.end(), this->_info + "\n Possible parameters are: ", [](auto && str, auto && param) { static auto first = true; auto && ret = str + (first ? " " : ", ") + param; first = false; return ret; }); } }; class ParameterAccessRightException : public ParameterException { public: ParameterAccessRightException(const std::string & name, const std::string & perm) : ParameterException(name, "Parameter " + name + " is not " + perm) {} }; class ParameterWrongTypeException : public ParameterException { public: ParameterWrongTypeException(const std::string & name, const std::type_info & wrong_type, const std::type_info & type) : ParameterException(name, "Parameter " + name + " type error, cannot convert " + debug::demangle(type.name()) + " to " + debug::demangle(wrong_type.name())) {} }; } // namespace debug /* -------------------------------------------------------------------------- */ template const ParameterTyped & Parameter::getParameterTyped() const { try { const auto & tmp = dynamic_cast &>(*this); return tmp; } catch (std::bad_cast &) { AKANTU_CUSTOM_EXCEPTION( debug::ParameterWrongTypeException(name, typeid(T), this->type())); } } /* -------------------------------------------------------------------------- */ template ParameterTyped & Parameter::getParameterTyped() { try { auto & tmp = dynamic_cast &>(*this); return tmp; } catch (std::bad_cast &) { AKANTU_CUSTOM_EXCEPTION( debug::ParameterWrongTypeException(name, typeid(T), this->type())); } } /* ------------------------------------------------------------------------ */ template void Parameter::set(const V & value) { if (!(isWritable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "writable")); ParameterTyped & typed_param = getParameterTyped(); typed_param.setTyped(value); } /* ------------------------------------------------------------------------ */ inline void Parameter::setAuto(__attribute__((unused)) const ParserParameter & value) { if (!(isParsable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "parsable")); } /* -------------------------------------------------------------------------- */ template const T & Parameter::get() const { if (!(isReadable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "readable")); const ParameterTyped & typed_param = getParameterTyped(); return typed_param.getTyped(); } /* -------------------------------------------------------------------------- */ template T & Parameter::get() { ParameterTyped & typed_param = getParameterTyped(); if (!(isReadable()) || !(this->isWritable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "accessible")); return typed_param.getTyped(); } /* -------------------------------------------------------------------------- */ template inline Parameter::operator T() const { return this->get(); } /* -------------------------------------------------------------------------- */ template ParameterTyped::ParameterTyped(std::string name, std::string description, ParameterAccessType param_type, T & param) : Parameter(name, description, param_type), param(param) {} /* -------------------------------------------------------------------------- */ template template void ParameterTyped::setTyped(const V & value) { param = value; } /* -------------------------------------------------------------------------- */ template inline void ParameterTyped::setAuto(const ParserParameter & value) { Parameter::setAuto(value); param = static_cast(value); } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped::setAuto(const ParserParameter & value) { Parameter::setAuto(value); param = value.getValue(); } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped>::setAuto(const ParserParameter & in_param) { Parameter::setAuto(in_param); Vector tmp = in_param; - for (UInt i = 0; i < param.size(); ++i) { - param(i) = tmp(i); + if (param.size() == 0) { + param = tmp; + } else { + for (UInt i = 0; i < param.size(); ++i) { + param(i) = tmp(i); + } } } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped>::setAuto(const ParserParameter & in_param) { Parameter::setAuto(in_param); Matrix tmp = in_param; - for (UInt i = 0; i < param.rows(); ++i) { - for (UInt j = 0; j < param.cols(); ++j) { - param(i, j) = tmp(i, j); + if (param.size() == 0) { + param = tmp; + } else { + for (UInt i = 0; i < param.rows(); ++i) { + for (UInt j = 0; j < param.cols(); ++j) { + param(i, j) = tmp(i, j); + } } } } /* -------------------------------------------------------------------------- */ template const T & ParameterTyped::getTyped() const { return param; } /* -------------------------------------------------------------------------- */ template T & ParameterTyped::getTyped() { return param; } /* -------------------------------------------------------------------------- */ template inline void ParameterTyped::printself(std::ostream & stream) const { Parameter::printself(stream); stream << param << "\n"; } /* -------------------------------------------------------------------------- */ template class ParameterTyped> : public Parameter { public: ParameterTyped(std::string name, std::string description, ParameterAccessType param_type, std::vector & param) : Parameter(name, description, param_type), param(param) {} /* ------------------------------------------------------------------------ */ template void setTyped(const V & value) { param = value; } - void setAuto(const ParserParameter & value) override { Parameter::setAuto(value); + void setAuto(const ParserParameter & value) override { + Parameter::setAuto(value); param.clear(); - const std::vector & tmp = value; - for(auto && z : tmp) { + const std::vector & tmp = value; + for (auto && z : tmp) { param.emplace_back(z); } } std::vector & getTyped() { return param; } const std::vector & getTyped() const { return param; } void printself(std::ostream & stream) const override { Parameter::printself(stream); stream << "[ "; for (auto && v : param) stream << v << " "; stream << "]\n"; } inline const std::type_info & type() const override { return typeid(std::vector); } private: /// Value of parameter std::vector & param; }; /* ------o-------------------------------------------------------------------- */ template <> inline void ParameterTyped::printself(std::ostream & stream) const { Parameter::printself(stream); stream << std::boolalpha << param << "\n"; } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::registerParam(std::string name, T & variable, ParameterAccessType type, const std::string & description) { auto it = params.find(name); if (it != params.end()) AKANTU_CUSTOM_EXCEPTION(debug::ParameterException( name, "Parameter named " + name + " already registered.")); auto * param = new ParameterTyped(name, description, type, variable); params[name] = param; } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::registerParam(std::string name, T & variable, const T & default_value, ParameterAccessType type, const std::string & description) { variable = default_value; registerParam(name, variable, type, description); } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::setMixed(const std::string & name, const V & value) { auto it = params.find(name); if (it == params.end()) { if (consisder_sub) { for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) { it->second->setMixed(name, value); } } else { AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this)); } } else { Parameter & param = *(it->second); param.set(value); } } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::set(const std::string & name, const T & value) { this->template setMixed(name, value); } /* -------------------------------------------------------------------------- */ template T & ParameterRegistry::get(const std::string & name) { auto it = params.find(name); if (it == params.end()) { if (consisder_sub) { for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) { try { return it->second->get(name); } catch (...) { } } } // nothing was found not even in sub registries AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this)); } Parameter & param = *(it->second); return param.get(); } /* -------------------------------------------------------------------------- */ const Parameter & ParameterRegistry::get(const std::string & name) const { auto it = params.find(name); if (it == params.end()) { if (consisder_sub) { for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) { try { return it->second->get(name); } catch (...) { } } } // nothing was found not even in sub registries AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this)); } Parameter & param = *(it->second); return param; } /* -------------------------------------------------------------------------- */ namespace { namespace details { template struct CastHelper { static R convert(const T &) { throw std::bad_cast(); } }; template struct CastHelper::value>> { static R convert(const T & val) { return val; } }; } // namespace details } // namespace template inline ParameterTyped::operator Real() const { if (not isReadable()) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "accessible")); return details::CastHelper::convert(param); } } // namespace akantu #endif /* __AKANTU_PARAMETER_REGISTRY_TMPL_HH__ */ diff --git a/src/io/parser/parser_tmpl.hh b/src/io/parser/parser_tmpl.hh index c13fc4500..c65c8d82c 100644 --- a/src/io/parser/parser_tmpl.hh +++ b/src/io/parser/parser_tmpl.hh @@ -1,118 +1,111 @@ /** * @file parser_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Apr 22 2015 * * @brief Implementation of the parser templated methods * * @section LICENSE * * Copyright (©) 2014, 2015 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 /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template -inline ParserParameter::operator T() const { +template inline ParserParameter::operator T() const { T t; std::stringstream sstr(value); sstr >> t; - if(sstr.bad()) + if (sstr.bad()) AKANTU_EXCEPTION("No known conversion of a ParserParameter \"" - << name << "\" to the type " - << typeid(T).name()); + << name << "\" to the type " << typeid(T).name()); return t; } /* -------------------------------------------------------------------------- */ -template<> -inline ParserParameter::operator const char *() const { +template <> inline ParserParameter::operator const char *() const { return value.c_str(); } /* -------------------------------------------------------------------------- */ -template<> -inline ParserParameter::operator Real() const { +template <> inline ParserParameter::operator Real() const { return Parser::parseReal(value, *parent_section); } -/* --------------------------------------------------------- ----------------- */ -template<> -inline ParserParameter::operator bool() const { +/* --------------------------------------------------------- ----------------- + */ +template <> inline ParserParameter::operator bool() const { bool b; std::stringstream sstr(value); sstr >> std::boolalpha >> b; - if(sstr.fail()) { + if (sstr.fail()) { sstr.clear(); sstr >> std::noboolalpha >> b; } return b; } /* -------------------------------------------------------------------------- */ -template<> -inline ParserParameter::operator std::vector() const { +template <> inline ParserParameter::operator std::vector() const { std::vector tmp; - std::regex pieces_regex("\\[([^ ,]+)(\\s*,\\s*([^ ,]+))*\\]"); - std::smatch pieces_match; - if (std::regex_match(value, pieces_match, pieces_regex)) { - for(size_t i = 1; i < pieces_match.size(); i+=2) { - tmp.push_back(pieces_match[i].str()); - } + auto string = + std::regex_replace(value, std::regex("[[:space:]]|\\[|\\]"), ""); + std::smatch sm; + while (std::regex_search(string, sm, std::regex("[^,]+"))) { + tmp.push_back(sm.str()); + string = sm.suffix(); } return tmp; } - -/* --------------------------------------------------------- ----------------- */ -template<> -inline ParserParameter::operator Vector() const { +/* --------------------------------------------------------- ----------------- + */ +template <> inline ParserParameter::operator Vector() const { return Parser::parseVector(value, *parent_section); } -/* --------------------------------------------------------- ----------------- */ -template<> -inline ParserParameter::operator Vector() const { +/* --------------------------------------------------------- ----------------- + */ +template <> inline ParserParameter::operator Vector() const { Vector tmp = Parser::parseVector(value, *parent_section); Vector tmp_uint(tmp.size()); - for (UInt i=0; i -inline ParserParameter::operator Matrix() const { +/* --------------------------------------------------------- ----------------- + */ +template <> inline ParserParameter::operator Matrix() const { return Parser::parseMatrix(value, *parent_section); } /* -------------------------------------------------------------------------- */ -template<> -inline ParserParameter::operator RandomParameter() const { +template <> inline ParserParameter::operator RandomParameter() const { return Parser::parseRandomParameter(value, *parent_section); } } // akantu diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index 4a99d46a9..0e4b72f0c 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,433 +1,436 @@ /** * @file element_type_map.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief storage class by element type * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_memory.hh" #include "aka_named_argument.hh" #include "element.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_HH__ namespace akantu { class FEEngine; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(all_ghost_types); DECLARE_NAMED_ARGUMENT(default_value); DECLARE_NAMED_ARGUMENT(element_kind); DECLARE_NAMED_ARGUMENT(ghost_type); DECLARE_NAMED_ARGUMENT(nb_component); DECLARE_NAMED_ARGUMENT(with_nb_element); DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element); DECLARE_NAMED_ARGUMENT(spatial_dimension); DECLARE_NAMED_ARGUMENT(do_not_default); } // namespace template class ElementTypeMap; /* -------------------------------------------------------------------------- */ /* ElementTypeMapBase */ /* -------------------------------------------------------------------------- */ /// Common non templated base class for the ElementTypeMap class class ElementTypeMapBase { public: virtual ~ElementTypeMapBase() = default; }; /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template class ElementTypeMap : public ElementTypeMapBase { public: ElementTypeMap(); ~ElementTypeMap() override; inline static std::string printType(const SupportType & type, const GhostType & ghost_type); /*! Tests whether a type is present in the object * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return true if the type is present. */ inline bool exists(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline const Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. THIS METHOD IS * NOT ARRAY SAFE, when using ElementTypeMapArray, use setArray instead * @param data to insert * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const Stored & insert, const SupportType & type, const GhostType & ghost_type = _not_ghost); /// print helper virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ /*! iterator allows to iterate over type-data pairs of the map. The interface * expects the SupportType to be ElementType. */ using DataMap = std::map; class type_iterator : private std::iterator { public: using value_type = const SupportType; using pointer = const SupportType *; using reference = const SupportType &; protected: using DataMapIterator = typename ElementTypeMap::DataMap::const_iterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); type_iterator() = default; inline reference operator*(); inline reference operator*() const; inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other) const; inline bool operator!=(const type_iterator & other) const; type_iterator & operator=(const type_iterator & other); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; /// helper class to use in range for constructions class ElementTypesIteratorHelper { public: using Container = ElementTypeMap; using iterator = typename Container::type_iterator; ElementTypesIteratorHelper(const Container & container, UInt dim, GhostType ghost_type, ElementKind kind) : container(std::cref(container)), dim(dim), ghost_type(ghost_type), kind(kind) {} template ElementTypesIteratorHelper(const Container & container, use_named_args_t, pack &&... _pack) : ElementTypesIteratorHelper( container, OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), OPTIONAL_NAMED_ARG(ghost_type, _not_ghost), OPTIONAL_NAMED_ARG(element_kind, _ek_regular)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) = default; iterator begin() { return container.get().firstType(dim, ghost_type, kind); } iterator end() { return container.get().lastType(dim, ghost_type, kind); } private: std::reference_wrapper container; UInt dim; GhostType ghost_type; ElementKind kind; }; private: ElementTypesIteratorHelper elementTypesImpl(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(use_named_args, std::forward(_pack)...); } template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(std::forward(_pack)...); } public: /*! Get an iterator to the beginning of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the first stored data matching the filters * or an iterator to the end of the map if none match*/ inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! Get an iterator to the end of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the last stored data matching the filters * or an iterator to the end of the map if none match */ inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; protected: /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline DataMap & getData(GhostType ghost_type); /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline const DataMap & getData(GhostType ghost_type) const; /* ------------------------------------------------------------------------ */ protected: DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ElementTypeMapArray : public ElementTypeMap *, SupportType>, public Memory { public: using type = T; using array_type = Array; protected: using parent = ElementTypeMap *, SupportType>; using DataMap = typename parent::DataMap; public: using type_iterator = typename parent::type_iterator; /// standard assigment (copy) operator void operator=(const ElementTypeMapArray &) = delete; ElementTypeMapArray(const ElementTypeMapArray &) = delete; /*! Constructor * @param id optional: identifier (string) * @param parent_id optional: parent identifier. for organizational purposes * only * @param memory_id optional: choose a specific memory, defaults to memory 0 */ ElementTypeMapArray(const ID & id = "by_element_type_array", const ID & parent_id = "no_parent", const MemoryID & memory_id = 0) : parent(), Memory(parent_id + ":" + id, memory_id), name(id){}; /*! allocate memory for a new array * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param ghost_type whether to add the field to the data map or the * ghost_data map * @return a reference to the allocated array */ inline Array & alloc(UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value = T()); /*! allocate memory for a new array in both the data and the ghost_data map * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map*/ inline void alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value = T()); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a reference to the array */ inline const Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /// access the data of an element, this combine the map and array accessor inline const T & operator()(const Element & element, UInt component = 0) const; /// access the data of an element, this combine the map and array accessor inline T & operator()(const Element & element, UInt component = 0); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a const reference to the array */ inline Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @param vect the vector to include into the map * @return stored data corresponding to type. */ inline void setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect); /*! frees all memory related to the data*/ inline void free(); /*! set all values in the ElementTypeMap to zero*/ inline void clear(); + /*! set all values in the ElementTypeMap to value */ + inline void set(const T & value); + /*! deletes and reorders entries in the stored arrays * @param new_numbering a ElementTypeMapArray of new indices. UInt(-1) * indicates * deleted entries. */ inline void onElementsRemoved(const ElementTypeMapArray & new_numbering); /// text output helper void printself(std::ostream & stream, int indent = 0) const override; /*! set the id * @param id the new name */ inline void setID(const ID & id) { this->id = id; } ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const { ElementTypeMap nb_components; for (auto & type : this->elementTypes(dim, ghost_type, kind)) { UInt nb_comp = (*this)(type, ghost_type).getNbComponent(); nb_components(type, ghost_type) = nb_comp; } return nb_components; } /* ------------------------------------------------------------------------ */ /* more evolved allocators */ /* ------------------------------------------------------------------------ */ public: /// initialize the arrays in accordance to a functor template void initialize(const Func & f, const T & default_value, bool do_not_default); /// initialize with sizes and number of components in accordance of a mesh /// content template void initialize(const Mesh & mesh, pack &&... _pack); /// initialize with sizes and number of components in accordance of a fe /// engine content (aka integration points) template void initialize(const FEEngine & fe_engine, pack &&... _pack); /* ------------------------------------------------------------------------ */ /* Accesssors */ /* ------------------------------------------------------------------------ */ public: /// get the name of the internal field AKANTU_GET_MACRO(Name, name, ID); /// name of the elment type map: e.g. connectivity, grad_u ID name; }; /// to store data Array by element type using ElementTypeMapReal = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapInt = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapUInt = ElementTypeMapArray; /// Map of data of type UInt stored in a mesh using UIntDataMap = std::map *>; using ElementTypeMapUIntDataMap = ElementTypeMap; } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_HH__ */ diff --git a/src/mesh/element_type_map_tmpl.hh b/src/mesh/element_type_map_tmpl.hh index 44d39fe21..c74978ee8 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,648 +1,660 @@ /** * @file element_type_map_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief implementation of template functions of the ElementTypeMap and * ElementTypeMapArray classes * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_static_if.hh" #include "element_type_map.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "element_type_conversion.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template inline std::string ElementTypeMap::printType(const SupportType & type, const GhostType & ghost_type) { std::stringstream sstr; sstr << "(" << ghost_type << ":" << type << ")"; return sstr.str(); } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::exists( const SupportType & type, const GhostType & ghost_type) const { return this->getData(ghost_type).find(type) != this->getData(ghost_type).end(); } /* -------------------------------------------------------------------------- */ template inline const Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMap::printType(type, ghost_type) << " in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); return it->second; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) { return this->getData(ghost_type)[type]; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap:: operator()(const Stored & insert, const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it != this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("Element of type " << ElementTypeMap::printType(type, ghost_type) << " already in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); } else { auto & data = this->getData(ghost_type); const auto & res = data.insert(std::pair(type, insert)); it = res.first; } return it->second; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ template inline const typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) const { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ /// Works only if stored is a pointer to a class with a printself method template void ElementTypeMap::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> [" << std::endl; for (auto gt : ghost_types) { const DataMap & data = getData(gt); for (auto & pair : data) { stream << space << space << ElementTypeMap::printType(pair.first, gt) << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template ElementTypeMap::ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ template ElementTypeMap::~ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ /* ElementTypeMapArray */ /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray::alloc( UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value) { std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; Array * tmp; auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { std::stringstream sstr; sstr << this->id << ":" << type << ghost_id; tmp = &(Memory::alloc(sstr.str(), size, nb_component, default_value)); std::stringstream sstrg; sstrg << ghost_type; // tmp->setTag(sstrg.str()); this->getData(ghost_type)[type] = tmp; } else { AKANTU_DEBUG_INFO( "The vector " << this->id << this->printType(type, ghost_type) << " already exists, it is resized instead of allocated."); tmp = it->second; it->second->resize(size); } return *tmp; } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value) { this->alloc(size, nb_component, type, _not_ghost, default_value); this->alloc(size, nb_component, type, _ghost, default_value); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::free() { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & pair : data) { dealloc(pair.second->getID()); } data.clear(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::clear() { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->clear(); } } } +/* -------------------------------------------------------------------------- */ +template +inline void ElementTypeMapArray::set(const T & value) { + for (auto gt : ghost_types) { + auto & data = this->getData(gt); + for (auto & vect : data) { + vect.second->set(value); + } + } +} + /* -------------------------------------------------------------------------- */ template inline const Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this const ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class(\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class (\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect) { auto it = this->getData(ghost_type).find(type); if (AKANTU_DEBUG_TEST(dblWarning) && it != this->getData(ghost_type).end() && it->second != &vect) { AKANTU_DEBUG_WARNING( "The Array " << this->printType(type, ghost_type) << " is already registred, this call can lead to a memory leak."); } this->getData(ghost_type)[type] = &(const_cast &>(vect)); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::onElementsRemoved( const ElementTypeMapArray & new_numbering) { for (auto gt : ghost_types) { for (auto & type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { auto support_type = convertType(type); if (this->exists(support_type, gt)) { const auto & renumbering = new_numbering(type, gt); if (renumbering.size() == 0) continue; auto & vect = this->operator()(support_type, gt); auto nb_component = vect.getNbComponent(); Array tmp(renumbering.size(), nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.size(); ++i) { UInt new_i = renumbering(i); if (new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i * nb_component, nb_component * sizeof(T)); ++new_size; } } tmp.resize(new_size); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; for (UInt g = _not_ghost; g <= _ghost; ++g) { auto gt = (GhostType)g; const DataMap & data = this->getData(gt); typename DataMap::const_iterator it; for (it = data.begin(); it != data.end(); ++it) { stream << space << space << ElementTypeMapArray::printType(it->first, gt) << " [" << std::endl; it->second->printself(stream, indent + 3); stream << space << space << " ]" << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /* SupportType Iterator */ /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek) : list_begin(list_begin), list_end(list_end), dim(dim), kind(ek) {} /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( const type_iterator & it) : list_begin(it.list_begin), list_end(it.list_end), dim(it.dim), kind(it.kind) {} /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator:: operator=(const type_iterator & it) { if (this != &it) { list_begin = it.list_begin; list_end = it.list_end; dim = it.dim; kind = it.kind; } return *this; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() const { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator::operator++() { ++list_begin; while ((list_begin != list_end) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(list_begin->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(list_begin->first))))) ++list_begin; return *this; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator ElementTypeMap::type_iterator::operator++(int) { type_iterator tmp(*this); operator++(); return tmp; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator==(const type_iterator & other) const { return this->list_begin == other.list_begin; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator!=(const type_iterator & other) const { return this->list_begin != other.list_begin; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl(UInt dim, GhostType ghost_type, ElementKind kind) const { return ElementTypesIteratorHelper(*this, dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ template template typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl( const use_named_args_t & unused, pack &&... _pack) const { return ElementTypesIteratorHelper(*this, unused, _pack...); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::firstType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator b, e; b = getData(ghost_type).begin(); e = getData(ghost_type).end(); // loop until the first valid type while ((b != e) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(b->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(b->first))))) ++b; return typename ElementTypeMap::type_iterator(b, e, dim, kind); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::lastType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator e; e = getData(ghost_type).end(); return typename ElementTypeMap::type_iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const ElementTypeMap & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ class ElementTypeMapArrayInializer { public: ElementTypeMapArrayInializer(UInt spatial_dimension = _all_dimensions, UInt nb_component = 1, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular) : spatial_dimension(spatial_dimension), nb_component(nb_component), ghost_type(ghost_type), element_kind(element_kind) {} const GhostType & ghostType() const { return ghost_type; } protected: UInt spatial_dimension; UInt nb_component; GhostType ghost_type; ElementKind element_kind; }; /* -------------------------------------------------------------------------- */ class MeshElementTypeMapArrayInializer : public ElementTypeMapArrayInializer { public: MeshElementTypeMapArrayInializer( const Mesh & mesh, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular, bool with_nb_element = false, bool with_nb_nodes_per_element = false) : ElementTypeMapArrayInializer(spatial_dimension, nb_component, ghost_type, element_kind), mesh(mesh), with_nb_element(with_nb_element), with_nb_nodes_per_element(with_nb_nodes_per_element) {} decltype(auto) elementTypes() const { return mesh.elementTypes(this->spatial_dimension, this->ghost_type, this->element_kind); } virtual UInt size(const ElementType & type) const { if (with_nb_element) return mesh.getNbElement(type, this->ghost_type); return 0; } UInt nbComponent(const ElementType & type) const { if (with_nb_nodes_per_element) return (this->nb_component * mesh.getNbNodesPerElement(type)); return this->nb_component; } protected: const Mesh & mesh; bool with_nb_element; bool with_nb_nodes_per_element; }; /* -------------------------------------------------------------------------- */ class FEEngineElementTypeMapArrayInializer : public MeshElementTypeMapArrayInializer { public: FEEngineElementTypeMapArrayInializer( const FEEngine & fe_engine, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular); UInt size(const ElementType & type) const override; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; ElementTypesIteratorHelper elementTypes() const; protected: const FEEngine & fe_engine; }; /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const Func & f, const T & default_value, bool do_not_default) { for (auto & type : f.elementTypes()) { if (not this->exists(type, f.ghostType())) if (do_not_default) { auto & array = this->alloc(0, f.nbComponent(type), type, f.ghostType()); array.resize(f.size(type)); } else { this->alloc(f.size(type), f.nbComponent(type), type, f.ghostType(), default_value); } else { auto & array = this->operator()(type, f.ghostType()); - array.resize(f.size(type)); if (not do_not_default) - array.set(default_value); + array.resize(f.size(type), default_value); + else + array.resize(f.size(type)); } } } /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const Mesh & mesh, pack &&... _pack) { GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); bool all_ghost_types = requested_ghost_type == _casper; for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; this->initialize( MeshElementTypeMapArrayInializer( mesh, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, mesh.getSpatialDimension()), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular), OPTIONAL_NAMED_ARG(with_nb_element, false), OPTIONAL_NAMED_ARG(with_nb_nodes_per_element, false)), OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const FEEngine & fe_engine, pack &&... _pack) { bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, true); GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _not_ghost); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; this->initialize(FEEngineElementTypeMapArrayInializer( fe_engine, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, UInt(-2)), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular)), OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template inline T & ElementTypeMapArray:: operator()(const Element & element, UInt component) { return this->operator()(element.type, element.ghost_type)(element.element, component); } /* -------------------------------------------------------------------------- */ template inline const T & ElementTypeMapArray:: operator()(const Element & element, UInt component) const { return this->operator()(element.type, element.ghost_type)(element.element, component); } } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ */ diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index e72a74112..877443b78 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,476 +1,535 @@ /** * @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: Fri Jan 22 2016 * * @brief class handling meshes * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_io.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ +#include "communicator.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "mesh_utils_distribution.hh" #include "node_synchronizer.hh" -#include "communicator.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumper_field.hh" #include "dumper_internal_material_field.hh" #endif /* -------------------------------------------------------------------------- */ #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), nodes_type(0, 1, id + ":nodes_type"), connectivities("connectivities", id, memory_id), ghosts_counters("ghosts_counters", id, memory_id), normals("normals", id, memory_id), spatial_dimension(spatial_dimension), lower_bounds(spatial_dimension, 0.), upper_bounds(spatial_dimension, 0.), size(spatial_dimension, 0.), local_lower_bounds(spatial_dimension, 0.), local_upper_bounds(spatial_dimension, 0.), mesh_data("mesh_data", id, memory_id), communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, - const ID & id, const MemoryID & memory_id) +Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id, + const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id, memory_id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ 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; 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); + + for_each_elements( + mesh, + [&](auto && element) { + Vector barycenter(spatial_dimension); + mesh.getBarycenter(element, barycenter); + auto norm_barycenter = barycenter.norm(); + + const auto & element_to_facet = mesh_facets->getElementToSubelement( + element.type, element.ghost_type); + + Vector barycenter_facet(spatial_dimension); + auto nb_facet = + mesh_facets->getNbElement(element.type, element.ghost_type); + + auto range = arange(nb_facet); + + // this is a spacial search coded the most inefficient way. + auto facet = + std::find_if(range.begin(), range.end(), [&](auto && facet) { + if (element_to_facet(facet)[1] == ElementNull) + return false; + + mesh_facets->getBarycenter( + Element{element.type, facet, element.ghost_type}, + barycenter_facet); + auto norm_distance = barycenter_facet.distance(barycenter); + + return (norm_distance < + (Math::getTolerance() * norm_barycenter)); + }); + if (facet == range.end()) + AKANTU_EXCEPTION("The element " + << element + << " did not find its associated facet in the " + "mesh_facets! Try to decrease math tolerance."); + + // set physical name + phys_data(Element{element.type, *facet, element.ghost_type}) = + mesh_phys_data(element); + }, + _spatial_dimension = spatial_dimension - 1); + + mesh_facets->createGroupsFromMeshData("physical_names"); } AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO mesh_io; mesh_io.read(filename, *this, mesh_io_type); type_iterator it = this->firstType(spatial_dimension, _not_ghost, _ek_not_defined); type_iterator last = this->lastType(spatial_dimension, _not_ghost, _ek_not_defined); if (it == last) AKANTU_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_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::write(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO mesh_io; mesh_io.write(filename, *this, mesh_io_type); } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Mesh [" << std::endl; stream << space << " + id : " << getID() << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent + 2); stream << space << " + connectivities [" << std::endl; connectivities.printself(stream, indent + 2); stream << space << " ]" << std::endl; GroupManager::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox() { AKANTU_DEBUG_IN(); for (UInt k = 0; k < spatial_dimension; ++k) { local_lower_bounds(k) = std::numeric_limits::max(); local_upper_bounds(k) = -std::numeric_limits::max(); } for (UInt i = 0; i < nodes->size(); ++i) { // if(!isPureGhostNode(i)) for (UInt k = 0; k < spatial_dimension; ++k) { local_lower_bounds(k) = std::min(local_lower_bounds[k], (*nodes)(i, k)); local_upper_bounds(k) = std::max(local_upper_bounds[k], (*nodes)(i, k)); } } if (this->is_distributed) { Matrix reduce_bounds(spatial_dimension, 2); for (UInt k = 0; k < spatial_dimension; ++k) { reduce_bounds(k, 0) = local_lower_bounds(k); reduce_bounds(k, 1) = -local_upper_bounds(k); } communicator->allReduce(reduce_bounds, SynchronizerOperation::_min); for (UInt k = 0; k < spatial_dimension; ++k) { lower_bounds(k) = reduce_bounds(k, 0); upper_bounds(k) = -reduce_bounds(k, 1); } } else { this->lower_bounds = this->local_lower_bounds; this->upper_bounds = this->local_upper_bounds; } size = upper_bounds - lower_bounds; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : global_connectivity.elementTypes(_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(); if (not nodes_global_ids && is_mesh_facets) { std::transform( local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [& node_ids = *mesh_parent->nodes_global_ids](UInt l)->UInt { return node_ids(l); }); } else { std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [& node_ids = *nodes_global_ids](UInt l)->UInt { return node_ids(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(); #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_DEBUG_ERROR("Cannot distribute a mesh without a partitioning tool"); } #endif 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); } } } } /* -------------------------------------------------------------------------- */ 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 80da370c0..fbe7d74ee 100644 --- a/src/mesh/mesh.hh +++ b/src/mesh/mesh.hh @@ -1,601 +1,604 @@ /** * @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: Thu Jan 14 2016 * * @brief the class representing the meshes * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_event_handler_manager.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element.hh" #include "element_class.hh" #include "element_type_map.hh" #include "group_manager.hh" #include "mesh_data.hh" #include "mesh_events.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class Communicator; class ElementSynchronizer; class NodeSynchronizer; } // 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} Mesh::type_iterator it = mesh.firstType(dim, ghost_type); Mesh::type_iterator end = mesh.lastType(dim, ghost_type); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); const Array & conn = mesh.getConnectivity(*it); for(UInt e = 0; e < nb_element; ++e) { ... } } @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 knowing its /// ID // Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id, // const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); ~Mesh() override; /// @typedef ConnectivityTypeList list of the types present in a Mesh using ConnectivityTypeList = std::set; /// read the mesh from a file void read(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); /// write the mesh to a file void write(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); private: /// initialize the connectivity to NULL and other stuff void init(); /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// patitionate the mesh among the processors involved in their computation virtual void distribute(Communicator & communicator); virtual void distribute(); /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Array & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; // /// extract coordinates of nodes from a reversed element // inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, // UInt * connectivity, // UInt n_nodes); /// convert a element to a linearized element inline UInt elementToLinearized(const Element & elem) const; /// convert a linearized element to an element inline Element linearizedToElement(UInt linearized_element) const; /// update the types offsets array for the conversions // inline void updateTypesOffsets(const GhostType & ghost_type); /// 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(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the id of the mesh AKANTU_GET_MACRO(ID, Memory::id, const ID &); /// get the id of the mesh AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &); /// get the spatial dimension of the mesh = number of component of the /// coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Array aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Array &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array &); /// get the normals for the elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt); /// get the Array of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array &); AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global id of a node inline UInt getNodeLocalId(UInt global_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Array AKANTU_GET_MACRO(NodesType, nodes_type, const Array &); protected: AKANTU_GET_MACRO_NOT_CONST(NodesType, nodes_type, Array &); public: inline NodeType getNodeType(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; AKANTU_GET_MACRO(LowerBounds, lower_bounds, const Vector &); AKANTU_GET_MACRO(UpperBounds, upper_bounds, const Vector &); AKANTU_GET_MACRO(LocalLowerBounds, local_lower_bounds, const Vector &); AKANTU_GET_MACRO(LocalUpperBounds, local_upper_bounds, const Vector &); /// get the connectivity Array for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); AKANTU_GET_MACRO(Connectivities, connectivities, const ElementTypeMapArray &); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get the number of element for a given ghost_type and a given dimension inline UInt getNbElement(const UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_not_defined) const; // /// get the connectivity list either for the elements or the ghost elements // inline const ConnectivityTypeList & // getConnectivityTypeList(const GhostType & ghost_type = _not_ghost) const; /// compute the barycenter of a given element private: inline void getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type = _not_ghost) const; public: inline void getBarycenter(const Element & element, Vector & barycenter) const; /// get the element connected to a subelement 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 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); /// 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); /// Facets mesh accessor inline const Mesh & getMeshFacets() const; inline Mesh & getMeshFacets(); /// Parent mesh accessor inline const Mesh & getMeshParent() const; inline bool isMeshFacets() const { return this->is_mesh_facets; } /// defines is the mesh is distributed or not inline bool isDistributed() const { return this->is_distributed; } #ifndef SWIG /// return the dumper from a group and and a dumper name DumperIOHelper & getGroupDumper(const std::string & dumper_name, const std::string & group_name); #endif /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type, UInt t); /// 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; /// 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); /// 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); /// get the number of nodes in the given element list static inline UInt getNbNodesPerElementList(const Array & elements); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ using type_iterator = ElementTypeMapArray::type_iterator; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; template ElementTypesIteratorHelper elementTypes(pack &&... _pack) const; inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.firstType(dim, ghost_type, kind); } inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.lastType(dim, ghost_type, kind); } AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer, const ElementSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer, ElementSynchronizer &); AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer, const NodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer, NodeSynchronizer &); // AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator // &); #ifndef SWIG AKANTU_GET_MACRO(Communicator, *communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &); #endif /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshAccessor; AKANTU_GET_MACRO(NodesPointer, *nodes, Array &); /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline Array & getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Array and create it if necessary inline Array & getNodesTypePointer(); /// get a pointer to the connectivity Array for the given type and create it /// if necessary inline Array & getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the ghost element counter inline Array & getGhostsCounters(const ElementType & type, const GhostType & ghost_type = _ghost) { AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost, "No ghost counter for _not_ghost elements"); return ghosts_counters(type, ghost_type); } /// get a pointer to the element_to_subelement Array for the given type and /// create it if necessary inline Array> & getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Array for the given type and /// create it if necessary inline Array & getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); AKANTU_GET_MACRO_NOT_CONST(MeshData, mesh_data, MeshData &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of the nodes coordinates std::shared_ptr> nodes; /// global node ids std::unique_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 Array nodes_type; /// global number of nodes; UInt nb_global_nodes{0}; /// all class of elements present in this mesh (for heterogenous meshes) ElementTypeMapArray connectivities; /// count the references on ghost elements ElementTypeMapArray ghosts_counters; /// map to normals for all class of elements present in this mesh ElementTypeMapArray normals; /// list of all existing types in the mesh // ConnectivityTypeList type_set; /// the spatial dimension of this mesh UInt spatial_dimension{0}; /// types offsets // Array types_offsets; // /// list of all existing types in the mesh // ConnectivityTypeList ghost_type_set; // /// ghost types offsets // Array ghost_types_offsets; /// min of coordinates Vector lower_bounds; /// max of coordinates Vector upper_bounds; /// size covered by the mesh on each direction Vector size; /// local min of coordinates Vector local_lower_bounds; /// local max of coordinates Vector local_upper_bounds; /// Extra data loaded from the mesh file MeshData mesh_data; /// facets' mesh std::unique_ptr mesh_facets; /// parent mesh (this is set for mesh_facets meshes) const Mesh * mesh_parent{nullptr}; /// defines if current mesh is mesh_facets or not bool is_mesh_facets{false}; /// defines if the mesh is centralized or distributed bool is_distributed{false}; /// Communicator on which mesh is distributed Communicator * communicator; /// Element synchronizer std::unique_ptr element_synchronizer; /// Node synchronizer std::unique_ptr node_synchronizer; using NodesToElements = std::vector>>; /// This info is stored to simplify the dynamic changes NodesToElements nodes_to_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ #include "element_type_map_tmpl.hh" #include "mesh_inline_impl.cc" #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/mesh/mesh_data.hh b/src/mesh/mesh_data.hh index 144c399c5..029af4fe9 100644 --- a/src/mesh/mesh_data.hh +++ b/src/mesh/mesh_data.hh @@ -1,150 +1,153 @@ /** * @file mesh_data.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Thu Nov 05 2015 * * @brief Stores generic data loaded from the mesh file * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_DATA_HH__ #define __AKANTU_MESH_DATA_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "element_type_map.hh" #include #include /* -------------------------------------------------------------------------- */ namespace akantu { #define AKANTU_MESH_DATA_TYPES \ ((_tc_int, Int))((_tc_uint, UInt))((_tc_real, Real))( \ (_tc_element, Element))((_tc_std_string, std::string))( \ (_tc_std_vector_element, std::vector)) #define AKANTU_MESH_DATA_TUPLE_FIRST_ELEM(s, data, elem) \ BOOST_PP_TUPLE_ELEM(2, 0, elem) enum MeshDataTypeCode : int { BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_MESH_DATA_TUPLE_FIRST_ELEM, , AKANTU_MESH_DATA_TYPES)), _tc_unknown }; class MeshData : public Memory { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: using TypeCode = MeshDataTypeCode; using ElementalDataMap = std::map; using StringVector = std::vector; using TypeCodeMap = std::map; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshData(const ID & id = "mesh_data", const ID & parent_id = "", const MemoryID & memory_id = 0); ~MeshData() override; /* ------------------------------------------------------------------------ */ /* Methods and accessors */ /* ------------------------------------------------------------------------ */ public: /// Register new elemental data (and alloc data) with check if the name is /// new template void registerElementalData(const ID & name); inline void registerElementalData(const ID & name, TypeCode type); - /// Get an existing elemental data array + /// tells if the given array exists template bool hasDataArray(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; + /// tells if the given data exists + bool hasData(const ID & data_name) const; + /// Get an existing elemental data array template const Array & getElementalDataArray(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; template Array & getElementalDataArray(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// Get an elemental data array, if it does not exist: allocate it template Array & getElementalDataArrayAlloc(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost, UInt nb_component = 1); /// get the names of the data stored in elemental_data inline void getTagNames(StringVector & tags, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get the type of the data stored in elemental_data template TypeCode getTypeCode() const; inline TypeCode getTypeCode(const ID & name) const; template inline UInt getNbComponentTemplated(const ID & name, const ElementType & el_type, const GhostType & ghost_type) const; inline UInt getNbComponent(const ID & name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// Get an existing elemental data template const ElementTypeMapArray & getElementalData(const ID & name) const; template ElementTypeMapArray & getElementalData(const ID & name); private: /// Register new elemental data (add alloc data) template ElementTypeMapArray * allocElementalData(const ID & name); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Map when elemental data is stored as ElementTypeMap ElementalDataMap elemental_data; /// Map when elementalType of the data stored in elemental_data TypeCodeMap typecode_map; }; } // akantu #include "mesh_data_tmpl.hh" #undef AKANTU_MESH_DATA_TUPLE_FIRST_ELEM #endif /* __AKANTU_MESH_DATA_HH__ */ diff --git a/src/mesh/mesh_data_tmpl.hh b/src/mesh/mesh_data_tmpl.hh index d915f1b07..2ec591f4f 100644 --- a/src/mesh/mesh_data_tmpl.hh +++ b/src/mesh/mesh_data_tmpl.hh @@ -1,272 +1,276 @@ /** * @file mesh_data_tmpl.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Thu Nov 05 2015 * * @brief Stores generic data loaded from the mesh file * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_data.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_DATA_TMPL_HH__ #define __AKANTU_MESH_DATA_TMPL_HH__ namespace akantu { #define MESH_DATA_GET_TYPE(r, data, type) \ template <> \ inline MeshDataTypeCode \ MeshData::getTypeCode() const { \ return BOOST_PP_TUPLE_ELEM(2, 0, type); \ } /* -------------------------------------------------------------------------- */ // get the type of the data stored in elemental_data template inline MeshDataTypeCode MeshData::getTypeCode() const { AKANTU_DEBUG_ERROR("Type " << debug::demangle(typeid(T).name()) << "not implemented by MeshData."); } /* -------------------------------------------------------------------------- */ BOOST_PP_SEQ_FOR_EACH(MESH_DATA_GET_TYPE, void, AKANTU_MESH_DATA_TYPES) #undef MESH_DATA_GET_TYPE inline MeshDataTypeCode MeshData::getTypeCode(const ID & name) const { auto it = typecode_map.find(name); if (it == typecode_map.end()) AKANTU_EXCEPTION("No dataset named " << name << " found."); return it->second; } /* -------------------------------------------------------------------------- */ // Register new elemental data templated (and alloc data) with check if the // name is new template void MeshData::registerElementalData(const ID & name) { auto it = elemental_data.find(name); if (it == elemental_data.end()) { allocElementalData(name); } else { AKANTU_DEBUG_INFO("Data named " << name << " already registered."); } } /* -------------------------------------------------------------------------- */ // Register new elemental data of a given MeshDataTypeCode with check if the // name is new #define AKANTU_MESH_DATA_CASE_MACRO(r, name, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ registerElementalData(name); \ break; \ } inline void MeshData::registerElementalData(const ID & name, MeshDataTypeCode type) { switch (type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_MESH_DATA_CASE_MACRO, name, AKANTU_MESH_DATA_TYPES) default: AKANTU_DEBUG_ERROR("Type " << type << "not implemented by MeshData."); } } #undef AKANTU_MESH_DATA_CASE_MACRO /* -------------------------------------------------------------------------- */ /// Register new elemental data (and alloc data) template ElementTypeMapArray * MeshData::allocElementalData(const ID & name) { auto * dataset = new ElementTypeMapArray(name, id, memory_id); elemental_data[name] = dataset; typecode_map[name] = getTypeCode(); return dataset; } /* -------------------------------------------------------------------------- */ template const ElementTypeMapArray & MeshData::getElementalData(const ID & name) const { auto it = elemental_data.find(name); if (it == elemental_data.end()) AKANTU_EXCEPTION("No dataset named " << name << " found."); return dynamic_cast &>(*(it->second)); } /* -------------------------------------------------------------------------- */ // Get an existing elemental data template ElementTypeMapArray & MeshData::getElementalData(const ID & name) { auto it = elemental_data.find(name); if (it == elemental_data.end()) AKANTU_EXCEPTION("No dataset named " << name << " found."); return dynamic_cast &>(*(it->second)); } /* -------------------------------------------------------------------------- */ -// Get an existing elemental data array for an element type template bool MeshData::hasDataArray(const ID & name, const ElementType & elem_type, const GhostType & ghost_type) const { auto it = elemental_data.find(name); if (it == elemental_data.end()) return false; auto & elem_map = dynamic_cast &>(*(it->second)); return elem_map.exists(elem_type, ghost_type); } /* -------------------------------------------------------------------------- */ -// Get an existing elemental data array for an element type +inline bool MeshData::hasData(const ID & name) const { + auto it = elemental_data.find(name); + return (it != elemental_data.end()); +} + +/* -------------------------------------------------------------------------- */ template const Array & MeshData::getElementalDataArray(const ID & name, const ElementType & elem_type, const GhostType & ghost_type) const { auto it = elemental_data.find(name); if (it == elemental_data.end()) { AKANTU_EXCEPTION("Data named " << name << " not registered for type: " << elem_type << " - ghost_type:" << ghost_type << "!"); } return dynamic_cast &>(*(it->second))( elem_type, ghost_type); } template Array & MeshData::getElementalDataArray(const ID & name, const ElementType & elem_type, const GhostType & ghost_type) { auto it = elemental_data.find(name); if (it == elemental_data.end()) { AKANTU_EXCEPTION("Data named " << name << " not registered for type: " << elem_type << " - ghost_type:" << ghost_type << "!"); } return dynamic_cast &>(*(it->second))(elem_type, ghost_type); } /* -------------------------------------------------------------------------- */ // Get an elemental data array, if it does not exist: allocate it template Array & MeshData::getElementalDataArrayAlloc(const ID & name, const ElementType & elem_type, const GhostType & ghost_type, UInt nb_component) { auto it = elemental_data.find(name); ElementTypeMapArray * dataset; if (it == elemental_data.end()) { dataset = allocElementalData(name); } else { dataset = dynamic_cast *>(it->second); } AKANTU_DEBUG_ASSERT( getTypeCode() == typecode_map.find(name)->second, "Function getElementalDataArrayAlloc called with the wrong type!"); if (!(dataset->exists(elem_type, ghost_type))) { dataset->alloc(0, nb_component, elem_type, ghost_type); } return (*dataset)(elem_type, ghost_type); } /* -------------------------------------------------------------------------- */ #define AKANTU_MESH_DATA_CASE_MACRO(r, name, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ nb_comp = getNbComponentTemplated( \ name, el_type, ghost_type); \ break; \ } inline UInt MeshData::getNbComponent(const ID & name, const ElementType & el_type, const GhostType & ghost_type) const { auto it = typecode_map.find(name); UInt nb_comp(0); if (it == typecode_map.end()) { AKANTU_EXCEPTION("Could not determine the type held in dataset " << name << " for type: " << el_type << " - ghost_type:" << ghost_type << "."); } MeshDataTypeCode type = it->second; switch (type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_MESH_DATA_CASE_MACRO, name, AKANTU_MESH_DATA_TYPES) default: AKANTU_DEBUG_ERROR( "Could not call the correct instance of getNbComponentTemplated."); break; } return nb_comp; } #undef AKANTU_MESH_DATA_CASE_MACRO /* -------------------------------------------------------------------------- */ template inline UInt MeshData::getNbComponentTemplated(const ID & name, const ElementType & el_type, const GhostType & ghost_type) const { return getElementalDataArray(name, el_type, ghost_type).getNbComponent(); } /* -------------------------------------------------------------------------- */ // get the names of the data stored in elemental_data #define AKANTU_MESH_DATA_CASE_MACRO(r, name, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ ElementTypeMapArray * dataset; \ dataset = \ dynamic_cast *>( \ it->second); \ exists = dataset->exists(el_type, ghost_type); \ break; \ } inline void MeshData::getTagNames(StringVector & tags, const ElementType & el_type, const GhostType & ghost_type) const { tags.clear(); bool exists(false); auto it = elemental_data.begin(); auto it_end = elemental_data.end(); for (; it != it_end; ++it) { MeshDataTypeCode type = getTypeCode(it->first); switch (type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_MESH_DATA_CASE_MACRO, , AKANTU_MESH_DATA_TYPES) default: AKANTU_DEBUG_ERROR( "Could not determine the proper type to (dynamic-)cast."); break; } if (exists) { tags.push_back(it->first); } } } #undef AKANTU_MESH_DATA_CASE_MACRO /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_MESH_DATA_TMPL_HH__ */ diff --git a/src/mesh/mesh_inline_impl.cc b/src/mesh/mesh_inline_impl.cc index 3ecf0af8b..e78c5748c 100644 --- a/src/mesh/mesh_inline_impl.cc +++ b/src/mesh/mesh_inline_impl.cc @@ -1,636 +1,641 @@ /** * @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: Thu Jan 21 2016 * * @brief Implementation of the inline functions of the mesh class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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" #include "aka_iterators.hh" #include "element_class.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_INLINE_IMPL_CC__ #define __AKANTU_MESH_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline ElementKind Element::kind() const { return Mesh::getKind(type); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template Mesh::ElementTypesIteratorHelper Mesh::elementTypes(pack &&... _pack) const { return connectivities.elementTypes(_pack...); } /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh) : new_numbering(mesh.getNbNodes(), 1, "new_numbering") {} /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh, const ID & new_numbering_id) : new_numbering(new_numbering_id, mesh.getID(), mesh.getMemoryID()) {} /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewElementsEvent & event) { this->nodes_to_elements.resize(nodes->size()); for (const auto & elem : event.getList()) { const Array & conn = connectivities(elem.type, elem.ghost_type); UInt nb_nodes_per_elem = this->getNbNodesPerElement(elem.type); for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = conn(elem.element, n); if (not nodes_to_elements[node]) nodes_to_elements[node] = std::make_unique>(); nodes_to_elements[node]->insert(elem); } } EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewNodesEvent & event) { this->computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedElementsEvent & event) { this->connectivities.onElementsRemoved(event.getNewNumbering()); this->fillNodesToElements(); this->computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedNodesEvent & event) { const auto & new_numbering = event.getNewNumbering(); this->removeNodesFromArray(*nodes, new_numbering); if (nodes_global_ids) this->removeNodesFromArray(*nodes_global_ids, new_numbering); if (nodes_type.size() != 0) this->removeNodesFromArray(nodes_type, new_numbering); if (not nodes_to_elements.empty()) { std::vector>> tmp( nodes_to_elements.size()); auto it = nodes_to_elements.begin(); UInt new_nb_nodes = 0; for (auto new_i : new_numbering) { if (new_i != UInt(-1)) { tmp[new_i] = std::move(*it); ++new_nb_nodes; } ++it; } tmp.resize(new_nb_nodes); std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin()); } computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template inline void Mesh::removeNodesFromArray(Array & vect, const Array & new_numbering) { Array tmp(vect.size(), vect.getNbComponent()); UInt nb_component = vect.getNbComponent(); UInt new_nb_nodes = 0; for (UInt i = 0; i < new_numbering.size(); ++i) { UInt new_i = new_numbering(i); if (new_i != UInt(-1)) { T * to_copy = vect.storage() + i * nb_component; std::uninitialized_copy(to_copy, to_copy + nb_component, tmp.storage() + new_i * nb_component); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if (not nodes_global_ids) { nodes_global_ids = std::make_unique>( 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 (nodes_type.size() == 0) { nodes_type.resize(nodes->size()); nodes_type.set(_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); array.set(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::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); } /* -------------------------------------------------------------------------- */ 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 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(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type) const { UInt * conn_val = getConnectivity(type, ghost_type).storage(); UInt nb_nodes_per_element = getNbNodesPerElement(type); auto * local_coord = new Real[spatial_dimension * nb_nodes_per_element]; UInt offset = element * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { memcpy(local_coord + n * spatial_dimension, nodes->storage() + conn_val[offset + n] * spatial_dimension, spatial_dimension * sizeof(Real)); } Math::barycenter(local_coord, nb_nodes_per_element, spatial_dimension, barycenter); delete[] local_coord; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(const Element & element, Vector & barycenter) const { getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type); } /* -------------------------------------------------------------------------- */ 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; } /* -------------------------------------------------------------------------- */ 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.size() ? (nodes_type(n) == _nt_pure_gost) : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return nodes_type.size() ? (nodes_type(n) == _nt_master) || (nodes_type(n) == _nt_normal) : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return nodes_type.size() ? nodes_type(n) == _nt_normal : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return nodes_type.size() ? nodes_type(n) == _nt_master : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return nodes_type.size() ? nodes_type(n) >= 0 : false; } /* -------------------------------------------------------------------------- */ inline NodeType Mesh::getNodeType(UInt local_id) const { return nodes_type.size() ? nodes_type(local_id) : _nt_normal; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeLocalId(UInt global_id) const { AKANTU_DEBUG_ASSERT(nodes_global_ids != nullptr, "The global ids are note set."); return nodes_global_ids->find(global_id); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->size(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElementList(const Array & elements) { UInt nb_nodes_per_element = 0; UInt nb_nodes = 0; ElementType current_element_type = _not_defined; Array::const_iterator el_it = elements.begin(); Array::const_iterator el_end = elements.end(); for (; el_it != el_end; ++el_it) { const Element & el = *el_it; if (el.type != current_element_type) { current_element_type = el.type; nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type); } nb_nodes += nb_nodes_per_element; } return nb_nodes; } /* -------------------------------------------------------------------------- */ inline Mesh & Mesh::getMeshFacets() { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshFacets() const { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshParent() const { if (!this->mesh_parent) AKANTU_SILENT_EXCEPTION( "No parent mesh is defined! This is only valid in a mesh_facets"); return *this->mesh_parent; } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_MESH_INLINE_IMPL_CC__ */ diff --git a/src/mesh_utils/cohesive_element_inserter.cc b/src/mesh_utils/cohesive_element_inserter.cc index d71b9bfe9..2c9af99e0 100644 --- a/src/mesh_utils/cohesive_element_inserter.cc +++ b/src/mesh_utils/cohesive_element_inserter.cc @@ -1,394 +1,279 @@ /** * @file cohesive_element_inserter.cc * * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Sun Oct 04 2015 * * @brief Cohesive element inserter functions * * @section LICENSE * * Copyright (©) 2014, 2015 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 "cohesive_element_inserter.hh" #include "communicator.hh" #include "element_group.hh" #include "global_ids_updater.hh" #include "mesh_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { -CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, bool is_extrinsic, - const ID & id) +CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, const ID & id) : Parsable(ParserType::_cohesive_inserter), id(id), mesh(mesh), mesh_facets(mesh.initMeshFacets()), insertion_facets("insertion_facets", id), insertion_limits(mesh.getSpatialDimension(), 2), check_facets("check_facets", id) { - this->registerParam("groups", physical_groups, _pat_parsable, + this->registerParam("cohesive_surfaces", physical_groups, _pat_parsable, "List of groups to consider for insertion"); this->registerParam("bounding_box", insertion_limits, _pat_parsable, "Global limit for insertion"); - init(is_extrinsic); + UInt spatial_dimension = mesh.getSpatialDimension(); + + MeshUtils::resetFacetToDouble(mesh_facets); + + /// init insertion limits + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + insertion_limits(dim, 0) = std::numeric_limits::max() * Real(-1.); + insertion_limits(dim, 1) = std::numeric_limits::max(); + } + + insertion_facets.initialize(mesh_facets, + _spatial_dimension = spatial_dimension - 1, + _with_nb_element = true, _default_value = false); + } /* -------------------------------------------------------------------------- */ CohesiveElementInserter::~CohesiveElementInserter() = default; /* -------------------------------------------------------------------------- */ -void CohesiveElementInserter::init(bool is_extrinsic) { - AKANTU_DEBUG_IN(); - - UInt spatial_dimension = mesh.getSpatialDimension(); +void CohesiveElementInserter::parseSection(const ParserSection & section) { + Parsable::parseSection(section); - MeshUtils::resetFacetToDouble(mesh_facets); + if (is_extrinsic) + limitCheckFacets(this->check_facets); +} - /// initialize facet insertion array - insertion_facets.initialize(mesh_facets, - _spatial_dimension = (spatial_dimension - 1), - _with_nb_element = true); +/* -------------------------------------------------------------------------- */ +void CohesiveElementInserter::limitCheckFacets() { + limitCheckFacets(this->check_facets); +} - /// init insertion limits - for (UInt dim = 0; dim < spatial_dimension; ++dim) { - insertion_limits(dim, 0) = std::numeric_limits::max() * (-1.); - insertion_limits(dim, 1) = std::numeric_limits::max(); - } +/* -------------------------------------------------------------------------- */ +void CohesiveElementInserter::setLimit(SpacialDirection axis, Real first_limit, + Real second_limit) { + AKANTU_DEBUG_ASSERT( + axis < mesh.getSpatialDimension(), + "You are trying to limit insertion in a direction that doesn't exist"); - if (is_extrinsic) { - initFacetsCheck(); - } + insertion_limits(axis, 0) = std::min(first_limit, second_limit); + insertion_limits(axis, 1) = std::max(first_limit, second_limit); +} - AKANTU_DEBUG_OUT(); +/* -------------------------------------------------------------------------- */ +UInt CohesiveElementInserter::insertIntrinsicElements() { + limitCheckFacets(insertion_facets); + return insertElements(); } /* -------------------------------------------------------------------------- */ -void CohesiveElementInserter::initFacetsCheck() { +void CohesiveElementInserter::limitCheckFacets( + ElementTypeMapArray & check_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); check_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, - _with_nb_element = true, - _default_value = true); + _with_nb_element = true, _default_value = true); + check_facets.set(true); + // remove the pure ghost elements for_each_elements( mesh_facets, [&](auto && facet) { const auto & element_to_facet = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); auto & left = element_to_facet[0]; auto & right = element_to_facet[1]; if (right == ElementNull || (left.ghost_type == _ghost && right.ghost_type == _ghost)) { check_facets(facet) = false; + return; + } + + if (left.kind() == _ek_cohesive or right.kind() == _ek_cohesive) { + check_facets(facet) = false; } }, _spatial_dimension = spatial_dimension - 1); - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void CohesiveElementInserter::limitCheckFacets() { - AKANTU_DEBUG_IN(); - - UInt spatial_dimension = mesh.getSpatialDimension(); Vector bary_facet(spatial_dimension); - + // set the limits to the bounding box for_each_elements(mesh_facets, [&](auto && facet) { auto & need_check = check_facets(facet); if (not need_check) return; mesh_facets.getBarycenter(facet, bary_facet); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension && bary_facet(coord_in_limit) > insertion_limits(coord_in_limit, 0) && bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit != spatial_dimension) need_check = false; }, _spatial_dimension = spatial_dimension - 1); - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void CohesiveElementInserter::setLimit(SpacialDirection axis, Real first_limit, - Real second_limit) { - AKANTU_DEBUG_ASSERT( - axis < mesh.getSpatialDimension(), - "You are trying to limit insertion in a direction that doesn't exist"); - - insertion_limits(axis, 0) = std::min(first_limit, second_limit); - insertion_limits(axis, 1) = std::max(first_limit, second_limit); -} - -/* -------------------------------------------------------------------------- */ -UInt CohesiveElementInserter::insertIntrinsicElements() { - AKANTU_DEBUG_IN(); - - UInt spatial_dimension = mesh.getSpatialDimension(); - - Vector bary_facet(spatial_dimension); - - for_each_elements( - mesh_facets, - [&](auto && facet) { - const auto & element_to_facet = mesh_facets.getElementToSubelement( - facet.type, facet.ghost_type)(facet.element); - if (element_to_facet[1] == ElementNull) - return; - - mesh_facets.getBarycenter(facet, bary_facet); - UInt coord_in_limit = 0; - - while (coord_in_limit < spatial_dimension && - bary_facet(coord_in_limit) > - insertion_limits(coord_in_limit, 0) && - bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) - ++coord_in_limit; + if (not mesh_facets.hasData("physical_names")) { + AKANTU_DEBUG_OUT(); + return; + } - if (coord_in_limit == spatial_dimension) - insertion_facets(facet) = true; - }, - _spatial_dimension = spatial_dimension - 1); + const auto & physical_ids = + mesh_facets.getData("physical_names"); - AKANTU_DEBUG_OUT(); + // set the limits to the physical surfaces + for_each_elements(mesh_facets, + [&](auto && facet) { + auto & need_check = check_facets(facet); + if (not need_check) + return; - return insertElements(); -} + const auto & physical_id = physical_ids(facet); + auto it = find(physical_groups.begin(), + physical_groups.end(), physical_id); -/* -------------------------------------------------------------------------- */ -void CohesiveElementInserter::insertIntrinsicElements( - const std::string & physname, UInt material_index) { - AKANTU_DEBUG_IN(); + need_check = (it != physical_groups.end()); + }, + _spatial_dimension = spatial_dimension - 1); - UInt spatial_dimension = mesh.getSpatialDimension(); - ElementTypeMapArray * phys_data; - try { - phys_data = &(mesh_facets.getData("physical_names")); - } catch (...) { - phys_data = &(mesh_facets.registerData("physical_names")); - phys_data->initialize(mesh_facets, - _spatial_dimension = spatial_dimension - 1, - _with_nb_element = true); - // mesh_facets.initElementTypeMapArray(*phys_data, 1, spatial_dimension - - // 1, - // false, _ek_regular, true); - } - Vector bary_facet(spatial_dimension); - mesh_facets.createElementGroup(physname); - - GhostType ghost_type = _not_ghost; - - for (auto && type_facet : - mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { - auto & f_insertion = insertion_facets(type_facet, ghost_type); - auto & element_to_facet = - mesh_facets.getElementToSubelement(type_facet, ghost_type); - - UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type); - UInt coord_in_limit = 0; - - auto & group = mesh.getElementGroup(physname); - auto & group_facet = mesh_facets.getElementGroup(physname); - - Vector bary_physgroup(spatial_dimension); - Real norm_bary; - for (auto && e : group.getElements(type_facet, ghost_type)) { - mesh.getBarycenter(Element{type_facet, e, ghost_type}, bary_physgroup); -#ifndef AKANTU_NDEBUG - bool find_a_partner = false; -#endif - norm_bary = bary_physgroup.norm(); - Array & material_id = (*phys_data)(type_facet, ghost_type); - - for (UInt f = 0; f < nb_facet; ++f) { - - if (element_to_facet(f)[1] == ElementNull) - continue; - - mesh_facets.getBarycenter(Element{type_facet, f, ghost_type}, - bary_facet); - - coord_in_limit = 0; - - while (coord_in_limit < spatial_dimension && - (std::abs(bary_facet(coord_in_limit) - - bary_physgroup(coord_in_limit)) / - norm_bary < - Math::getTolerance())) - ++coord_in_limit; - - if (coord_in_limit == spatial_dimension) { - f_insertion(f) = true; -#ifndef AKANTU_NDEBUG - find_a_partner = true; -#endif - group_facet.add(type_facet, f, ghost_type, false); - material_id(f) = material_index; - break; - } - } - AKANTU_DEBUG_ASSERT(find_a_partner, - "The element nO " - << e << " of physical group " << physname - << " did not find its associated facet!" - << " Try to decrease math tolerance. " - << std::endl); - } - } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertElements(bool only_double_facets) { CohesiveNewNodesEvent node_event; NewElementsEvent element_event; Array new_pairs(0, 2); UInt nb_new_elements = MeshUtils::insertCohesiveElements( mesh, mesh_facets, insertion_facets, new_pairs, element_event.getList(), only_double_facets); UInt nb_new_nodes = new_pairs.size(); node_event.getList().reserve(nb_new_nodes); node_event.getOldNodesList().reserve(nb_new_nodes); for (UInt n = 0; n < nb_new_nodes; ++n) { node_event.getList().push_back(new_pairs(n, 1)); node_event.getOldNodesList().push_back(new_pairs(n, 0)); } if (mesh.isDistributed()) { /// update nodes type updateNodesType(mesh, node_event); updateNodesType(mesh_facets, node_event); /// update global ids nb_new_nodes = this->updateGlobalIDs(node_event); /// compute total number of new elements const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_elements, SynchronizerOperation::_sum); } if (nb_new_nodes > 0) mesh.sendEvent(node_event); if (nb_new_elements > 0) { updateInsertionFacets(); // mesh.updateTypesOffsets(_not_ghost); mesh.sendEvent(element_event); MeshUtils::resetFacetToDouble(mesh_facets); } if (mesh.isDistributed()) { /// update global ids this->synchronizeGlobalIDs(node_event); } return nb_new_elements; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::updateInsertionFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); - for (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - - GhostType facet_gt = *gt; - Mesh::type_iterator it = - mesh_facets.firstType(spatial_dimension - 1, facet_gt); - Mesh::type_iterator last = - mesh_facets.lastType(spatial_dimension - 1, facet_gt); - - for (; it != last; ++it) { - ElementType facet_type = *it; - - Array & ins_facets = insertion_facets(facet_type, facet_gt); - - // this is the extrinsic case - if (check_facets.exists(facet_type, facet_gt)) { - Array & f_check = check_facets(facet_type, facet_gt); - - UInt nb_facets = f_check.size(); - - for (UInt f = 0; f < ins_facets.size(); ++f) { - if (ins_facets(f)) { - ++nb_facets; - ins_facets(f) = false; - f_check(f) = false; - } - } - - f_check.resize(nb_facets); - } else { // and this the intrinsic one - ins_facets.resize(mesh_facets.getNbElement(facet_type, facet_gt)); - ins_facets.set(false); + for (auto && facet_gt : ghost_types) { + for (auto && facet_type : + mesh_facets.elementTypes(spatial_dimension - 1, facet_gt)) { + auto & ins_facets = insertion_facets(facet_type, facet_gt); + + // this is the intrinsic case + if (not is_extrinsic) + continue; + + auto & f_check = check_facets(facet_type, facet_gt); + for (auto && pair : zip(ins_facets, f_check)) { + bool & ins = std::get<0>(pair); + bool & check = std::get<1>(pair); + if (ins) + ins = check = false; } } } - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void CohesiveElementInserter::printself(std::ostream & stream, - int indent) const { - std::string space; - for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) - ; - - stream << space << "CohesiveElementInserter [" << std::endl; - - stream << space << " + mesh [" << std::endl; - mesh.printself(stream, indent + 2); - stream << space << AKANTU_INDENT << "]" << std::endl; + // resize for the newly added facets + insertion_facets.initialize(mesh_facets, + _spatial_dimension = spatial_dimension - 1, + _with_nb_element = true, _default_value = false); - stream << space << " + mesh_facets [" << std::endl; - mesh_facets.printself(stream, indent + 2); - stream << space << AKANTU_INDENT << "]" << std::endl; + // resize for the newly added facets + if (is_extrinsic) { + check_facets.initialize(mesh_facets, + _spatial_dimension = spatial_dimension - 1, + _with_nb_element = true, _default_value = false); + } else { + insertion_facets.set(false); + } - stream << space << "]" << std::endl; + AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/mesh_utils/cohesive_element_inserter.hh b/src/mesh_utils/cohesive_element_inserter.hh index b5bb140ce..4bc5d7b13 100644 --- a/src/mesh_utils/cohesive_element_inserter.hh +++ b/src/mesh_utils/cohesive_element_inserter.hh @@ -1,204 +1,187 @@ /** * @file cohesive_element_inserter.hh * * @author Fabian Barras * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Fri Oct 02 2015 * * @brief Cohesive element inserter * * @section LICENSE * * Copyright (©) 2014, 2015 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 "data_accessor.hh" #include "mesh_utils.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COHESIVE_ELEMENT_INSERTER_HH__ #define __AKANTU_COHESIVE_ELEMENT_INSERTER_HH__ namespace akantu { class GlobalIdsUpdater; class FacetSynchronizer; } // akantu namespace akantu { class CohesiveElementInserter : public DataAccessor, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - CohesiveElementInserter(Mesh & mesh, bool is_extrinsic = false, + CohesiveElementInserter(Mesh & mesh, const ID & id = "cohesive_element_inserter"); ~CohesiveElementInserter() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// init function - void init(bool is_extrinsic); - /// set range limitation for intrinsic cohesive element insertion void setLimit(SpacialDirection axis, Real first_limit, Real second_limit); /// insert intrinsic cohesive elements in a predefined range UInt insertIntrinsicElements(); - /// preset insertion of intrinsic cohesive elements along - /// a predefined group of facet and assign them a defined material index. - /// insertElement() method has to be called to finalize insertion. - void insertIntrinsicElements(const std::string & physname, - UInt material_index); - /// insert extrinsic cohesive elements (returns the number of new /// cohesive elements) UInt insertElements(bool only_double_facets = false); - /// function to print the contain of the class - virtual void printself(std::ostream & stream, int indent = 0) const; - /// limit check facets to match given insertion limits void limitCheckFacets(); /// init parallel variables void initParallel(ElementSynchronizer & element_synchronizer); + void parseSection(const ParserSection & section) override; + protected: - /// init facets check - void initFacetsCheck(); + /// internal version of limitCheckFacets + void limitCheckFacets(ElementTypeMapArray & check_facets); /// update facet insertion arrays after facets doubling void updateInsertionFacets(); /// update nodes type and global ids for parallel simulations /// (locally, within each processor) UInt updateGlobalIDs(NewNodesEvent & node_event); /// synchronize the global ids among the processors in parallel simulations void synchronizeGlobalIDs(NewNodesEvent & node_event); /// update nodes type void updateNodesType(Mesh & mesh, NewNodesEvent & node_event); /// functions for parallel communications 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; template inline void packUnpackGroupedInsertionData(CommunicationBuffer & buffer, const Array & elements) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_NOT_CONST(InsertionFacetsByElement, insertion_facets, ElementTypeMapArray &); - AKANTU_GET_MACRO(InsertionFacetsByElement, insertion_facets, const ElementTypeMapArray &); - AKANTU_GET_MACRO_BY_ELEMENT_TYPE(InsertionFacets, insertion_facets, bool); - AKANTU_GET_MACRO(CheckFacets, check_facets, const ElementTypeMapArray &); + AKANTU_GET_MACRO(CheckFacets, check_facets, + const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(CheckFacets, check_facets, bool); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(CheckFacets, check_facets, bool); + AKANTU_GET_MACRO(MeshFacets, mesh_facets, const Mesh &); AKANTU_GET_MACRO_NOT_CONST(MeshFacets, mesh_facets, Mesh &); + AKANTU_SET_MACRO(IsExtrinsic, is_extrinsic, bool); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// object id ID id; /// main mesh where to insert cohesive elements Mesh & mesh; /// mesh containing facets Mesh & mesh_facets; /// list of facets where to insert elements ElementTypeMapArray insertion_facets; /// limits for element insertion Matrix insertion_limits; /// list of groups to consider for insertion, ignored if empty std::vector physical_groups; /// vector containing facets in which extrinsic cohesive elements can be /// inserted ElementTypeMapArray check_facets; /// global connectivity ids updater std::unique_ptr global_ids_updater; -}; - -/* -------------------------------------------------------------------------- */ -/* inline functions */ -/* -------------------------------------------------------------------------- */ -/// standard output stream operator -inline std::ostream & operator<<(std::ostream & stream, - const CohesiveElementInserter & _this) { - _this.printself(stream); - return stream; -} + /// is this inserter used in extrinsic + bool is_extrinsic{false}; +}; class CohesiveNewNodesEvent : public NewNodesEvent { public: CohesiveNewNodesEvent() = default; ~CohesiveNewNodesEvent() override = default; AKANTU_GET_MACRO_NOT_CONST(OldNodesList, old_nodes, Array &); AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array &); private: Array old_nodes; }; } // akantu #include "cohesive_element_inserter_inline_impl.cc" #endif /* __AKANTU_COHESIVE_ELEMENT_INSERTER_HH__ */ diff --git a/src/model/solid_mechanics/material_selector.hh b/src/model/solid_mechanics/material_selector.hh index fd9d6d75a..ccbf3b1b4 100644 --- a/src/model/solid_mechanics/material_selector.hh +++ b/src/model/solid_mechanics/material_selector.hh @@ -1,159 +1,163 @@ /** * @file material_selector.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Thu Dec 17 2015 * * @brief class describing how to choose a material for a given element * * @section LICENSE * * Copyright (©) 2014, 2015 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" #include "element.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ - #ifndef __AKANTU_MATERIAL_SELECTOR_HH__ #define __AKANTU_MATERIAL_SELECTOR_HH__ /* -------------------------------------------------------------------------- */ - namespace akantu { class SolidMechanicsModel; /** * main class to assign same or different materials for different * elements */ -class MaterialSelector { +class MaterialSelector : public std::enable_shared_from_this { public: MaterialSelector() = default; virtual ~MaterialSelector() = default; virtual inline UInt operator()(const Element & element) { if (fallback_selector) return (*fallback_selector)(element); return fallback_value; } inline void setFallback(UInt f) { fallback_value = f; } inline void setFallback(const std::shared_ptr & fallback_selector) { this->fallback_selector = fallback_selector; } + + inline void + setFallback(MaterialSelector & fallback_selector) { + this->fallback_selector = fallback_selector.shared_from_this(); + } + inline std::shared_ptr & getFallbackSelector() { return this->fallback_selector; } inline UInt getFallbackValue() { return this->fallback_value; } protected: UInt fallback_value{0}; std::shared_ptr fallback_selector; }; /* -------------------------------------------------------------------------- */ /** * class that assigns the first material to regular elements by default */ class DefaultMaterialSelector : public MaterialSelector { public: explicit DefaultMaterialSelector( const ElementTypeMapArray & material_index) : material_index(material_index) {} UInt operator()(const Element & element) override { if (not material_index.exists(element.type, element.ghost_type)) return MaterialSelector::operator()(element); const auto & mat_indexes = material_index(element.type, element.ghost_type); if (element.element < mat_indexes.size()) { auto && tmp_mat = mat_indexes(element.element); if (tmp_mat != UInt(-1)) return tmp_mat; } return MaterialSelector::operator()(element); } private: const ElementTypeMapArray & material_index; }; /* -------------------------------------------------------------------------- */ /** * Use elemental data to assign materials */ template class ElementDataMaterialSelector : public MaterialSelector { public: ElementDataMaterialSelector(const ElementTypeMapArray & element_data, const SolidMechanicsModel & model, UInt first_index = 1) : element_data(element_data), model(model), first_index(first_index) {} inline T elementData(const Element & element) { DebugLevel dbl = debug::getDebugLevel(); debug::setDebugLevel(dblError); T data = element_data(element.type, element.ghost_type)(element.element); debug::setDebugLevel(dbl); return data; } inline UInt operator()(const Element & element) override { return MaterialSelector::operator()(element); } protected: /// list of element with the specified data (i.e. tag value) const ElementTypeMapArray & element_data; /// the model that the materials belong const SolidMechanicsModel & model; /// first material index: equal to 1 if none specified UInt first_index; }; /* -------------------------------------------------------------------------- */ /** * class to use mesh data information to assign different materials * where name is the tag value: tag_0, tag_1 */ template class MeshDataMaterialSelector : public ElementDataMaterialSelector { public: MeshDataMaterialSelector(const std::string & name, const SolidMechanicsModel & model, UInt first_index = 1); }; } // namespace akantu #endif /* __AKANTU_MATERIAL_SELECTOR_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 8f54e94e7..de016d862 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,562 +1,562 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the fem object needed for boundary conditions void initFEEngineBoundary(); protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::force vector (external forces) AKANTU_GET_MACRO(Force, *external_force, Array &); /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); - AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, material_selector, std::shared_ptr); + AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); AKANTU_SET_MACRO(MaterialSelector, material_selector, std::shared_ptr); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// lumped mass array Array * mass{nullptr}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force{nullptr}; /// internal forces array Array * internal_force{nullptr}; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs{nullptr}; /// array of current position used during update residual Array * current_position{nullptr}; UInt current_position_release{0}; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; #ifdef SWIGPYTHON public: #endif /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; #ifdef SWIGPYTHON protected: #endif /// class defining of to choose a material std::shared_ptr material_selector; /// flag defining if the increment must be computed or not bool increment_flag; /// tells if the material are instantiated bool are_materials_instantiated; using flatten_internal_map = std::map, ElementTypeMapArray *>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc index 605d44ae9..704ade675 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc @@ -1,162 +1,162 @@ /** * @file material_selector_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Dec 11 2015 * @date last modification: Mon Dec 14 2015 * * @brief Material selector for cohesive elements * * @section LICENSE * * Copyright (©) 2015 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 "material_selector_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DefaultMaterialCohesiveSelector::DefaultMaterialCohesiveSelector( const SolidMechanicsModelCohesive & model) : facet_material(model.getFacetMaterial()), mesh(model.getMesh()) { // backward compatibility v3: to get the former behavior back when the user // creates its own selector this->fallback_selector = std::make_shared(model.getMaterialByElement()); } /* -------------------------------------------------------------------------- */ UInt DefaultMaterialCohesiveSelector::operator()(const Element & element) { if (Mesh::getKind(element.type) == _ek_cohesive) { try { const Array & cohesive_el_to_facet = mesh.getMeshFacets().getSubelementToElement(element.type, element.ghost_type); bool third_dimension = (mesh.getSpatialDimension() == 3); const Element & facet = cohesive_el_to_facet(element.element, third_dimension); if (facet_material.exists(facet.type, facet.ghost_type)) { return facet_material(facet.type, facet.ghost_type)(facet.element); } else { return fallback_value; } } catch (...) { return fallback_value; } } else if (Mesh::getSpatialDimension(element.type) == mesh.getSpatialDimension() - 1) { return facet_material(element.type, element.ghost_type)(element.element); } else { return MaterialSelector::operator()(element); } } /* -------------------------------------------------------------------------- */ MeshDataMaterialCohesiveSelector::MeshDataMaterialCohesiveSelector( const SolidMechanicsModelCohesive & model) - : mesh_facets(model.getMeshFacets()), - material_index(mesh_facets.getData("physical_names")) { + : model(model), mesh_facets(model.getMeshFacets()), + material_index(mesh_facets.getData("physical_names")) { third_dimension = (model.getSpatialDimension() == 3); // backward compatibility v3: to get the former behavior back when the user // creates its own selector this->fallback_selector = std::make_shared>("physical_names", model); } /* -------------------------------------------------------------------------- */ UInt MeshDataMaterialCohesiveSelector::operator()(const Element & element) { if (Mesh::getKind(element.type) == _ek_cohesive) { const auto & facet = mesh_facets.getSubelementToElement( element.type, element.ghost_type)(element.element, third_dimension); - auto material_id = material_index(facet); - if (material_id != UInt(-1)) - return material_id; - - else + try { + std::string material_name = this->material_index(facet); + return this->model.getMaterialIndex(material_name); + } catch(...) { return fallback_value; + } } else return MaterialSelector::operator()(element); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ MaterialCohesiveRulesSelector::MaterialCohesiveRulesSelector( const SolidMechanicsModelCohesive & model, const MaterialCohesiveRules & rules, ID mesh_data_id) // what we have here is the name of model and also // the name of different materials : model(model), mesh_data_id(std::move(mesh_data_id)), mesh(model.getMesh()), mesh_facets(model.getMeshFacets()), spatial_dimension(model.getSpatialDimension()), rules(rules) { // cohesive fallback this->fallback_selector = std::make_shared(model); // non cohesive fallback this->fallback_selector->setFallback( std::make_shared>("physical_names", model)); } /* -------------------------------------------------------------------------- */ UInt MaterialCohesiveRulesSelector::operator()(const Element & element) { if (mesh_facets.getSpatialDimension(element.type) == (spatial_dimension - 1)) { const std::vector & element_to_subelement = mesh_facets.getElementToSubelement(element.type, element.ghost_type)(element.element); // Array & facets_check = model.getFacetsCheck(); const Element & el1 = element_to_subelement[0]; const Element & el2 = element_to_subelement[1]; ID id1 = mesh.getData(mesh_data_id, el1.type, el1.ghost_type)(el1.element); ID id2 = id1; if (el2 != ElementNull) { id2 = mesh.getData(mesh_data_id, el2.type, el2.ghost_type)(el2.element); } auto rit = rules.find(std::make_pair(id1, id2)); if (rit == rules.end()) { rit = rules.find(std::make_pair(id2, id1)); } if (rit != rules.end()) { return model.getMaterialIndex(rit->second); } } return MaterialSelector::operator()(element); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh index ca6cb74c6..a5dd25995 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh @@ -1,97 +1,98 @@ /** * @file material_selector_cohesive.hh * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Dec 11 2015 * @date last modification: Mon Dec 14 2015 * * @brief Material selectors for cohesive elements * * @section LICENSE * * Copyright (©) 2015 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 "material_selector.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } namespace akantu { #ifndef __AKANTU_MATERIAL_SELECTOR_COHESIVE_HH__ #define __AKANTU_MATERIAL_SELECTOR_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ /** * class that assigns the first cohesive material by default to the * cohesive elements */ class DefaultMaterialCohesiveSelector : public MaterialSelector { public: DefaultMaterialCohesiveSelector(const SolidMechanicsModelCohesive & model); UInt operator()(const Element & element) override; private: const ElementTypeMapArray & facet_material; const Mesh & mesh; }; /* -------------------------------------------------------------------------- */ /// To be used with intrinsic elements inserted along mesh physical surfaces class MeshDataMaterialCohesiveSelector : public MaterialSelector { public: MeshDataMaterialCohesiveSelector(const SolidMechanicsModelCohesive & model); UInt operator()(const Element & element) override; protected: + const SolidMechanicsModelCohesive & model; const Mesh & mesh_facets; - const ElementTypeMapArray & material_index; + const ElementTypeMapArray & material_index; bool third_dimension; }; /// bulk1, bulk2 -> cohesive using MaterialCohesiveRules = std::map, ID>; /* -------------------------------------------------------------------------- */ class MaterialCohesiveRulesSelector : public MaterialSelector { public: MaterialCohesiveRulesSelector(const SolidMechanicsModelCohesive & model, const MaterialCohesiveRules & rules, ID mesh_data_id = "physical_names"); UInt operator()(const Element & element); private: const SolidMechanicsModelCohesive & model; ID mesh_data_id; const Mesh & mesh; const Mesh & mesh_facets; UInt spatial_dimension; MaterialCohesiveRules rules; }; #endif /* __AKANTU_MATERIAL_SELECTOR_COHESIVE_HH__ */ } // akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index bcd0a3548..02a5301da 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,745 +1,752 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Jan 13 2016 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "parser.hh" #include "shape_cohesive.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); this->model_type = ModelType::_solid_mechanics_model_cohesive; auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); - material_selector = tmp_material_selector; + this->material_selector = tmp_material_selector; #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements - cohesive_synchronizer = std::make_unique( + this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); - cohesive_synchronizer->split(synchronizer, [](auto && el) { + this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, _gst_material_id); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_boundary); } + this->inserter = std::make_unique( + this->mesh, id + ":cohesive_element_inserter"); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = dynamic_cast(options); this->is_extrinsic = smmc_options.extrinsic; - inserter = std::make_unique( - mesh, is_extrinsic, id + ":cohesive_element_inserter"); + inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = dynamic_cast(mesh_facets.getElementSynchronizer()); this->registerSynchronizer(synchronizer, _gst_smmc_facets); this->registerSynchronizer(synchronizer, _gst_smmc_facets_conn); synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique(mesh_facets); facet_stress_synchronizer->updateSchemes( [&](auto & scheme, auto & proc, auto & direction) { scheme.copy(const_cast(synchronizer) .getCommunications() .getScheme(proc, direction)); }); this->registerSynchronizer(*facet_stress_synchronizer, _gst_smmc_facets_stress); } inserter->initParallel(*cohesive_synchronizer); } ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material - UInt cohesive_index = 0; + UInt cohesive_index = UInt(-1); - while ((dynamic_cast(materials[cohesive_index].get()) == - nullptr) && - cohesive_index <= materials.size()) - ++cohesive_index; + for(auto && material : enumerate(materials)) { + if(dynamic_cast(std::get<1>(material).get())){ + cohesive_index = std::get<0>(material); + break; + } + } - AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(), - "No cohesive materials in the material input file"); + if(cohesive_index == UInt(-1)) + AKANTU_EXCEPTION("No cohesive materials in the material input file"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks if (is_extrinsic) { this->initExtrinsicMaterials(); } else { this->initIntrinsicMaterials(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initExtrinsicMaterials() { const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = material_selector->getFallbackValue()); for_each_elements(mesh_facets, [&](auto && element) { auto mat_index = (*material_selector)(element); auto & mat = dynamic_cast( *materials[mat_index]); facet_material(element) = mat_index; mat.addFacet(element); }, _spatial_dimension = spatial_dimension - 1); SolidMechanicsModel::initMaterials(); this->initAutomaticInsertion(); } /* -------------------------------------------------------------------------- */ #if 0 void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials( const std::string & cohesive_surfaces) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != nullptr) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif std::istringstream split(cohesive_surfaces); std::string physname; while (std::getline(split, physname, ',')) { AKANTU_DEBUG_INFO( "Pre-inserting cohesive elements along facets from physical surface: " << physname); insertElementsFromMeshData(physname); } synchronizeInsertionData(); SolidMechanicsModel::initMaterials(); auto && tmp = std::make_shared(*this); tmp->setFallback(material_selector->getFallbackValue()); tmp->setFallback(material_selector->getFallbackSelector()); material_selector = tmp; // UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::synchronizeInsertionData() { if (inserter->getMeshFacets().isDistributed()) { inserter->getMeshFacets().getElementSynchronizer().synchronizeOnce( *inserter, _gst_ce_groups); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicMaterials() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initMaterials(); - inserter->insertElements(); + + this->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.size() == 0) continue; type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateFacetStressSynchronizer() { if (facet_stress_synchronizer != nullptr) { const auto & rank_to_element = mesh.getElementSynchronizer().getElementToRank(); const auto & facet_checks = inserter->getCheckFacets(); const auto & element_to_facet = mesh.getElementToSubelement(); UInt rank = mesh.getCommunicator().whoAmI(); facet_stress_synchronizer->updateSchemes( [&](auto & scheme, auto & proc, auto & /*direction*/) { UInt el = 0; for (auto && element : scheme) { if (not facet_checks(element)) continue; const auto & next_el = element_to_facet(element); UInt rank_left = rank_to_element(next_el[0]); UInt rank_right = rank_to_element(next_el[1]); if ((rank_left == rank and rank_right == proc) or (rank_left == proc and rank_right == rank)) { scheme[el] = element; ++el; } } scheme.resize(el); }); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); + this->inserter->limitCheckFacets(); + this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element Array & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); Array & el_q_facet = elements_quad_facets(type, elem_gt); ElementType facet_type = Mesh::getFacetType(type); UInt nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element global_facet_elem = facet_to_element(el, f); UInt global_facet = global_facet_elem.element; GhostType facet_gt = global_facet_elem.ghost_type; const Array & quad_f = quad_facets(facet_type, facet_gt); for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < Model::spatial_dimension; ++s) { el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s) = quad_f(global_facet * nb_quad_per_facet + q, s); } } } } } } /// loop over non cohesive materials for (auto && material : materials) { if (dynamic_cast(material.get())) continue; /// initialize the interpolation function material->initElementalFieldInterpolation(elements_quad_facets); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { auto & mat = dynamic_cast(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); if (isDefaultSolverExplicit()) { for (auto & material : materials) { try { auto & mat = dynamic_cast(*material); mat.computeEnergies(); } catch (std::bad_cast & bce) { } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { try { auto & mat __attribute__((unused)) = dynamic_cast(*material); } catch (std::bad_cast &) { /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } } #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData( debug::_dm_model_cohesive, "Interpolated stresses before", facet_stress); #endif this->synchronize(_gst_smmc_facets_stress); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses", facet_stress); #endif } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { interpolateStress(); for (auto & mat : materials) { try { auto & mat_cohesive = dynamic_cast(*mat); /// check which not ghost cohesive elements are to be created mat_cohesive.checkInsertion(); } catch (std::bad_cast &) { } } /// communicate data among processors this->synchronize(_gst_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); updateCohesiveSynchronizers(); SolidMechanicsModel::onElementsAdded(element_list, event); if (cohesive_synchronizer != nullptr) cohesive_synchronizer->computeAllBufferSizes(*this); if (is_extrinsic) resizeFacetStress(); /// if (method != _explicit_lumped_mass) { /// this->initSolver(); /// } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); // Array nodes_list(nb_new_nodes); // for (UInt n = 0; n < nb_new_nodes; ++n) // nodes_list(n) = doubled_nodes(n, 1); SolidMechanicsModel::onNodesAdded(new_nodes, event); UInt new_node, old_node; try { const auto & cohesive_event = dynamic_cast(event); const auto & old_nodes = cohesive_event.getOldNodesList(); auto copy = [this, &new_node, &old_node](auto & arr) { for (UInt s = 0; s < spatial_dimension; ++s) { arr(new_node, s) = arr(old_node, s); } }; for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; copy(*displacement); copy(*blocked_dofs); if (velocity) copy(*velocity); if (acceleration) copy(*acceleration); if (current_position) copy(*current_position); if (previous_displacement) copy(*previous_displacement); } // if (this->getDOFManager().hasMatrix("M")) { // this->assembleMass(old_nodes); // } // if (this->getDOFManager().hasLumpedMatrix("M")) { // this->assembleMassLumped(old_nodes); // } } catch (std::bad_cast &) { } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod &) { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ for (auto & mat : materials) { if (mat->isFiniteDeformation()) mat->computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); // for (auto && ghost_type : ghost_types) { // for (const auto & type : // mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { // UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); // UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") // .getNbIntegrationPoints(type, // ghost_type); // UInt new_size = nb_facet * nb_quadrature_points; // facet_stress(type, ghost_type).resize(new_size); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_material.cc b/src/model/solid_mechanics/solid_mechanics_model_material.cc index 1fcec1fff..ec3057ed4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_material.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_material.cc @@ -1,242 +1,254 @@ /** * @file solid_mechanics_model_material.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Nov 26 2010 * @date last modification: Mon Nov 16 2015 * * @brief instatiation of materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_factory.hh" #include "aka_math.hh" #include "material_non_local.hh" #include "mesh_iterators.hh" #include "non_local_manager.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material & SolidMechanicsModel::registerNewMaterial(const ParserSection & section) { std::string mat_name; std::string mat_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); mat_name = tmp; /** this can seam weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_DEBUG_ERROR( "A material of type \'" << mat_type << "\' in the input file has been defined without a name!"); } Material & mat = this->registerNewMaterial(mat_name, mat_type, opt_param); mat.parseSection(section); return mat; } /* -------------------------------------------------------------------------- */ Material & SolidMechanicsModel::registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(materials_names_to_id.find(mat_name) == materials_names_to_id.end(), "A material with this name '" << mat_name << "' has already been registered. " << "Please use unique names for materials"); UInt mat_count = materials.size(); materials_names_to_id[mat_name] = mat_count; std::stringstream sstr_mat; sstr_mat << this->id << ":" << mat_count << ":" << mat_type; ID mat_id = sstr_mat.str(); std::unique_ptr material = MaterialFactory::getInstance().allocate( mat_type, spatial_dimension, opt_param, *this, mat_id); materials.push_back(std::move(material)); return *(materials.back()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::instantiateMaterials() { - auto sub_sections = this->parser.getSubSections(ParserType::_material); + ParserSection model_section; + bool is_empty; + std::tie(model_section, is_empty) = this->getParserSection(); + + if(not is_empty) { + auto model_materials = model_section.getSubSections(ParserType::_material); + for (const auto & section : model_materials) { + this->registerNewMaterial(section); + } + } + auto sub_sections = this->parser.getSubSections(ParserType::_material); for (const auto & section : sub_sections) { - registerNewMaterial(section); + this->registerNewMaterial(section); } + + #ifdef AKANTU_DAMAGE_NON_LOCAL for (auto & material : materials) { if (dynamic_cast(material.get()) == nullptr) continue; this->non_local_manager = std::make_unique( *this, *this, id + ":non_local_manager", memory_id); break; } #endif are_materials_instantiated = true; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assignMaterialToElements( const ElementTypeMapArray * filter) { for_each_elements( mesh, [&](auto && element) { UInt mat_index = (*material_selector)(element); AKANTU_DEBUG_ASSERT( mat_index < materials.size(), "The material selector returned an index that does not exists"); material_index(element) = mat_index; }, _element_filter = filter, _ghost_type = _not_ghost); if (non_local_manager) non_local_manager->synchronize(*this, _gst_material_id); for_each_elements( mesh, [&](auto && element) { auto mat_index = material_index(element); auto index = materials[mat_index]->addElement(element); material_local_numbering(element) = index; }, _element_filter = filter, _ghost_type = _not_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initMaterials() { AKANTU_DEBUG_ASSERT(materials.size() != 0, "No material to initialize !"); if (!are_materials_instantiated) instantiateMaterials(); this->assignMaterialToElements(); // synchronize the element material arrays this->synchronize(_gst_material_id); for (auto & material : materials) { /// init internals properties material->initMaterial(); } this->synchronize(_gst_smm_init_mat); if (this->non_local_manager) { this->non_local_manager->initialize(); } } /* -------------------------------------------------------------------------- */ Int SolidMechanicsModel::getInternalIndexFromID(const ID & id) const { AKANTU_DEBUG_IN(); auto it = materials.begin(); auto end = materials.end(); for (; it != end; ++it) if ((*it)->getID() == id) { AKANTU_DEBUG_OUT(); return (it - materials.begin()); } AKANTU_DEBUG_OUT(); return -1; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::reassignMaterial() { AKANTU_DEBUG_IN(); std::vector> element_to_add(materials.size()); std::vector> element_to_remove(materials.size()); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & mat_indexes = material_index(type, ghost_type); for (UInt el = 0; el < nb_element; ++el) { element.element = el; UInt old_material = mat_indexes(el); UInt new_material = (*material_selector)(element); if (old_material != new_material) { element_to_add[new_material].push_back(element); element_to_remove[old_material].push_back(element); } } } } UInt mat_index = 0; for (auto mat_it = materials.begin(); mat_it != materials.end(); ++mat_it, ++mat_index) { (*mat_it)->removeElements(element_to_remove[mat_index]); (*mat_it)->addElements(element_to_add[mat_index]); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::applyEigenGradU( const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type) { AKANTU_DEBUG_ASSERT(prescribed_eigen_grad_u.size() == spatial_dimension * spatial_dimension, "The prescribed grad_u is not of the good size"); for (auto & material : materials) { if (material->getName() == material_name) material->applyEigenGradU(prescribed_eigen_grad_u, ghost_type); } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/3d_spherical_inclusion.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/3d_spherical_inclusion.geo index 6f87f12e1..ad6b42ee9 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/3d_spherical_inclusion.geo +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/3d_spherical_inclusion.geo @@ -1,100 +1,100 @@ -dx = 0.1; +dx = 1; w = 2; radius = 0.2; Point(1) = {0,0,0,dx}; Point(2) = {0,.3,0,dx}; Point(3) = {0,-.3,0,dx}; Point(4) = {w,0,0,dx}; Point(5) = {w,.3,0,dx}; Point(6) = {w,-.3,0,dx}; Point(7) = {0,0,w,dx}; Point(8) = {0,.3,w,dx}; Point(9) = {0,-.3,w,dx}; Point(10) = {w,0,w,dx}; Point(11) = {w,.3,w,dx}; Point(12) = {w,-.3,w,dx}; Point(13) = {0.5*w,0,0.5*w,dx}; Point(14) = {0.5*w+radius,0,0.5*w+radius,dx}; Point(15) = {0.5*w-radius,0,0.5*w+radius,dx}; Point(16) = {0.5*w+radius,0,0.5*w-radius,dx}; Point(17) = {0.5*w-radius,0,0.5*w-radius,dx}; Line(1) = {1, 2}; Line(2) = {2, 5}; Line(3) = {5, 4}; Line(4) = {1, 4}; Line(5) = {1, 3}; Line(6) = {6, 4}; Line(7) = {3, 6}; Line(8) = {8, 11}; Line(9) = {7, 10}; Line(10) = {9, 12}; Line(11) = {8, 7}; Line(12) = {11, 10}; Line(13) = {10, 12}; Line(14) = {7, 9}; Line(15) = {2, 8}; Line(16) = {1, 7}; Line(17) = {3, 9}; Line(18) = {5, 11}; Line(19) = {4, 10}; Line(20) = {6, 12}; Line Loop(21) = {18, 12, -19, -3}; Plane Surface(22) = {21}; Line Loop(23) = {19, 13, -20, 6}; Plane Surface(24) = {23}; Line Loop(25) = {11, -16, 1, 15}; Plane Surface(26) = {25}; Line Loop(27) = {14, -17, -5, 16}; Plane Surface(28) = {27}; Line Loop(29) = {8, 12, -9, -11}; Plane Surface(30) = {29}; Line Loop(31) = {9, 13, -10, -14}; Plane Surface(32) = {31}; Line Loop(33) = {2, 3, -4, 1}; Plane Surface(34) = {33}; Line Loop(35) = {7, 6, -4, 5}; Plane Surface(36) = {35}; Line Loop(37) = {18, -8, -15, 2}; Plane Surface(38) = {37}; Line Loop(39) = {7, 20, -10, -17}; Plane Surface(40) = {39}; Circle(41) = {16, 13, 14}; Circle(42) = {14, 13, 15}; Circle(43) = {15, 13, 17}; Circle(44) = {17, 13, 16}; Translate {0.5, 0, 0.5} { Duplicata { Line{44, 43, 42, 41}; } } Translate {-0.5, 0, 0.5} { Duplicata { Line{44, 43, 42, 41}; } } Translate {0.5, 0, -0.5} { Duplicata { Line{44, 43, 42, 41}; } } Translate {-0.5, 0, -0.5} { Duplicata { Line{44, 43, 42, 41}; } } Line Loop(61) = {9, -19, -4, 16}; Line Loop(62) = {52, 51, 50, 49}; Line Loop(63) = {57, 60, 59, 58}; Line Loop(64) = {44, 41, 42, 43}; Line Loop(65) = {55, 54, 53, 56}; Line Loop(66) = {48, 47, 46, 45}; Plane Surface(67) = {61, 62, 63, 64, 65, 66}; Plane Surface(68) = {62}; Plane Surface(69) = {63}; Plane Surface(70) = {64}; Plane Surface(71) = {66}; Plane Surface(72) = {65}; Surface Loop(73) = {26, 30, 38, 22, 34, 67, 68, 69, 70, 72, 71}; Volume(74) = {73}; Surface Loop(75) = {32, 24, 40, 36, 28, 67, 68, 69, 70, 72, 71}; Volume(76) = {75}; Physical Surface("interface") = {67}; Physical Surface("coh1") = {68}; Physical Surface("coh2") = {69}; Physical Surface("coh3") = {70}; Physical Surface("coh4") = {71}; Physical Surface("coh5") = {72}; Physical Volume("bulk") = {74, 76}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat index 4f6fc41e8..f39b45b53 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat @@ -1,65 +1,68 @@ -material elastic [ - name = bulk - rho = 2500 - nu = 0.29 - E = 70e9 -# finite_deformation = 1 -] +model solid_mechanics_model_cohesive [ + cohesive_inserter [ + cohesive_surfaces = [coh1, coh2, coh3, \ + coh4, coh5, interface] + ] -material cohesive_exponential [ - name = coh1 - sigma_c = 1.5e6 - beta = 1 - delta_c = 1e-4 - exponential_penalty = true - contact_tangent = 1.0 -] + material elastic [ + name = bulk + rho = 2500 + nu = 0.29 + E = 70e9 + # finite_deformation = 1 + ] -material cohesive_exponential [ - name = coh2 - sigma_c = 1.5e6 - beta = 1 - delta_c = 1e-4 - exponential_penalty = true - contact_tangent = 1.0 -] + material cohesive_exponential [ + name = coh1 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 + ] -material cohesive_exponential [ - name = coh3 - sigma_c = 1.5e6 - beta = 1 - delta_c = 1e-4 - exponential_penalty = true - contact_tangent = 1.0 -] + material cohesive_exponential [ + name = coh2 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 + ] -material cohesive_exponential [ - name = coh4 - sigma_c = 1.5e6 - beta = 1 - delta_c = 1e-4 - exponential_penalty = true - contact_tangent = 1.0 -] + material cohesive_exponential [ + name = coh3 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 + ] -material cohesive_exponential [ - name = coh5 - sigma_c = 1.5e6 - beta = 1 - delta_c = 1e-4 - exponential_penalty = true - contact_tangent = 1.0 -] + material cohesive_exponential [ + name = coh4 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 + ] -material cohesive_exponential [ - name = interface - sigma_c = 1.5e6 - beta = 1 - delta_c = 1e-4 - exponential_penalty = true - contact_tangent = 1.0 -] + material cohesive_exponential [ + name = coh5 + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 + ] -mesh parameters [ - cohesive_surfaces = coh1,coh2,coh3,coh4,coh5,interface + material cohesive_exponential [ + name = interface + sigma_c = 1.5e6 + beta = 1 + delta_c = 1e-4 + exponential_penalty = true + contact_tangent = 1.0 + ] ] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc index 7dd425c39..d06b84f83 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc @@ -1,107 +1,112 @@ /** * @file test_cohesive_insertion_along_physical_surfaces.cc * * @author Fabian Barras * * @date creation: Fri Aug 07 2015 * * @brief Test intrinsic insertion of cohesive elements along physical surfaces * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "material_cohesive.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("input_file.dat", argc, argv); Math::setTolerance(1e-15); const UInt spatial_dimension = 3; Mesh mesh(spatial_dimension); mesh.read("3d_spherical_inclusion.msh"); SolidMechanicsModelCohesive model(mesh); + + auto && material_selector = std::make_shared(model); + material_selector->setFallback(model.getMaterialSelector()); + model.setMaterialSelector(material_selector); + model.initFull(_analysis_method = _static); std::vector surfaces_name = {"interface", "coh1", "coh2", "coh3", "coh4", "coh5"}; UInt nb_surf = surfaces_name.size(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for (; it != end; ++it) { for (UInt i = 0; i < nb_surf; ++i) { UInt expected_insertion = mesh.getElementGroup(surfaces_name[i]) .getElements(mesh.getFacetType(*it)) .size(); UInt inserted_elements = model.getMaterial(surfaces_name[i]).getElementFilter()(*it).size(); AKANTU_DEBUG_ASSERT((expected_insertion == inserted_elements), std::endl << "!!! Mismatch in insertion of surface named " << surfaces_name[i] << " --> " << inserted_elements << " inserted elements out of " << expected_insertion << std::endl); } } /*std::string paraview_folder = "paraview/test_intrinsic_insertion_along_physical_surfaces/"; model.setDirectory(paraview_folder); model.setBaseName("bulk"); model.addDumpField("partitions"); model.dump(); model.setDirectoryToDumper("cohesive elements", paraview_folder); model.setBaseNameToDumper("cohesive elements", "one_cohesive_element"); model.addDumpFieldToDumper("cohesive elements", "partitions"); model.addDumpFieldToDumper("cohesive elements", "material_index"); model.dump("cohesive elements"); */ - model.assembleStiffnessMatrix(); + //model.assembleStiffnessMatrix(); - finalize(); + //finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc index eecab945d..8a5e96bf6 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc @@ -1,182 +1,180 @@ /** * @file test_cohesive_intrinsic.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Mon Jan 18 2016 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model_cohesive.hh" #include "dumper_paraview.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; static void updateDisplacement(SolidMechanicsModelCohesive &, Array &, ElementType, Real); int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); debug::setDebugLevel(dblWarning); const UInt spatial_dimension = 2; const UInt max_steps = 350; const ElementType type = _triangle_6; Mesh mesh(spatial_dimension); mesh.read("triangle.msh"); std::cout << mesh << std::endl; SolidMechanicsModelCohesive model(mesh); + model.getElementInserter().setLimit(_x, -0.26, -0.24); /// model initialization model.initFull(); - model.getElementInserter().setLimit(_x, -0.26, -0.24); - model.insertIntrinsicElements(); - mesh.write("mesh_cohesive.msh"); Real time_step = model.getStableTimeStep() * 0.8; model.setTimeStep(time_step); // std::cout << "Time step: " << time_step << std::endl; model.assembleMassLumped(); Array & boundary = model.getBlockedDOFs(); // const Array & residual = model.getResidual(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_element = mesh.getNbElement(type); /// boundary conditions for (UInt dim = 0; dim < spatial_dimension; ++dim) { for (UInt n = 0; n < nb_nodes; ++n) { boundary(n, dim) = true; } } model.assembleInternalForces(); model.setBaseName("intrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("internal_force"); model.addDumpField("stress"); model.addDumpField("strain"); model.addDumpField("external_force"); model.dump(); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_triangle"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump("cohesive elements"); /// update displacement Array elements; Vector bary(spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter({type, el, _not_ghost}, bary); if (bary(0) > -0.25) elements.push_back(el); } Real increment = 0.01; updateDisplacement(model, elements, type, increment); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); updateDisplacement(model, elements, type, increment); if (s % 1 == 0) { model.dump(); model.dump("cohesive elements"); std::cout << "passing step " << s << "/" << max_steps << ", Ed = " << model.getEnergy("dissipated") << std::endl; } } Real Ed = model.getEnergy("dissipated"); Real Edt = 2 * sqrt(2); std::cout << Ed << " " << Edt << std::endl; if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic was passed!" << std::endl; return EXIT_SUCCESS; } static void updateDisplacement(SolidMechanicsModelCohesive & model, Array & elements, ElementType type, Real increment) { Mesh & mesh = model.getFEEngine().getMesh(); UInt nb_element = elements.size(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type); Array & displacement = model.getDisplacement(); Array update(nb_nodes); update.clear(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity(elements(el), n); if (!update(node)) { displacement(node, 0) += increment; // displacement(node, 1) += increment; update(node) = true; } } } } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_quadrangle.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_quadrangle.cc index 296998f43..c90210373 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_quadrangle.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_quadrangle.cc @@ -1,221 +1,209 @@ /** * @file test_cohesive_intrinsic_quadrangle.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Thu Dec 11 2014 * * @brief Intrinsic cohesive elements' test for quadrangles * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; static void updateDisplacement(SolidMechanicsModelCohesive &, Array &, ElementType, Real); int main(int argc, char *argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 2; const UInt max_steps = 350; const ElementType type = _quadrangle_4; Mesh mesh(spatial_dimension); mesh.read("quadrangle.msh"); // debug::setDebugLevel(dblDump); // std::cout << mesh << std::endl; // debug::setDebugLevel(dblWarning); SolidMechanicsModelCohesive model(mesh); - + model.getElementInserter().setLimit(_x, -0.01, 0.01); /// model initialization model.initFull(); - model.getElementInserter().setLimit(_x, -0.01, 0.01); - model.insertIntrinsicElements(); - - // mesh.write("mesh_cohesive_quadrangle.msh"); - - // debug::setDebugLevel(dblDump); - // std::cout << mesh << std::endl; - // debug::setDebugLevel(dblWarning); - - Real time_step = model.getStableTimeStep()*0.8; model.setTimeStep(time_step); - // std::cout << "Time step: " << time_step << std::endl; - model.assembleMassLumped(); Array & boundary = model.getBlockedDOFs(); // const Array & residual = model.getResidual(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_element = mesh.getNbElement(type); /// boundary conditions for (UInt dim = 0; dim < spatial_dimension; ++dim) { for (UInt n = 0; n < nb_nodes; ++n) { boundary(n, dim) = true; } } model.assembleInternalForces(); model.setBaseName("intrinsic_quadrangle"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity" ); model.addDumpField("acceleration"); model.addDumpField("internal_force" ); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpField("external_force"); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_quadrangle"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump(); model.dump("cohesive elements"); /// update displacement Array elements; Vector bary(spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter({type, el, _not_ghost}, bary); if (bary(_x) > 0.) elements.push_back(el); } Real increment = 0.01; updateDisplacement(model, elements, type, increment); // for (UInt n = 0; n < nb_nodes; ++n) { // if (position(n, 1) + displacement(n, 1) > 0) { // if (position(n, 0) == 0) { // displacement(n, 1) -= 0.25; // } // if (position(n, 0) == 1) { // displacement(n, 1) += 0.25; // } // } // } // std::ofstream edis("edis.txt"); // std::ofstream erev("erev.txt"); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); updateDisplacement(model, elements, type, increment); if(s % 1 == 0) { model.dump(); model.dump("cohesive elements"); std::cout << "passing step " << s << "/" << max_steps << std::endl; } // // update displacement // for (UInt n = 0; n < nb_nodes; ++n) { // if (position(n, 1) + displacement(n, 1) > 0) { // displacement(n, 0) -= 0.01; // } // } // Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); // Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); // edis << s << " " // << Ed << std::endl; // erev << s << " " // << Er << std::endl; } // edis.close(); // erev.close(); Real Ed = model.getEnergy("dissipated"); Real Edt = 1; std::cout << Ed << " " << Edt << std::endl; if (Ed < Edt * 0.999 || Ed > Edt * 1.001) { std::cout << "The dissipated energy is incorrect" << std::endl; return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic_quadrangle was passed!" << std::endl; return EXIT_SUCCESS; } static void updateDisplacement(SolidMechanicsModelCohesive & model, Array & elements, ElementType type, Real increment) { Mesh & mesh = model.getFEEngine().getMesh(); UInt nb_element = elements.size(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type); Array & displacement = model.getDisplacement(); Array update(nb_nodes); update.clear(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity(elements(el), n); if (!update(node)) { displacement(node, 0) += increment; // displacement(node, 1) += increment; update(node) = true; } } } } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc index e833c8387..f81f1e0ef 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc @@ -1,355 +1,354 @@ /** * @file test_cohesive_intrinsic_tetrahedron.cc * * @author Marco Vocialta * * @date creation: Tue Aug 27 2013 * @date last modification: Thu Oct 15 2015 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; class Checker { public: Checker(const SolidMechanicsModelCohesive & model, const Array & elements, ElementType type); void check(const Vector & opening, const Matrix & rotation) { checkTractions(opening, rotation); checkEquilibrium(); computeEnergy(opening); } void updateDisplacement(const Vector & increment); protected: void checkTractions(const Vector & opening, const Matrix & rotation); void checkEquilibrium(); void checkResidual(const Matrix & rotation); void computeEnergy(const Vector & opening); private: std::set nodes_to_check; const SolidMechanicsModelCohesive & model; ElementType type; // const Array & elements; const Material & mat_cohesive; Real sigma_c; const Real beta; const Real G_c; const Real delta_0; const Real kappa; Real delta_c; const UInt spatial_dimension; const ElementType type_facet; const ElementType type_cohesive; const Array & traction; const Array & damage; const UInt nb_quad_per_el; const UInt nb_element; const Real beta2_kappa2; const Real beta2_kappa; Vector theoretical_traction; Vector traction_old; Vector opening_old; Real Ed; }; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("material_tetrahedron.dat", argc, argv); // debug::setDebugLevel(dblDump); const UInt spatial_dimension = 3; const UInt max_steps = 60; const Real increment_constant = 0.01; Math::setTolerance(1.e-12); const ElementType type = _tetrahedron_10; Mesh mesh(spatial_dimension); mesh.read("tetrahedron.msh"); SolidMechanicsModelCohesive model(mesh); + model.getElementInserter().setLimit(_x, -0.01, 0.01); /// model initialization model.initFull(); - model.getElementInserter().setLimit(_x, -0.01, 0.01); - model.insertIntrinsicElements(); Array & boundary = model.getBlockedDOFs(); boundary.set(true); UInt nb_element = mesh.getNbElement(type); model.setBaseName("intrinsic_tetrahedron"); model.addDumpFieldVector("displacement"); model.addDumpField("internal_force"); model.dump(); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_tetrahedron"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump("cohesive elements"); /// find elements to displace Array elements; Vector bary(spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter({type, el, _not_ghost}, bary); if (bary(_x) > 0.01) elements.push_back(el); } /// rotate mesh Real angle = 1.; // clang-format off Matrix rotation{ {std::cos(angle), std::sin(angle) * -1., 0.}, {std::sin(angle), std::cos(angle), 0.}, {0., 0., 1.}}; // clang-format on Vector increment_tmp{increment_constant, 2. * increment_constant, 3. * increment_constant}; Vector increment = rotation * increment_tmp; auto & position = mesh.getNodes(); auto position_it = position.begin(spatial_dimension); auto position_end = position.end(spatial_dimension); for (; position_it != position_end; ++position_it) { auto & pos = *position_it; pos = rotation * pos; } model.dump(); model.dump("cohesive elements"); /// find nodes to check Checker checker(model, elements, type); checker.updateDisplacement(increment); Real theoretical_Ed = 0; Vector opening(spatial_dimension, 0.); Vector opening_old(spatial_dimension, 0.); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); model.dump(); model.dump("cohesive elements"); opening += increment_tmp; checker.check(opening, rotation); checker.updateDisplacement(increment); } model.dump(); model.dump("cohesive elements"); Real Ed = model.getEnergy("dissipated"); theoretical_Ed *= 4.; std::cout << Ed << " " << theoretical_Ed << std::endl; if (!Math::are_float_equal(Ed, theoretical_Ed) || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; finalize(); return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic_tetrahedron was passed!" << std::endl; return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void Checker::updateDisplacement(const Vector & increment) { Mesh & mesh = model.getFEEngine().getMesh(); const auto & connectivity = mesh.getConnectivity(type); auto & displacement = model.getDisplacement(); Array update(displacement.size()); update.clear(); auto conn_it = connectivity.begin(connectivity.getNbComponent()); auto conn_end = connectivity.begin(connectivity.getNbComponent()); for (; conn_it != conn_end; ++conn_it) { const auto & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) { UInt node = conn(n); if (!update(node)) { Vector node_disp(displacement.storage() + node * spatial_dimension, spatial_dimension); node_disp += increment; update(node) = true; } } } } /* -------------------------------------------------------------------------- */ Checker::Checker(const SolidMechanicsModelCohesive & model, const Array & elements, ElementType type) : model(model), type(std::move(type)), // elements(elements), mat_cohesive(model.getMaterial(1)), sigma_c(mat_cohesive.get("sigma_c")), beta(mat_cohesive.get("beta")), G_c(mat_cohesive.get("G_c")), delta_0(mat_cohesive.get("delta_0")), kappa(mat_cohesive.get("kappa")), spatial_dimension(model.getSpatialDimension()), type_facet(Mesh::getFacetType(type)), type_cohesive(FEEngine::getCohesiveElementType(type_facet)), traction(mat_cohesive.getArray("tractions", type_cohesive)), damage(mat_cohesive.getArray("damage", type_cohesive)), nb_quad_per_el(model.getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive)), nb_element(model.getMesh().getNbElement(type_cohesive)), beta2_kappa2(beta * beta / kappa / kappa), beta2_kappa(beta * beta / kappa) { const Mesh & mesh = model.getMesh(); const auto & connectivity = mesh.getConnectivity(type); const auto & position = mesh.getNodes(); auto conn_it = connectivity.begin(connectivity.getNbComponent()); for (const auto & element : elements) { Vector conn_el(conn_it[element]); for (UInt n = 0; n < conn_el.size(); ++n) { UInt node = conn_el(n); if (Math::are_float_equal(position(node, _x), 0.)) nodes_to_check.insert(node); } } delta_c = 2 * G_c / sigma_c; sigma_c *= delta_c / (delta_c - delta_0); } /* -------------------------------------------------------------------------- */ void Checker::checkTractions(const Vector & opening, const Matrix & rotation) { auto normal_opening = opening * Vector{1., 0., 0.}; auto tangential_opening = opening - normal_opening; const Real normal_opening_norm = normal_opening.norm(); const Real tangential_opening_norm = tangential_opening.norm(); const Real delta = std::max(std::sqrt(tangential_opening_norm * tangential_opening_norm * beta2_kappa2 + normal_opening_norm * normal_opening_norm), 0.); Real theoretical_damage = std::min(delta / delta_c, 1.); theoretical_traction = (tangential_opening * beta2_kappa + normal_opening) * sigma_c / delta * (1. - theoretical_damage); // adjust damage theoretical_damage = std::max((delta - delta_0) / (delta_c - delta_0), 0.); theoretical_damage = std::min(theoretical_damage, 1.); Vector theoretical_traction_rotated = rotation * theoretical_traction; std::for_each( traction.begin(spatial_dimension), traction.end(spatial_dimension), [&theoretical_traction_rotated](auto && traction) { Real diff = Vector(theoretical_traction_rotated - traction).norm(); if (diff > 1e-14) throw std::domain_error("Tractions are incorrect"); }); std::for_each(damage.begin(), damage.end(), [&theoretical_damage](auto && damage) { if (not Math::are_float_equal(theoretical_damage, damage)) throw std::domain_error("Damage is incorrect"); }); } /* -------------------------------------------------------------------------- */ void Checker::computeEnergy(const Vector & opening) { /// compute energy auto Do = opening - opening_old; auto Dt = traction_old + theoretical_traction; Ed += .5 * Do.dot(Dt); opening_old = opening; traction_old = theoretical_traction; } /* -------------------------------------------------------------------------- */ void Checker::checkEquilibrium() { Vector residual_sum(spatial_dimension, 0.); const auto & residual = model.getInternalForce(); auto res_it = residual.begin(spatial_dimension); auto res_end = residual.end(spatial_dimension); for (; res_it != res_end; ++res_it) residual_sum += *res_it; if (!Math::are_float_equal(residual_sum.norm(), 0.)) throw std::domain_error("System is not in equilibrium!"); } /* -------------------------------------------------------------------------- */ void Checker::checkResidual(const Matrix & rotation) { Vector total_force(spatial_dimension, 0.); const auto & residual = model.getInternalForce(); for (auto node : nodes_to_check) { Vector res(residual.begin(spatial_dimension)[node]); total_force += res; } Vector theoretical_total_force(spatial_dimension); theoretical_total_force.mul(rotation, theoretical_traction); theoretical_total_force *= -1 * 2 * 2; for (UInt s = 0; s < spatial_dimension; ++s) { if (!Math::are_float_equal(total_force(s), theoretical_total_force(s))) { std::cout << "Total force isn't correct!" << std::endl; std::terminate(); } } } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron_fragmentation.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron_fragmentation.cc index 55a703625..5165b0dcd 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron_fragmentation.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron_fragmentation.cc @@ -1,127 +1,126 @@ /** * @file test_cohesive_intrinsic_tetrahedron_fragmentation.cc * * @author Marco Vocialta * * @date creation: Wed Oct 09 2013 * @date last modification: Thu Dec 11 2014 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); // debug::setDebugLevel(dblDump); ElementType type = _tetrahedron_10; const UInt spatial_dimension = 3; const UInt max_steps = 100; Mesh mesh(spatial_dimension); mesh.read("tetrahedron_full.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(); - model.insertIntrinsicElements(); Real time_step = model.getStableTimeStep() * 0.8; model.setTimeStep(time_step); // std::cout << "Time step: " << time_step << std::endl; model.assembleMassLumped(); model.assembleInternalForces(); model.setBaseName("intrinsic_tetrahedron_fragmentation"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("internal_force"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_tetrahedron_fragmentation"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump(); model.dump("cohesive elements"); /// update displacement UInt nb_element = mesh.getNbElement(type); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Vector bary(spatial_dimension); const Array & connectivity = mesh.getConnectivity(type); Array & displacement = model.getDisplacement(); Array update(nb_nodes); for (UInt s = 0; s < max_steps; ++s) { Real increment = s / 10.; update.clear(); for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter({type, el, _not_ghost}, bary); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity(el, n); if (!update(node)) { for (UInt dim = 0; dim < spatial_dimension; ++dim) { displacement(node, dim) = increment * bary(dim); update(node) = true; } } } } if (s % 10 == 0) { model.dump(); model.dump("cohesive elements"); } } if (nb_nodes != nb_element * Mesh::getNbNodesPerElement(type)) { std::cout << "Wrong number of nodes" << std::endl; finalize(); return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic_tetrahedron was passed!" << std::endl; return EXIT_SUCCESS; }