diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 1306d2f68..05f51fe94 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,712 +1,723 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Sat May 01 2021 * * @brief common type descriptions for akantu * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_COMMON_HH_ #define AKANTU_COMMON_HH_ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Constants */ /* -------------------------------------------------------------------------- */ namespace { [[gnu::unused]] constexpr Int _all_dimensions{ std::numeric_limits::max()}; #ifdef AKANTU_NDEBUG [[gnu::unused]] constexpr Real REAL_INIT_VALUE{0.}; #else [[gnu::unused]] constexpr Real REAL_INIT_VALUE{ std::numeric_limits::quiet_NaN()}; #endif } // namespace /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "aka_enum_macros.hh" /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of /// events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) \ (contact_mechanics_model) \ (coupler_solid_contact) \ (coupler_solid_cohesive_contact) \ (phase_field_model) \ (coupler_solid_phasefield) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { model, solid_mechanics_model, solid_mechanics_model_cohesive, heat_transfer_model, structural_mechanics_model, embedded_model, }; #endif /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3, _explicit_contact = 4, _implicit_contact = 5 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_NON_LINEAR_SOLVER_TYPES \ (linear) \ (newton_raphson) \ (newton_raphson_modified) \ (lumped) \ (gmres) \ (bfgs) \ (cg) \ (newton_raphson_contact) \ (auto) // clang-format on AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) #else /// Type of non linear resolution available in akantu enum class NonLinearSolverType { _linear, ///< No non linear convergence loop _newton_raphson, ///< Regular Newton-Raphson _newton_raphson_modified, ///< Newton-Raphson with initial tangent _lumped, ///< Case of lumped mass or equivalent matrix _gmres, _bfgs, _cg, _newton_raphson_contact, ///< Regular Newton-Raphson modified /// for contact problem _auto, ///< This will take a default value that make sense in case of /// model::getNewSolver }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_TIME_STEP_SOLVER_TYPE \ (static) \ (dynamic) \ (dynamic_lumped) \ (not_defined) // clang-format on AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) #else /// Type of time stepping solver enum class TimeStepSolverType { _static, ///< Static solution _dynamic, ///< Dynamic solver _dynamic_lumped, ///< Dynamic solver with lumped mass _not_defined, ///< For not defined cases }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_INTEGRATION_SCHEME_TYPE \ (pseudo_time) \ (forward_euler) \ (trapezoidal_rule_1) \ (backward_euler) \ (central_difference) \ (fox_goodwin) \ (trapezoidal_rule_2) \ (linear_acceleration) \ (newmark_beta) \ (generalized_trapezoidal) // clang-format on AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) #else /// Type of integration scheme enum class IntegrationSchemeType { _pseudo_time, ///< Pseudo Time _forward_euler, ///< GeneralizedTrapezoidal(0) _trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _backward_euler, ///< GeneralizedTrapezoidal(1) _central_difference, ///< NewmarkBeta(0, 1/2) _fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ (residual) \ (solution) \ (residual_mass_wgh) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) #else /// enum SolveConvergenceCriteria different convergence criteria enum class SolveConvergenceCriteria { _residual, ///< Use residual to test the convergence _solution, ///< Use solution to test the convergence _residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; #endif /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum MatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /// @enum Type of contact detection enum DetectionType { _explicit, _implicit }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_CONTACT_STATE \ (no_contact) \ (stick) \ (slip) // clang-format on AKANTU_CLASS_ENUM_DECLARE(ContactState, AKANTU_CONTACT_STATE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ContactState, AKANTU_CONTACT_STATE) AKANTU_CLASS_ENUM_INPUT_STREAM(ContactState, AKANTU_CONTACT_STATE) #else /// @enum no contact or stick or slip state enum class ContactState { _no_contact = 0, _stick = 1, _slip = 2, }; #endif /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_SYNCHRONIZATION_TAG \ (whatever) \ (update) \ (ask_nodes) \ (size) \ (smm_mass) \ (smm_for_gradu) \ (smm_boundary) \ (smm_uv) \ (smm_res) \ (smm_init_mat) \ (smm_stress) \ (smmc_facets) \ (smmc_facets_conn) \ (smmc_facets_stress) \ (smmc_damage) \ (giu_global_conn) \ (ce_groups) \ (ce_insertion_order) \ (gm_clusters) \ (htm_temperature) \ (htm_gradient_temperature) \ (htm_phi) \ (htm_gradient_phi) \ (pfm_damage) \ (pfm_driving) \ (pfm_history) \ (pfm_energy) \ (csp_damage) \ (csp_strain) \ (mnl_for_average) \ (mnl_weight) \ (nh_criterion) \ (test) \ (user_1) \ (user_2) \ (material_id) \ (for_dump) \ (cf_nodal) \ (cf_incr) \ (solver_solution) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) #else /// @enum SynchronizationTag type of synchronizations enum class SynchronizationTag { //--- Generic tags --- _whatever, _update, _ask_nodes, _size, //--- SolidMechanicsModel tags --- _smm_mass, ///< synchronization of the SolidMechanicsModel.mass _smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _smm_uv, ///< synchronization of the nodal velocities and displacement _smm_res, ///< synchronization of the nodal residual _smm_init_mat, ///< synchronization of the data to initialize materials _smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _smmc_facets, ///< synchronization of facet data to setup facet synch _smmc_facets_conn, ///< synchronization of facet global connectivity _smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups _ce_insertion_order, ///< synchronization of the order of insertion of /// cohesive elements // --- GroupManager tags --- _gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _htm_temperature, ///< synchronization of the nodal temperature _htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- PhaseFieldModel tags --- _pfm_damage, ///< synchronization of the nodal damage _pfm_driving, ///< synchronization of the driving forces to /// compute the internal _pfm_history, ///< synchronization of the damage history to /// compute the internal _pfm_energy, ///< synchronization of the damage energy /// density to compute the internal // --- CouplerSolidPhaseField tags --- _csp_damage, ///< synchronization of the damage from phase /// model to solid model _csp_strain, ///< synchronization of the strain from solid /// model to phase model // --- LevelSet tags --- _htm_phi, ///< synchronization of the nodal level set value phi _htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _mnl_for_average, ///< synchronization of data to average in non local /// material _mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _nh_criterion, // --- General tags --- _test, ///< Test tag _user_1, ///< tag for user simulations _user_2, ///< tag for user simulations _material_id, ///< synchronization of the material ids _for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _cf_nodal, ///< synchronization of disp, velo, and current position _cf_incr, ///< synchronization of increment // --- Solver tags --- _solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; #endif /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } std::ostream & operator<<(std::ostream & stream, NodeFlag flag); } // namespace akantu AKANTU_ENUM_HASH(GhostType) namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; namespace { constexpr ghost_type_t ghost_types{_casper}; } /// 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 /* -------------------------------------------------------------------------- */ #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_AUTO(name, variable) \ +#define AKANTU_GET_MACRO_AUTO(name, variable) \ inline decltype(auto) get##name() const { return (variable); } +#define AKANTU_GET_MACRO_AUTO_NOT_CONST(name, variable) \ + inline decltype(auto) get##name() { return (variable); } + #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_DEREF_PTR(name, ptr) \ inline const auto & get##name() const { \ if (not(ptr)) { \ AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ } \ return (*(ptr)); \ } #define AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(name, ptr) \ inline auto & get##name() { \ if (not(ptr)) { \ AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ } \ return (*(ptr)); \ } +#define AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(name, ptr) \ + inline decltype(auto) get##name() { \ + if (not ptr) { \ + AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ + } \ + return (*ptr); \ + } + #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name(const support & el_type, \ GhostType ghost_type = _not_ghost) \ con { /* NOLINT */ \ return variable(el_type, ghost_type); \ } // NOLINT #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ struct TensorTrait {}; struct TensorProxyTrait {}; } // namespace akantu /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ namespace aka { /* ------------------------------------------------------------------------ */ template using is_scalar = std::is_arithmetic; /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> bool is_of_type(T && t) { return ( dynamic_cast>::value, std::add_const_t, R>>>(&t) != nullptr); } /* -------------------------------------------------------------------------- */ template bool is_of_type(std::unique_ptr & t) { return ( dynamic_cast::value, std::add_const_t, R>>>( t.get()) != nullptr); } /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { static_assert( disjunction< std::is_base_of, std::decay_t>, // down-cast std::is_base_of, std::decay_t> // up-cast >::value, "Type T and R are not valid for a as_type conversion"); return dynamic_cast>::value, std::add_const_t, R>>>(t); } /* -------------------------------------------------------------------------- */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { return &as_type(*t); } /* -------------------------------------------------------------------------- */ template decltype(auto) as_type(const std::shared_ptr & t) { return std::dynamic_pointer_cast(t); } } // namespace aka #include "aka_common_inline_impl.hh" #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); #define AKANTU_CURRENT_FUNCTION \ (std::string(__func__) + "():" + std::to_string(__LINE__)) } // namespace akantu /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic * number is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif // AKANTU_COMMON_HH_ diff --git a/src/common/aka_element_classes_info_inline_impl.hh b/src/common/aka_element_classes_info_inline_impl.hh index 1407320ff..617b31c40 100644 --- a/src/common/aka_element_classes_info_inline_impl.hh +++ b/src/common/aka_element_classes_info_inline_impl.hh @@ -1,226 +1,226 @@ /** * @file aka_element_classes_info_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Thu Jun 18 2015 * @date last modification: Tue Sep 29 2020 * * @brief Implementation of the streaming fonction for the element classes * enums * * * @section LICENSE * * Copyright (©) 2015-2021 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 "aka_tuple_tools.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_HH_ #define AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ // BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING #define AKANTU_BOOST_CASE_MACRO(r, macro, _type) \ case _type: { \ macro(_type); \ break; \ } #define AKANTU_BOOST_LIST_SWITCH(macro1, list1, var) \ do { \ switch (var) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1) \ default: { \ - AKANTU_ERROR("Type (" << var << ") not handled by this function"); \ + AKANTU_ERROR("Type (" << var << ") not handled by this function"); \ } \ } \ } while (0) #define AKANTU_BOOST_LIST_SWITCH_NO_DEFAULT(macro1, list1, var) \ do { \ switch (var) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1) \ case _not_defined: \ break; \ case _max_element_type: \ break; \ } \ } while (0) #define AKANTU_BOOST_LIST_SWITCH_CONSTEXPR(macro1, list1, var) \ do { \ switch (var) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1) \ default: { \ - macro1(_not_defined); \ + macro1(_not_defined); \ } \ } \ } while (0) #define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1) \ AKANTU_BOOST_LIST_SWITCH(macro1, list1, type) #define AKANTU_BOOST_ELEMENT_SWITCH_NO_DEFAULT(macro1, list1) \ AKANTU_BOOST_LIST_SWITCH_NO_DEFAULT(macro1, list1, type) #define AKANTU_BOOST_ELEMENT_SWITCH_CONSTEXPR(macro1, list1) \ AKANTU_BOOST_LIST_SWITCH_CONSTEXPR(macro1, list1, type) #define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(macro) \ AKANTU_BOOST_ELEMENT_SWITCH_NO_DEFAULT(macro, AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_BOOST_ALL_ELEMENT_SWITCH_CONSTEXPR(macro) \ AKANTU_BOOST_ELEMENT_SWITCH_CONSTEXPR(macro, AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_BOOST_LIST_MACRO(r, macro, type) macro(type) #define AKANTU_BOOST_APPLY_ON_LIST(macro, list) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list) #define AKANTU_BOOST_ALL_ELEMENT_LIST(macro) \ AKANTU_BOOST_APPLY_ON_LIST(macro, AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_GET_ELEMENT_LIST(kind) AKANTU##kind##_ELEMENT_TYPE #define AKANTU_BOOST_KIND_ELEMENT_SWITCH(macro, kind) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, AKANTU_GET_ELEMENT_LIST(kind)) // BOOST_PP_SEQ_TO_LIST does not exists in Boost < 1.49 #define AKANTU_GENERATE_KIND_LIST(seq) \ BOOST_PP_TUPLE_TO_LIST(BOOST_PP_SEQ_SIZE(seq), BOOST_PP_SEQ_TO_TUPLE(seq)) #define AKANTU_ELEMENT_KIND_BOOST_LIST \ AKANTU_GENERATE_KIND_LIST(AKANTU_ELEMENT_KIND) #define AKANTU_BOOST_ALL_KIND_LIST(macro, list) \ BOOST_PP_LIST_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list) #define AKANTU_BOOST_ALL_KIND(macro) \ AKANTU_BOOST_ALL_KIND_LIST(macro, AKANTU_ELEMENT_KIND_BOOST_LIST) #define AKANTU_BOOST_ALL_KIND_SWITCH(macro) \ AKANTU_BOOST_LIST_SWITCH(macro, AKANTU_ELEMENT_KIND, kind) /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ AKANTU_ENUM_OUTPUT_STREAM( ElementType, AKANTU_ALL_ELEMENT_TYPE(_not_defined)(_max_element_type)) AKANTU_ENUM_INPUT_STREAM(ElementType, AKANTU_ALL_ELEMENT_TYPE) AKANTU_ENUM_OUTPUT_STREAM(InterpolationType, AKANTU_INTERPOLATION_TYPES) AKANTU_ENUM_INPUT_STREAM(InterpolationType, AKANTU_INTERPOLATION_TYPES) AKANTU_ENUM_OUTPUT_STREAM(ElementKind, AKANTU_ELEMENT_KIND) AKANTU_ENUM_INPUT_STREAM(ElementKind, AKANTU_ELEMENT_KIND) template <::akantu::ElementType t> using element_type_t = std::integral_constant<::akantu::ElementType, t>; #define OP_CAT(s, data, elem) BOOST_PP_CAT(_element_type, elem) // creating a type instead of a using helps to debug #define AKANTU_DECLARE_ELEMENT_TYPE_STRUCT(r, data, elem) \ struct BOOST_PP_CAT(_element_type, elem) \ : public element_type_t<::akantu::elem> {}; BOOST_PP_SEQ_FOR_EACH(AKANTU_DECLARE_ELEMENT_TYPE_STRUCT, _, AKANTU_ALL_ELEMENT_TYPE) #undef AKANTU_DECLARE_ELEMENT_TYPE_STRUCT template struct ElementTypes; template <> struct ElementTypes<_ek_regular> { using type = std::tuple; }; #if defined(AKANTU_COHESIVE_ELEMENT) template <> struct ElementTypes<_ek_cohesive> { using type = std::tuple; }; #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) template <> struct ElementTypes<_ek_structural> { using type = std::tuple; }; #endif #undef OP_CAT template using ElementTypes_t = typename ElementTypes::type; #define OP_CAT(s, data, elem) ElementTypes_t using AllElementTypes = tuple::cat_t; + BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ELEMENT_KIND))>; #undef OP_CAT - + namespace details { template struct visit_tuple_impl { template static constexpr decltype(auto) visit(const Tuple &, Function && function, const DynamicType & type) { using integral_type = std::tuple_element_t; if (integral_type::value == type) { return function(integral_type{}); } else { return visit_tuple_impl::visit( Tuple{}, std::forward(function), type); } } }; template <> struct visit_tuple_impl<0> { template - static constexpr decltype(auto) visit(const Tuple &, Function && function, - const DynamicType & /*type*/) { - return function(element_type_t<_not_defined>{}); + static constexpr auto + visit(const Tuple &, Function && function, const DynamicType & /*type*/) + -> decltype(function(std::tuple_element_t<0, Tuple>{})) { + throw; } }; } // namespace details template constexpr decltype(auto) tuple_dispatch(Function && function, const DynamicType & type) { return details::visit_tuple_impl::value>::visit( Tuple{}, std::forward(function), type); } - } // namespace akantu #endif /* AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_HH_ */ diff --git a/src/common/aka_view_iterators.hh b/src/common/aka_view_iterators.hh index 42daf6542..74628e51d 100644 --- a/src/common/aka_view_iterators.hh +++ b/src/common/aka_view_iterators.hh @@ -1,585 +1,597 @@ /** * @file aka_view_iterators.hh * * @author Nicolas Richart * * @date creation Thu Nov 15 2018 * * @brief View iterators * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" //#include "aka_types.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_VIEW_ITERATORS_HH__ #define __AKANTU_AKA_VIEW_ITERATORS_HH__ namespace akantu { template class TensorBase; } namespace akantu { /* -------------------------------------------------------------------------- */ /* Iterators */ /* -------------------------------------------------------------------------- */ namespace detail { template struct IteratorHelper { static constexpr Int dim = 0; }; template struct IteratorHelper> { private: using T = typename Derived::Scalar; static constexpr Int m = Derived::RowsAtCompileTime; static constexpr Int n = Derived::ColsAtCompileTime; public: static constexpr Int dim = Eigen::MatrixBase::IsVectorAtCompileTime ? 1 : 2; using pointer = T *; using proxy = Eigen::Map>; using const_proxy = Eigen::Map>; }; template struct IteratorHelper> { private: using T = typename Derived::Scalar; static constexpr Int m = Derived::RowsAtCompileTime; static constexpr Int n = Derived::ColsAtCompileTime; public: static constexpr Int dim = Derived::IsVectorAtCompileTime and m != 1 ? 1 : 2; using pointer = T *; using proxy = Eigen::Map; using const_proxy = Eigen::Map; }; template struct IteratorHelper> { static constexpr Int dim = _dim; using pointer = T *; using proxy = TensorProxy; using const_proxy = TensorProxy; }; template struct IteratorHelper> { static constexpr Int dim = _dim; using pointer = T *; using proxy = TensorProxy; using const_proxy = TensorProxy; }; /* -------------------------------------------------------------------------- */ template >::dim> class internal_view_iterator { protected: using helper = IteratorHelper>; using internal_value_type = IR; using internal_pointer = IR *; using scalar_pointer = typename helper::pointer; using proxy_t = typename helper::proxy; using const_proxy_t = typename helper::const_proxy; static constexpr int dim_ = dim; public: using pointer = proxy_t *; using value_type = proxy_t; using reference = proxy_t &; using const_reference = const_proxy_t &; using difference_type = Int; using iterator_category = std::random_access_iterator_tag; private: template constexpr auto get_new_const_proxy(scalar_pointer data, std::index_sequence) const { return const_proxy_t(data, dims[I]...); } constexpr auto get_new_const_proxy(scalar_pointer data) const { return get_new_const_proxy(data, std::make_index_sequence()); } template constexpr auto get_new_proxy(scalar_pointer data, std::index_sequence) { return proxy_t(data, dims[I]...); } constexpr auto get_new_proxy(scalar_pointer data) { return get_new_proxy(data, std::make_index_sequence()); } template constexpr void reset_proxy(T & t, scalar_pointer data, std::index_sequence) const { new (&t) T(data, dims[I]...); } template constexpr auto reset_proxy(T & t) const { return reset_proxy(t, this->ret_ptr, std::make_index_sequence()); } protected: template ().data()), decltype(std::declval().data())>::value> * = nullptr> explicit internal_view_iterator( internal_view_iterator & it) : dims(it.dims), _offset(it._offset), initial(it.initial), ret_ptr(it.ret_ptr), proxy(get_new_proxy(it.ret_ptr)), const_proxy(get_new_const_proxy(it.ret_ptr)) {} template friend class internal_view_iterator; template using valid_args_t = std::enable_if_t< aka::conjunction, std::is_enum>...>::value and dim == sizeof...(Args), int>; public: template = 0> internal_view_iterator(scalar_pointer data, Ns... ns) : dims({Int(ns)...}), _offset(detail::product_all(std::forward(ns)...)), initial(data), ret_ptr(data), proxy(data, ns...), const_proxy(data, ns...) {} // Static 1x1 matrix cannot be differenciated from vector in eigen template , std::enable_if_t * = nullptr> constexpr internal_view_iterator(scalar_pointer data, Idx rows) : dims({rows, 1}), _offset(rows), initial(data), ret_ptr(data), proxy(data, rows, 1), const_proxy(data, rows, 1) { AKANTU_DEBUG_ASSERT(rows == 1, "1x1 Matrix"); } // Static matrix again hard to distinguish from vectors template , std::enable_if_t<(RD::RowsAtCompileTime > 1) and RD::ColsAtCompileTime == 1> * = nullptr> constexpr internal_view_iterator(scalar_pointer data, Idx rows, Idx cols) : dims({rows}), _offset(rows), initial(data), ret_ptr(data), proxy(data, rows, 1), const_proxy(data, rows, 1) { AKANTU_DEBUG_ASSERT(cols == 1, "nx1 Matrix"); } template * = nullptr> internal_view_iterator() : proxy(nullptr, 0), const_proxy(nullptr, 0) {} template * = nullptr> internal_view_iterator() : proxy(nullptr, 0, 0), const_proxy(nullptr, 0, 0) {} internal_view_iterator(const internal_view_iterator & it) : proxy(get_new_proxy(it.ret_ptr)), const_proxy(get_new_const_proxy(it.ret_ptr)) { if (this != &it) { this->dims = it.dims; this->_offset = it._offset; this->initial = it.initial; this->ret_ptr = it.ret_ptr; } } internal_view_iterator(internal_view_iterator && it) = default; virtual ~internal_view_iterator() = default; template ().data()), decltype(std::declval().data())>::value> * = nullptr> inline internal_view_iterator & operator=(const internal_view_iterator & it) { this->_offset = it._offset; this->initial = it.initial; this->ret_ptr = it.ret_ptr; reset_proxy(this->proxy); reset_proxy(this->const_proxy); return *this; } inline internal_view_iterator & operator=(const internal_view_iterator & it) { if (this != &it) { this->_offset = it._offset; this->initial = it.initial; this->ret_ptr = it.ret_ptr; reset_proxy(this->proxy); reset_proxy(this->const_proxy); } return *this; } public: UInt getCurrentIndex() { return (this->ret_ptr - this->initial) / this->_offset; } inline reference operator*() { reset_proxy(proxy); return proxy; } inline const_reference operator*() const { reset_proxy(const_proxy); return proxy; } inline pointer operator->() { reset_proxy(proxy); return &proxy; } inline daughter & operator++() { ret_ptr += _offset; return static_cast(*this); } inline daughter & operator--() { ret_ptr -= _offset; return static_cast(*this); } inline daughter & operator+=(UInt n) { ret_ptr += _offset * n; return static_cast(*this); } inline daughter & operator-=(UInt n) { ret_ptr -= _offset * n; return static_cast(*this); } inline auto operator[](UInt n) { return get_new_proxy(ret_ptr + n * _offset); } inline auto operator[](UInt n) const { return get_new_const_proxy(ret_ptr + n * _offset); } inline bool operator==(const internal_view_iterator & other) const { return this->ret_ptr == other.ret_ptr; } inline bool operator!=(const internal_view_iterator & other) const { return this->ret_ptr != other.ret_ptr; } inline bool operator<(const internal_view_iterator & other) const { return this->ret_ptr < other.ret_ptr; } inline bool operator<=(const internal_view_iterator & other) const { return this->ret_ptr <= other.ret_ptr; } inline bool operator>(const internal_view_iterator & other) const { return this->ret_ptr > other.ret_ptr; } inline bool operator>=(const internal_view_iterator & other) const { return this->ret_ptr >= other.ret_ptr; } inline daughter operator+(difference_type n) { daughter tmp(static_cast(*this)); tmp += n; return tmp; } inline daughter operator-(difference_type n) { daughter tmp(static_cast(*this)); tmp -= n; return tmp; } inline difference_type operator-(const internal_view_iterator & b) { return (this->ret_ptr - b.ret_ptr) / _offset; } inline scalar_pointer data() const { return ret_ptr; } inline difference_type offset() const { return _offset; } inline decltype(auto) getDims() const { return dims; } protected: std::array dims; difference_type _offset{0}; scalar_pointer initial{nullptr}; scalar_pointer ret_ptr{nullptr}; proxy_t proxy; const_proxy_t const_proxy; }; /* ------------------------------------------------------------------------ */ /** * Specialization for scalar types */ template class internal_view_iterator { public: using value_type = R; using pointer = R *; using reference = R &; using const_reference = const R &; using difference_type = std::ptrdiff_t; using iterator_category = std::random_access_iterator_tag; static constexpr int dim_ = 0; protected: using internal_value_type = IR; using internal_pointer = IR *; public: internal_view_iterator(pointer data = nullptr) : ret(data), initial(data){}; internal_view_iterator(const internal_view_iterator & it) = default; internal_view_iterator(internal_view_iterator && it) = default; virtual ~internal_view_iterator() = default; inline internal_view_iterator & operator=(const internal_view_iterator & it) = default; UInt getCurrentIndex() { return (this->ret - this->initial); }; inline reference operator*() { return *ret; } inline const_reference operator*() const { return *ret; } inline pointer operator->() { return ret; }; inline daughter & operator++() { ++ret; return static_cast(*this); } inline daughter & operator--() { --ret; return static_cast(*this); } inline daughter & operator+=(const UInt n) { ret += n; return static_cast(*this); } inline daughter & operator-=(const UInt n) { ret -= n; return static_cast(*this); } inline reference operator[](const UInt n) { return ret[n]; } inline bool operator==(const internal_view_iterator & other) const { return ret == other.ret; } inline bool operator!=(const internal_view_iterator & other) const { return ret != other.ret; } inline bool operator<(const internal_view_iterator & other) const { return ret < other.ret; } inline bool operator<=(const internal_view_iterator & other) const { return ret <= other.ret; } inline bool operator>(const internal_view_iterator & other) const { return ret > other.ret; } inline bool operator>=(const internal_view_iterator & other) const { return ret >= other.ret; } inline daughter operator-(difference_type n) { return daughter(ret - n); } inline daughter operator+(difference_type n) { return daughter(ret + n); } inline difference_type operator-(const internal_view_iterator & b) { return ret - b.ret; } inline pointer data() const { return ret; } inline decltype(auto) getDims() const { return dims; } protected: std::array dims; pointer ret{nullptr}; pointer initial{nullptr}; }; } // namespace detail /* -------------------------------------------------------------------------- */ template class view_iterator; template class const_view_iterator : public detail::internal_view_iterator, R> { public: using parent = detail::internal_view_iterator; using value_type = typename parent::value_type; using pointer = typename parent::pointer; using reference = typename parent::reference; using difference_type = typename parent::difference_type; using iterator_category = typename parent::iterator_category; protected: + template static inline auto convert_helper(const view_iterator & it, std::index_sequence) { return const_view_iterator(it.data(), it.getDims()[I]...); } + template + static inline auto convert_helper(const const_view_iterator & it, + std::index_sequence) { + return const_view_iterator(it.data(), it.getDims()[I]...); + } + template static inline auto convert(const view_iterator & it) { return convert_helper(it, std::make_index_sequence()); } + template + static inline auto convert(const const_view_iterator & it) { + return convert_helper(it, std::make_index_sequence()); + } + public: const_view_iterator() : parent(){}; const_view_iterator(const const_view_iterator & it) = default; const_view_iterator(const_view_iterator && it) = default; template const_view_iterator(P * data, Ns... ns) : parent(data, ns...) {} const_view_iterator & operator=(const const_view_iterator & it) = default; template ::value> * = nullptr> - const_view_iterator(const const_view_iterator & it) : parent(it) {} + const_view_iterator(const const_view_iterator & it) : parent(convert(it)) {} template ::value and std::is_convertible::value> * = nullptr> const_view_iterator & operator=(const const_view_iterator & it) { return dynamic_cast(parent::operator=(it)); } template ::value> * = nullptr> const_view_iterator operator=(const view_iterator & it) { return convert(it); } }; template ::value> struct ConstConverterIteratorHelper { protected: template static inline auto convert_helper(const view_iterator & it, std::index_sequence) { return const_view_iterator(it.data(), it.getDims()[I]...); } public: static inline auto convert(const view_iterator & it) { return convert_helper( it, std::make_index_sequence< std::tuple_size::value>()); } }; template struct ConstConverterIteratorHelper { static inline auto convert(const view_iterator & it) { return const_view_iterator(it.data()); } }; template class view_iterator : public detail::internal_view_iterator> { public: using parent = detail::internal_view_iterator; using value_type = typename parent::value_type; using pointer = typename parent::pointer; using reference = typename parent::reference; using difference_type = typename parent::difference_type; using iterator_category = typename parent::iterator_category; public: view_iterator() : parent(){}; view_iterator(const view_iterator & it) = default; view_iterator(view_iterator && it) = default; template view_iterator(P * data, Ns... ns) : parent(data, ns...) {} view_iterator & operator=(const view_iterator & it) = default; operator const_view_iterator() { return ConstConverterIteratorHelper::convert(*this); } }; namespace { template struct ViewIteratorHelper { using type = TensorProxy; }; template struct ViewIteratorHelper<0, T> { using type = T; }; template struct ViewIteratorHelper<1, T> { using type = Eigen::Map>; }; template struct ViewIteratorHelper<1, const T> { using type = Eigen::Map>; }; template struct ViewIteratorHelper<2, T> { using type = Eigen::Map>; }; template struct ViewIteratorHelper<2, const T> { using type = Eigen::Map>; }; template using ViewIteratorHelper_t = typename ViewIteratorHelper::type; } // namespace // #include // namespace std { // template // struct iterator_traits<::akantu::types::details::vector_iterator> { // protected: // using iterator = ::akantu::types::details::vector_iterator; // public: // using iterator_category = typename iterator::iterator_category; // using value_type = typename iterator::value_type; // using difference_type = typename iterator::difference_type; // using pointer = typename iterator::pointer; // using reference = typename iterator::reference; // }; // } // namespace std } // namespace akantu #endif /* !__AKANTU_AKA_VIEW_ITERATORS_HH__ */ diff --git a/src/fe_engine/element_class.hh b/src/fe_engine/element_class.hh index 6be5ecbb2..10c984aff 100644 --- a/src/fe_engine/element_class.hh +++ b/src/fe_engine/element_class.hh @@ -1,485 +1,491 @@ /** * @file element_class.hh * * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Dec 11 2020 * * @brief Declaration of the ElementClass main class and the * Integration and Interpolation elements * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_CLASS_HH_ #define AKANTU_ELEMENT_CLASS_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /// default element class structure template struct ElementClassProperty { static const GeometricalType geometrical_type{_gt_not_defined}; static const InterpolationType interpolation_type{_itp_not_defined}; static const ElementKind element_kind{_ek_regular}; static const Int spatial_dimension{0}; static const GaussIntegrationType gauss_itegration_type{_git_not_defined}; static const Int polynomial_degree{0}; }; #if !defined(DOXYGEN) /// Macro to generate the element class structures for different element types #define AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(elem_type, geom_type, \ interp_type, elem_kind, sp, \ gauss_int_type, min_int_order) \ template <> struct ElementClassProperty { \ static const GeometricalType geometrical_type{geom_type}; \ static const InterpolationType interpolation_type{interp_type}; \ static const ElementKind element_kind{elem_kind}; \ static const Int spatial_dimension{sp}; \ static const GaussIntegrationType gauss_integration_type{gauss_int_type}; \ static const Int polynomial_degree{min_int_order}; \ } #else #define AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(elem_type, geom_type, \ interp_type, elem_kind, sp, \ gauss_int_type, min_int_order) #endif /* -------------------------------------------------------------------------- */ /* Geometry */ /* -------------------------------------------------------------------------- */ /// Default GeometricalShape structure template struct GeometricalShape { static const GeometricalShapeType shape{_gst_point}; }; /// Templated GeometricalShape with function contains template struct GeometricalShapeContains { /// Check if the point (vector in 2 and 3D) at natural coordinate coord template static inline bool contains(const Eigen::MatrixBase & coord); }; #if !defined(DOXYGEN) /// Macro to generate the GeometricalShape structures for different geometrical /// types #define AKANTU_DEFINE_SHAPE(geom_type, geom_shape) \ template <> struct GeometricalShape { \ static const GeometricalShapeType shape{geom_shape}; \ } AKANTU_DEFINE_SHAPE(_gt_hexahedron_20, _gst_square); AKANTU_DEFINE_SHAPE(_gt_hexahedron_8, _gst_square); AKANTU_DEFINE_SHAPE(_gt_pentahedron_15, _gst_prism); AKANTU_DEFINE_SHAPE(_gt_pentahedron_6, _gst_prism); AKANTU_DEFINE_SHAPE(_gt_point, _gst_point); AKANTU_DEFINE_SHAPE(_gt_quadrangle_4, _gst_square); AKANTU_DEFINE_SHAPE(_gt_quadrangle_8, _gst_square); AKANTU_DEFINE_SHAPE(_gt_segment_2, _gst_square); AKANTU_DEFINE_SHAPE(_gt_segment_3, _gst_square); AKANTU_DEFINE_SHAPE(_gt_tetrahedron_10, _gst_triangle); AKANTU_DEFINE_SHAPE(_gt_tetrahedron_4, _gst_triangle); AKANTU_DEFINE_SHAPE(_gt_triangle_3, _gst_triangle); AKANTU_DEFINE_SHAPE(_gt_triangle_6, _gst_triangle); #endif /* -------------------------------------------------------------------------- */ template struct GeometricalElementProperty {}; template struct ElementClassExtraGeometryProperties {}; /* -------------------------------------------------------------------------- */ /// Templated GeometricalElement with function getInradius template ::shape> class GeometricalElement { using geometrical_property = GeometricalElementProperty; public: /// compute the in-radius: \todo should be renamed for characteristic length template static inline Real getInradius(const Eigen::MatrixBase & /*X*/) { AKANTU_TO_IMPLEMENT(); } /// compute the normal to the element template static inline void getNormal(const Eigen::MatrixBase & /*X*/, Eigen::MatrixBase & /*n*/) { AKANTU_TO_IMPLEMENT(); } /// true if the natural coordinates are in the element template static inline bool contains(const Eigen::MatrixBase & coord); public: static constexpr auto getSpatialDimension() { return geometrical_property::spatial_dimension; } static constexpr auto getNbNodesPerElement() { return geometrical_property::nb_nodes_per_element; } static inline constexpr auto getNbFacetTypes() { return geometrical_property::nb_facet_types; }; static inline constexpr Int getNbFacetsPerElement(Idx t); static inline constexpr Int getNbFacetsPerElement(); static inline constexpr decltype(auto) getFacetLocalConnectivityPerElement(Idx t = 0); template ::value, std::enable_if_t<(t < size)> * = nullptr> static inline constexpr decltype(auto) getFacetLocalConnectivityPerElement(); template ::value, std::enable_if_t * = nullptr> static inline constexpr decltype(auto) getFacetLocalConnectivityPerElement(); }; /* -------------------------------------------------------------------------- */ /* Interpolation */ /* -------------------------------------------------------------------------- */ /// default InterpolationProperty structure template struct InterpolationProperty {}; #if !defined(DOXYGEN) /// Macro to generate the InterpolationProperty structures for different /// interpolation types #define AKANTU_DEFINE_INTERPOLATION_TYPE_PROPERTY(itp_type, itp_kind, \ nb_nodes, ndim) \ template <> struct InterpolationProperty { \ static constexpr InterpolationKind kind{itp_kind}; \ static constexpr Int nb_nodes_per_element{nb_nodes}; \ static constexpr Int natural_space_dimension{ndim}; \ } #else #define AKANTU_DEFINE_INTERPOLATION_TYPE_PROPERTY(itp_type, itp_kind, \ nb_nodes, ndim) #endif /* -------------------------------------------------------------------------- */ /// Generic (templated by the enum InterpolationType which specifies the order /// and the dimension of the interpolation) class handling the elemental /// interpolation template ::kind> class InterpolationElement { public: using interpolation_property = InterpolationProperty; /// compute the shape values for a given set of points in natural coordinates template ::value> * = nullptr> static inline void computeShapes(const Eigen::MatrixBase & natural_coord, const Eigen::MatrixBase & shapes); /// compute the shape values for a given point in natural coordinates template ::value> * = nullptr> static inline void computeShapes(const Eigen::MatrixBase &, Eigen::MatrixBase &) { AKANTU_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ template static inline void computeDNDS(const Eigen::MatrixBase & Xs, Tensor3Base & dnds); /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template static inline void computeDNDS(const Eigen::MatrixBase & /*Xs*/, Eigen::MatrixBase & /*dNdS*/) { AKANTU_TO_IMPLEMENT(); } /** * compute @f$ @f$ **/ static inline void computeD2NDS2(const Matrix & natural_coord, Tensor3 & d2nds2); /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the * second variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template static inline void computeD2NDS2(const vector_type &, matrix_type &) { AKANTU_TO_IMPLEMENT(); } /// compute jacobian (or integration variable change factor) for a given point /// in the case of spatial_dimension != natural_space_dimension template static inline Real computeSpecialJacobian(const Eigen::MatrixBase &) { AKANTU_TO_IMPLEMENT(); } /// interpolate a field given (arbitrary) natural coordinates template static inline decltype(auto) interpolateOnNaturalCoordinates( const Eigen::MatrixBase & natural_coords, const Eigen::MatrixBase & nodal_values) { using interpolation = InterpolationProperty; Eigen::Matrix shapes; computeShapes(natural_coords, shapes); return interpolate(nodal_values, shapes); } /// interpolate a field given the shape functions on the interpolation point template static inline auto interpolate(const Eigen::MatrixBase & nodal_values, const Eigen::MatrixBase & shapes); /// interpolate a field given the shape functions on the interpolations points template static inline void interpolate(const Eigen::MatrixBase & nodal_values, const Eigen::MatrixBase & Ns, const Eigen::MatrixBase & interpolated); /// compute the gradient of a given field on the given natural coordinates template static inline void gradientOnNaturalCoordinates(const Eigen::MatrixBase & natural_coords, const Eigen::MatrixBase & f, const Eigen::MatrixBase & dfds); public: static constexpr auto getShapeSize() { return InterpolationProperty::nb_nodes_per_element; } static constexpr auto getShapeDerivativesSize() { return (InterpolationProperty::nb_nodes_per_element * InterpolationProperty::natural_space_dimension); } static constexpr auto getNaturalSpaceDimension() { return InterpolationProperty::natural_space_dimension; } static constexpr auto getNbNodesPerInterpolationElement() { return InterpolationProperty::nb_nodes_per_element; } }; /* -------------------------------------------------------------------------- */ /* Integration */ /* -------------------------------------------------------------------------- */ template struct GaussIntegrationTypeData { /// quadrature points in natural coordinates static Real quad_positions[]; /// weights for the Gauss integration static Real quad_weights[]; }; template ::polynomial_degree> class GaussIntegrationElement { public: static constexpr Int getNbQuadraturePoints(); static constexpr decltype(auto) getQuadraturePoints(); static constexpr decltype(auto) getWeights(); }; /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ template ::element_kind> class ElementClass : public GeometricalElement< ElementClassProperty::geometrical_type>, public InterpolationElement< ElementClassProperty::interpolation_type> { protected: using geometrical_element = GeometricalElement::geometrical_type>; using interpolation_element = InterpolationElement< ElementClassProperty::interpolation_type>; using element_property = ElementClassProperty; using interpolation_property = typename interpolation_element::interpolation_property; public: /** * compute @f$ J = \frac{\partial x_j}{\partial s_i} @f$ the variation of real * coordinates along with variation of natural coordinates on a given point in * natural coordinates */ template static inline decltype(auto) computeJMat(const Eigen::MatrixBase & dnds, const Eigen::MatrixBase & node_coords); /** * compute the Jacobian matrix by computing the variation of real coordinates * along with variation of natural coordinates on a given set of points in * natural coordinates */ template static inline void computeJMat(const Tensor3Base & dnds, const Eigen::MatrixBase & node_coords, Tensor3Base & J); /// compute the jacobians of a serie of natural coordinates template static inline void computeJacobian(const Eigen::MatrixBase & natural_coords, const Eigen::MatrixBase & node_coords, Eigen::MatrixBase & jacobians); /// compute jacobian (or integration variable change factor) for a set of /// points template static inline void computeJacobian(const Tensor3Base & J, Eigen::MatrixBase & jacobians); /// compute jacobian (or integration variable change factor) for a given point template static inline Real computeJacobian(const Eigen::MatrixBase & J); /// compute shape derivatives (input is dxds) for a set of points static inline void computeShapeDerivatives(const Tensor3Base & J, const Tensor3Base & dnds, Tensor3Base & shape_deriv); /// compute shape derivatives (input is dxds) for a given point template static inline void computeShapeDerivatives(const Eigen::MatrixBase & J, const Eigen::MatrixBase & dnds, Eigen::MatrixBase & shape_deriv); /// compute the normal of a surface defined by the function f template static inline void computeNormalsOnNaturalCoordinates(const Eigen::MatrixBase & coord, const Eigen::MatrixBase & f, Eigen::MatrixBase & normals); /// get natural coordinates from real coordinates template static inline void inverseMap(const Eigen::MatrixBase & real_coords, const Eigen::MatrixBase & node_coords, Eigen::MatrixBase & natural_coords, Real tolerance = 1e-10); /// get natural coordinates from real coordinates template * = nullptr> static inline void inverseMap(const Eigen::MatrixBase & real_coords, const Eigen::MatrixBase & node_coords, const Eigen::MatrixBase & natural_coords_, UInt max_iterations = 100, Real tolerance = 1e-10); public: static constexpr auto getKind() { return element_kind; } static constexpr auto getSpatialDimension() { return ElementClassProperty::spatial_dimension; } using element_class_extra_geom_property = ElementClassExtraGeometryProperties; static constexpr decltype(auto) getP1ElementType() { return element_class_extra_geom_property::p1_type; } static constexpr decltype(auto) getFacetType(UInt t = 0) { return element_class_extra_geom_property::facet_type[t]; } static constexpr decltype(auto) getFacetTypes(); }; /* -------------------------------------------------------------------------- */ } // namespace akantu /* -------------------------------------------------------------------------- */ #include "geometrical_element_property.hh" #include "interpolation_element_tmpl.hh" /* -------------------------------------------------------------------------- */ #include "element_class_tmpl.hh" /* -------------------------------------------------------------------------- */ +namespace akantu { +AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_not_defined, _gt_not_defined, + _itp_not_defined, _ek_not_defined, 0, + _git_not_defined, 0); +} // namespace akantu + #include "element_class_hexahedron_8_inline_impl.cc" #include "element_class_pentahedron_6_inline_impl.cc" /* keep order */ #include "element_class_hexahedron_20_inline_impl.cc" #include "element_class_pentahedron_15_inline_impl.cc" #include "element_class_point_1_inline_impl.cc" #include "element_class_quadrangle_4_inline_impl.cc" #include "element_class_quadrangle_8_inline_impl.cc" #include "element_class_segment_2_inline_impl.cc" #include "element_class_segment_3_inline_impl.cc" #include "element_class_tetrahedron_10_inline_impl.cc" #include "element_class_tetrahedron_4_inline_impl.cc" #include "element_class_triangle_3_inline_impl.cc" #include "element_class_triangle_6_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "element_class_structural.hh" #endif #if defined(AKANTU_COHESIVE_ELEMENT) #include "cohesive_element.hh" #endif #if defined(AKANTU_IGFEM) #include "element_class_igfem.hh" #endif #endif /* AKANTU_ELEMENT_CLASS_HH_ */ diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh index b6c8bac0a..d4a28b5b0 100644 --- a/src/fe_engine/fe_engine_template.hh +++ b/src/fe_engine/fe_engine_template.hh @@ -1,450 +1,466 @@ /** * @file fe_engine_template.hh * * @author Guillaume Anciaux * @author Sébastien Hartmann * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri May 14 2021 * * @brief templated class that calls integration and shape objects * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fe_engine.hh" //#include "integrator.hh" //#include "shape_functions.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_FE_ENGINE_TEMPLATE_HH_ #define AKANTU_FE_ENGINE_TEMPLATE_HH_ namespace akantu { class Integrator; class ShapeFunctions; } // namespace akantu namespace akantu { class DOFManager; namespace fe_engine { namespace details { template struct AssembleLumpedTemplateHelper; template struct AssembleFieldMatrixHelper; } // namespace details } // namespace fe_engine template struct AssembleFieldMatrixStructHelper; struct DefaultIntegrationOrderFunctor { template static inline constexpr int getOrder() { return ElementClassProperty::polynomial_degree; } }; /* -------------------------------------------------------------------------- */ template