diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index d2695c476..2fe7ca644 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,667 +1,668 @@ /** * @file aka_common.hh * * @author Nicolas Richart * * @date Mon Jun 14 19:12:20 2010 * * @brief common type descriptions for akantu * * @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 . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ }; #if defined(WIN32) # define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ typedef double Real; typedef unsigned int UInt; typedef unsigned long long UInt64; typedef signed int Int; typedef std::string ID; static const Real UINT_INIT_VALUE = 0; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = 0; #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ typedef UInt MemoryID; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ typedef std::string Surface; typedef std::pair SurfacePair; typedef std::list< SurfacePair > SurfacePairList; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; /// @boost sequence of element to loop on in global tasks #define AKANTU_ek_regular_ELEMENT_TYPE \ (_point_1) \ (_segment_2) \ (_segment_3) \ (_triangle_3) \ (_triangle_6) \ (_quadrangle_4) \ (_quadrangle_8) \ (_tetrahedron_4) \ (_tetrahedron_10) \ (_pentahedron_6) \ (_hexahedron_8) #if defined(AKANTU_STRUCTURAL_MECHANICS) #define AKANTU_ek_structural_ELEMENT_TYPE \ (_bernoulli_beam_2) \ (_bernoulli_beam_3) \ (_kirchhoff_shell) #else #define AKANTU_ek_structural_ELEMENT_TYPE #endif #if defined(AKANTU_COHESIVE_ELEMENT) # define AKANTU_ek_cohesive_ELEMENT_TYPE \ (_cohesive_2d_4) \ (_cohesive_2d_6) \ (_cohesive_1d_2) \ (_cohesive_3d_6) \ (_cohesive_3d_12) #else # define AKANTU_ek_cohesive_ELEMENT_TYPE #endif #if defined(AKANTU_IGFEM) #define AKANTU_ek_igfem_ELEMENT_TYPE \ (_igfem_triangle_3) #else #define AKANTU_ek_igfem_ELEMENT_TYPE #endif #define AKANTU_ALL_ELEMENT_TYPE \ AKANTU_ek_regular_ELEMENT_TYPE \ AKANTU_ek_cohesive_ELEMENT_TYPE \ AKANTU_ek_structural_ELEMENT_TYPE \ AKANTU_ek_igfem_ELEMENT_TYPE #define AKANTU_NOT_STRUCTURAL_ELEMENT_TYPE \ AKANTU_ek_regular_ELEMENT_TYPE \ AKANTU_ek_cohesive_ELEMENT_TYPE \ AKANTU_ek_igfem_ELEMENT_TYPE /// @enum ElementType type of elements enum ElementType { _not_defined, _point_1, _segment_2, ///< first order segment _segment_3, ///< second order segment _triangle_3, ///< first order triangle _triangle_6, ///< second order triangle _tetrahedron_4, ///< first order tetrahedron _tetrahedron_10, ///< second order tetrahedron _quadrangle_4, ///< first order quadrangle _quadrangle_8, ///< second order quadrangle _hexahedron_8, ///< first order hexahedron _pentahedron_6, ///< first order pentahedron #if defined (AKANTU_STRUCTURAL_MECHANICS) _bernoulli_beam_2, ///< Bernoulli beam 2D _bernoulli_beam_3, ///< Bernoulli beam 3D _kirchhoff_shell, ///< Kirchhoff shell #endif #if defined(AKANTU_COHESIVE_ELEMENT) _cohesive_2d_4, ///< first order 2D cohesive _cohesive_2d_6, ///< second order 2D cohesive _cohesive_1d_2, ///< first order 1D cohesive _cohesive_3d_6, ///< first order 3D cohesive _cohesive_3d_12, ///< second order 3D cohesive #endif #if defined(AKANTU_IGFEM) _igfem_triangle_3, ///< first order triangle for IGFEM #endif _max_element_type }; /// @enum GeometricalType type of element potentially contained in a Mesh enum GeometricalType { _gt_point, ///< point @remark only for some algorithm to be generic like mesh partitioning */ _gt_segment_2, ///< 2 nodes segment _gt_segment_3, ///< 3 nodes segment _gt_triangle_3, ///< 3 nodes triangle _gt_triangle_6, ///< 6 nodes triangle _gt_quadrangle_4, ///< 4 nodes quadrangle _gt_quadrangle_8, ///< 8 nodes quadrangle _gt_tetrahedron_4, ///< 4 nodes tetrahedron _gt_tetrahedron_10, ///< 10 nodes tetrahedron _gt_hexahedron_8, ///< 8 nodes hexahedron _gt_pentahedron_6, ///< 6 nodes pentahedron #if defined(AKANTU_COHESIVE_ELEMENT) _gt_cohesive_2d_4, ///< 4 nodes 2D cohesive _gt_cohesive_2d_6, ///< 6 nodes 2D cohesive _gt_cohesive_1d_2, ///< 2 nodes 1D cohesive _gt_cohesive_3d_6, ///< 6 nodes 3D cohesive _gt_cohesive_3d_12, ///< 12 nodes 3D cohesive #endif _gt_not_defined }; /// @enum InterpolationType type of elements enum InterpolationType { _itp_lagrange_point_1, ///< zeroth (!) order lagrangian point (for compatibility purposes) _itp_lagrange_segment_2, ///< first order lagrangian segment _itp_lagrange_segment_3, ///< second order lagrangian segment _itp_lagrange_triangle_3, ///< first order lagrangian triangle _itp_lagrange_triangle_6, ///< second order lagrangian triangle _itp_lagrange_quadrangle_4, ///< first order lagrangian quadrangle _itp_serendip_quadrangle_8, /**< second order serendipian quadrangle @remark used insted of the 9 node lagrangian element */ _itp_lagrange_tetrahedron_4, ///< first order lagrangian tetrahedron _itp_lagrange_tetrahedron_10, ///< second order lagrangian tetrahedron _itp_lagrange_hexahedron_8, ///< first order lagrangian hexahedron _itp_lagrange_pentahedron_6, ///< first order lagrangian pentahedron #if defined(AKANTU_STRUCTURAL_MECHANICS) _itp_bernoulli_beam, ///< Bernoulli beam _itp_kirchhoff_shell, ///< Kirchhoff shell #endif _itp_not_defined }; //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type); #define AKANTU_REGULAR_KIND (_ek_regular) #ifdef AKANTU_COHESIVE_ELEMENT # define AKANTU_COHESIVE_KIND (_ek_cohesive) #else # define AKANTU_COHESIVE_KIND #endif #ifdef AKANTU_STRUCTURAL_MECHANICS # define AKANTU_STRUCTURAL_KIND (_ek_structural) #else # define AKANTU_STRUCTURAL_KIND #endif #ifdef AKANTU_IGFEM # define AKANTU_IGFEM_KIND (_ek_igfem) #else # define AKANTU_IGFEM_KIND #endif #define AKANTU_ELEMENT_KIND \ AKANTU_REGULAR_KIND \ AKANTU_COHESIVE_KIND \ AKANTU_STRUCTURAL_KIND \ AKANTU_IGFEM_KIND enum ElementKind { BOOST_PP_SEQ_ENUM(AKANTU_ELEMENT_KIND), _ek_not_defined }; /// 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 AnalysisMethod type of solving method used to solve the equation of motion enum AnalysisMethod { _static, _implicit_dynamic, _explicit_lumped_mass, _explicit_lumped_capacity, _explicit_consistent_mass }; /// enum SolveConvergenceMethod different resolution algorithms enum SolveConvergenceMethod { _scm_newton_raphson_tangent, ///< Newton-Raphson with tangent matrix _scm_newton_raphson_tangent_modified, ///< Newton-Raphson with constant tangent matrix _scm_newton_raphson_tangent_not_computed ///< Newton-Raphson with constant tangent matrix (user has to assemble it) }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence - _scc_increment ///< Use increment to test the convergence + _scc_increment, ///< Use increment 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 }; /// myfunction(double * position, double * stress/force, double * normal, unsigned int material_id) typedef void (*BoundaryFunction)(double *,double *, double*, unsigned int); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_traction, _bft_traction_local }; /// @enum SparseMatrixType type of sparse matrix used enum SparseMatrixType { _unsymmetric, _symmetric }; /* -------------------------------------------------------------------------- */ /* Contact */ /* -------------------------------------------------------------------------- */ typedef ID ContactID; typedef ID ContactSearchID; typedef ID ContactNeighborStructureID; enum ContactType { _ct_not_defined = 0, _ct_2d_expli = 1, _ct_3d_expli = 2, _ct_rigid = 3 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1, _cnst_2d_grid = 2 }; /* -------------------------------------------------------------------------- */ /* Friction */ /* -------------------------------------------------------------------------- */ typedef ID FrictionID; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- 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 //--- CohesiveElementInserter tags --- _gst_ce_inserter, //< synchronization of global nodes id of newly inserted cohesive elements //--- 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 --- /// synchronization of the nodal level set value phi _gst_htm_phi, /// synchronization of the element gradient phi _gst_htm_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 //--- General tags --- _gst_test, //< Test tag _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 }; /// 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 { typedef GhostType type; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; typedef safe_enum ghost_type_t; /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type); /// @enum SynchronizerOperation reduce operation that the synchronizer can perform enum SynchronizerOperation { _so_sum, _so_min, _so_max, _so_null }; /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ // int 2 type construct template struct Int2Type { enum { result = d }; }; // type 2 type construct template class Type2Type { typedef T OriginalType; }; /* -------------------------------------------------------------------------- */ template struct is_scalar { enum{ value = false }; }; #define AKANTU_SPECIFY_IS_SCALAR(type) \ template<> \ struct is_scalar { \ enum { value = true }; \ } AKANTU_SPECIFY_IS_SCALAR(Real); AKANTU_SPECIFY_IS_SCALAR(UInt); AKANTU_SPECIFY_IS_SCALAR(Int); AKANTU_SPECIFY_IS_SCALAR(bool); template < typename T1, typename T2 > struct is_same { enum { value = false }; // is_same represents a bool. }; template < typename T > struct is_same { enum { value = true }; }; /* -------------------------------------------------------------------------- */ #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 (); /* -------------------------------------------------------------------------- */ /* * 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); /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "aka_fwd.hh" __BEGIN_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(); __END_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_DEBUG_ERROR("Type (" << var << ") not handled by this function"); \ } \ } \ } while(0) #define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1) \ AKANTU_BOOST_LIST_SWITCH(macro1, list1, type) #define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_ek_regular_ELEMENT_TYPE) #define AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_ek_cohesive_ELEMENT_TYPE) #define AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_ek_structural_ELEMENT_TYPE) #define AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_ek_igfem_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_BOOST_REGULAR_ELEMENT_LIST(macro) \ AKANTU_BOOST_APPLY_ON_LIST(macro, \ AKANTU_ek_regular_ELEMENT_TYPE) #define AKANTU_BOOST_STRUCTURAL_ELEMENT_LIST(macro) \ AKANTU_BOOST_APPLY_ON_LIST(macro, \ AKANTU_ek_structural_ELEMENT_TYPE) #define AKANTU_BOOST_COHESIVE_ELEMENT_LIST(macro) \ AKANTU_BOOST_APPLY_ON_LIST(macro, \ AKANTU_ek_cohesive_ELEMENT_TYPE) #define AKANTU_BOOST_IGFEM_ELEMENT_LIST(macro) \ AKANTU_BOOST_APPLY_ON_LIST(macro, \ AKANTU_ek_igfem_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) /// define kept for compatibility reasons (they are most probably not needed /// anymore) \todo check if they can be removed #define AKANTU_REGULAR_ELEMENT_TYPE AKANTU_ek_regular_ELEMENT_TYPE #define AKANTU_COHESIVE_ELEMENT_TYPE AKANTU_ek_cohesive_ELEMENT_TYPE #define AKANTU_STRUCTURAL_ELEMENT_TYPE AKANTU_ek_structural_ELEMENT_TYPE #define AKANTU_IGFEM_ELEMENT_TYPE AKANTU_ek_igfem_ELEMENT_TYPE #include "aka_common_inline_impl.cc" #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/common/aka_common_inline_impl.cc b/src/common/aka_common_inline_impl.cc index d0fc9781e..b113f8b1d 100644 --- a/src/common/aka_common_inline_impl.cc +++ b/src/common/aka_common_inline_impl.cc @@ -1,232 +1,232 @@ /** * @file aka_common_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Thu Dec 01 12:54:29 2011 * * @brief inline implementations of common akantu type descriptions * * @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 . * * @section DESCRIPTION * * All common things to be included in the projects files * */ #include #include #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type) { #define STRINGIFY(type) \ stream << BOOST_PP_STRINGIZE(type) switch(type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, \ STRINGIFY, \ AKANTU_ALL_ELEMENT_TYPE) case _not_defined: stream << "_not_defined"; break; case _max_element_type: stream << "_max_element_type"; break; } #undef STRINGIFY return stream; } /* -------------------------------------------------------------------------- */ //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementKind kind ) { #define STRINGIFY(kind) \ stream << BOOST_PP_STRINGIZE(kind) AKANTU_BOOST_ALL_KIND_SWITCH(STRINGIFY); #undef STRINGIFY return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for InterpolationType inline std::ostream & operator <<(std::ostream & stream, InterpolationType type) { switch(type) { case _itp_lagrange_point_1 : stream << "_itp_lagrange_point_1" ; break; case _itp_lagrange_segment_2 : stream << "_itp_lagrange_segment_2" ; break; case _itp_lagrange_segment_3 : stream << "_itp_lagrange_segment_3" ; break; case _itp_lagrange_triangle_3 : stream << "_itp_lagrange_triangle_3" ; break; case _itp_lagrange_triangle_6 : stream << "_itp_lagrange_triangle_6" ; break; case _itp_lagrange_quadrangle_4 : stream << "_itp_lagrange_quadrangle_4" ; break; case _itp_serendip_quadrangle_8 : stream << "_itp_serendip_quadrangle_8" ; break; case _itp_lagrange_tetrahedron_4 : stream << "_itp_lagrange_tetrahedron_4" ; break; case _itp_lagrange_tetrahedron_10 : stream << "_itp_lagrange_tetrahedron_10"; break; case _itp_lagrange_hexahedron_8 : stream << "_itp_lagrange_hexahedron_8" ; break; case _itp_lagrange_pentahedron_6 : stream << "_itp_lagrange_pentahedron_6" ; break; #if defined(AKANTU_STRUCTURAL_MECHANICS) case _itp_bernoulli_beam : stream << "_itp_bernoulli_beam" ; break; case _itp_kirchhoff_shell : stream << "_itp_kirchhoff_shell" ; break; #endif case _itp_not_defined : stream << "_itp_not_defined" ; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type) { switch(type) { case _not_ghost : stream << "not_ghost"; break; case _ghost : stream << "ghost" ; break; case _casper : stream << "Casper the friendly ghost"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SynchronizationTag inline std::ostream & operator <<(std::ostream & stream, SynchronizationTag type) { switch(type) { case _gst_smm_mass : stream << "_gst_smm_mass" ; break; case _gst_smm_for_gradu : stream << "_gst_smm_for_gradu" ; break; case _gst_smm_boundary : stream << "_gst_smm_boundary" ; break; case _gst_smm_uv : stream << "_gst_smm_uv" ; break; case _gst_smm_res : stream << "_gst_smm_res" ; break; case _gst_smm_init_mat : stream << "_gst_smm_init_mat" ; break; case _gst_smm_stress : stream << "_gst_smm_stress" ; break; case _gst_smmc_facets : stream << "_gst_smmc_facets" ; break; case _gst_smmc_facets_conn : stream << "_gst_smmc_facets_conn" ; break; case _gst_smmc_facets_stress : stream << "_gst_smmc_facets_stress" ; break; case _gst_smmc_damage : stream << "_gst_smmc_damage" ; break; case _gst_ce_inserter : stream << "_gst_ce_inserter" ; break; case _gst_gm_clusters : stream << "_gst_gm_clusters" ; break; case _gst_htm_capacity : stream << "_gst_htm_capacity" ; break; case _gst_htm_temperature : stream << "_gst_htm_temperature" ; break; case _gst_htm_gradient_temperature : stream << "_gst_htm_gradient_temperature"; break; case _gst_htm_phi : stream << "_gst_htm_phi" ; break; case _gst_htm_gradient_phi : stream << "_gst_htm_gradient_phi" ; break; case _gst_mnl_for_average : stream << "_gst_mnl_for_average" ; break; case _gst_mnl_weight : stream << "_gst_mnl_weight" ; break; case _gst_test : stream << "_gst_test" ; break; case _gst_material_id : stream << "_gst_material_id" ; break; case _gst_for_dump : stream << "_gst_for_dump" ; break; case _gst_cf_nodal : stream << "_gst_cf_nodal" ; break; case _gst_cf_incr : stream << "_gst_cf_incr" ; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SolveConvergenceCriteria inline std::ostream & operator <<(std::ostream & stream, SolveConvergenceCriteria criteria) { - switch(criteria) - { - case _scc_residual : stream << "_scc_residual" ; break; - case _scc_increment: stream << "_scc_increment"; break; - } + switch(criteria) { + case _scc_residual : stream << "_scc_residual" ; break; + case _scc_increment: stream << "_scc_increment"; break; + case _scc_residual_mass_wgh: stream << "_scc_residual_mass_wgh"; break; + } return stream; } /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str) { std::string lstr = str; std::transform(lstr.begin(), lstr.end(), lstr.begin(), (int(*)(int))tolower); return lstr; } /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim) { std::string trimed = to_trim; //left trim trimed.erase(trimed.begin(), std::find_if(trimed.begin(), trimed.end(), std::not1(std::ptr_fun(isspace)))); // right trim trimed.erase(std::find_if(trimed.rbegin(), trimed.rend(), std::not1(std::ptr_fun(isspace))).base(), trimed.end()); return trimed; } __END_AKANTU__ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template std::string printMemorySize(UInt size) { Real real_size = size * sizeof(T); UInt mult = (std::log(real_size) / std::log(2)) / 10; std::stringstream sstr; real_size /= Real(1 << (10 * mult)); sstr << std::setprecision(2) << std::fixed << real_size; std::string size_prefix; switch(mult) { case 0: sstr << ""; break; case 1: sstr << "Ki"; break; case 2: sstr << "Mi"; break; case 3: sstr << "Gi"; break; // I started on this type of machines // (32bit computers) (Nicolas) case 4: sstr << "Ti"; break; case 5: sstr << "Pi"; break; case 6: sstr << "Ei"; break; // theoritical limit of RAM of the current // computers in 2014 (64bit computers) (Nicolas) case 7: sstr << "Zi"; break; case 8: sstr << "Yi"; break; default: AKANTU_DEBUG_ERROR("The programmer in 2014 didn't thought so far (even wikipedia does not go further)." << " You have at least 1024 times more than a yobibit of RAM!!!" << " Just add the prefix corresponding in this switch case."); } sstr << "Byte"; return sstr.str(); } __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 39f1adbfd..f315bafea 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1989 +1,2032 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the SolidMechanicsModel class * * @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_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "integration_scheme_2nd_order.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" # include "dumper_iohelper_tmpl.hh" # include "dumper_iohelper_tmpl_homogenizing_field.hh" # include "dumper_iohelper_tmpl_material_internal_field.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ using std::cout; using std::endl; const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), Dumpable(), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), element_index_by_material("element index by material", id), material_selector(new DefaultMaterialSelector(element_index_by_material)), is_default_material_selector(true), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL), are_materials_instantiated(false) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEEngineObject("SolidMechanicsFEEngine", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->blocked_dofs = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->previous_displacement = NULL; materials.clear(); mesh.registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; delete solver; delete mass_matrix; delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; delete jacobian_matrix; delete synch_parallel; if(is_default_material_selector) { delete material_selector; material_selector = NULL; } AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); method = smm_options.analysis_method; // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } // initialize the materials if(this->parser->getLastParsedFile() != "") { instantiateMaterials(); } if(!smm_options.no_init_materials) { initMaterials(); } if(increment_flag) initBC(*this, *displacement, *increment, *force); else initBC(*this, *displacement, *force); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; synch_parallel = &createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnControlPoints(_not_ghost); fem_boundary.computeNormalsOnControlPoints(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) { AKANTU_DEBUG_IN(); //in case of switch from implicit to explicit if(!this->isExplicit()) method = analysis_method; if (integrator) delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); std::stringstream sstr; sstr << id << ":increment_acceleration"; increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":previous_displacement"; previous_displacement = &(alloc (sstr_disp_t.str(), nb_nodes, spatial_dimension, 0.)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); element_index_by_material.alloc(nb_element, 2, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); synch_registry->registerSynchronizer(*synch, _gst_for_dump); changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->storage(); Real * position_val = mesh.getNodes().storage(); Real * displacement_val = displacement->storage(); /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->storage(), force->storage(), nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); // f = f_ext - f_int // f = f_ext if(need_initialize) initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant non local stresses"); synch_registry->waitEndSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant stresses"); // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (isExplicit()) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Update the residual"); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->storage(); Real * accel_val = acceleration->storage(); Real * res_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } blocked_dofs_val++; res_val++; mass_val++; accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv if(velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_lumped_mass) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs, f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & blocked_dofs, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; blocked_dofs_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { if(previous_displacement) increment->copy(*previous_displacement); else increment->copy(*displacement); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); if(increment_flag) { Real * inc_val = increment->storage(); Real * dis_val = displacement->storage(); UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent(); for (UInt n = 0; n < nb_degree_of_freedom; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment_acceleration); if(previous_displacement) previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); } #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS if(solver) solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initJacobianMatrix() { #ifdef AKANTU_USE_MUMPS // @todo make it more flexible: this is an ugly patch to treat the case of non // fix profile of the K matrix delete jacobian_matrix; std::stringstream sstr_sti; sstr_sti << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id); std::stringstream sstr_solv; sstr_solv << id << ":solver"; delete solver; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); if(solver) solver->initialize(_solver_no_options); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; if (!increment) setIncrementFlagOn(); initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*blocked_dofs); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } // /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::solveStatic(Array & boundary_normal, Array & EulerAngles) { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_INFO("Solving Ku = f"); // AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, // "You should first initialize the implicit solver and assemble the stiffness matrix"); // UInt nb_nodes = displacement->getSize(); // UInt nb_degree_of_freedom = displacement->getNbComponent(); // // if(method != _static) // jacobian_matrix->copyContent(*stiffness_matrix); // Array * residual_rotated = new Array (nb_nodes, nb_degree_of_freedom, "residual_rotated"); // //stiffness_matrix->saveMatrix("stiffness_original.out"); // jacobian_matrix->applyBoundaryNormal(boundary_normal, EulerAngles, *residual, (*stiffness_matrix).getA(), *residual_rotated); // //jacobian_matrix->saveMatrix("stiffness_rotated_dir.out"); // jacobian_matrix->applyBoundary(*blocked_dofs); // solver->setRHS(*residual_rotated); // delete residual_rotated; // if (!increment) setIncrementFlagOn(); // solver->solve(*increment); // Matrix T(nb_degree_of_freedom, nb_degree_of_freedom); // Matrix small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); // Matrix T_small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); // Real * increment_val = increment->storage(); // Real * displacement_val = displacement->storage(); // bool * blocked_dofs_val = blocked_dofs->storage(); // for (UInt n = 0; n < nb_nodes; ++n) { // bool constrain_ij = false; // for (UInt j = 0; j < nb_degree_of_freedom; j++) { // if (boundary_normal(n, j)) { // constrain_ij = true; // break; // } // } // if (constrain_ij) { // if (nb_degree_of_freedom == 2) { // Real Theta = EulerAngles(n, 0); // T(0, 0) = cos(Theta); // T(0, 1) = -sin(Theta); // T(1, 1) = cos(Theta); // T(1, 0) = sin(Theta); // } else if (nb_degree_of_freedom == 3) { // Real Theta_x = EulerAngles(n, 0); // Real Theta_y = EulerAngles(n, 1); // Real Theta_z = EulerAngles(n, 2); // T(0, 0) = cos(Theta_y) * cos(Theta_z); // T(0, 1) = -cos(Theta_y) * sin(Theta_z); // T(0, 2) = sin(Theta_y); // T(1, 0) = cos(Theta_x) * sin(Theta_z) + cos(Theta_z) * sin(Theta_x) * sin(Theta_y); // T(1, 1) = cos(Theta_x) * cos(Theta_z) - sin(Theta_x) * sin(Theta_y) * sin(Theta_z); // T(1, 2) = -cos(Theta_y) * sin(Theta_x); // T(2, 0) = sin(Theta_x) * sin(Theta_z) - cos(Theta_x) * cos(Theta_z) * sin(Theta_y); // T(2, 1) = cos(Theta_z) * sin(Theta_x) + cos(Theta_x) * sin(Theta_y) * sin(Theta_z); // T(2, 2) = cos(Theta_x) * cos(Theta_y); // } // small_rhs.clear(); // T_small_rhs.clear(); // for (UInt j = 0; j < nb_degree_of_freedom; j++) // if(!(boundary_normal(n,j)) ) // small_rhs(j,j)=increment_val[j]; // T_small_rhs.mul(T,small_rhs); // for (UInt j = 0; j < nb_degree_of_freedom; j++){ // if(!(boundary_normal(n,j))){ // for (UInt k = 0; k < nb_degree_of_freedom; k++) // displacement_val[k]+=T_small_rhs(k,j); // } // } // displacement_val += nb_degree_of_freedom; // blocked_dofs_val += nb_degree_of_freedom; // increment_val += nb_degree_of_freedom; // } else { // for (UInt j = 0; j < nb_degree_of_freedom; j++, ++displacement_val, ++increment_val, ++blocked_dofs_val) { // if (!(*blocked_dofs_val)) { // *displacement_val += *increment_val; // } // } // } // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { if(!velocity_damping_matrix) velocity_damping_matrix = new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } blocked_dofs_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); if (norm[1] < Math::getTolerance()) { error = norm[0]; AKANTU_DEBUG_OUT(); // cout<<"Error 1: "< Math::getTolerance()) error = norm[0] / norm[1]; else error = norm[0]; //In case the total displacement is zero! // cout<<"Error 2: "< bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*blocked_dofs_val)) { norm += *residual_val * *residual_val; } blocked_dofs_val++; residual_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } +/* -------------------------------------------------------------------------- */ +template<> +bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance, + Real & norm) { + AKANTU_DEBUG_IN(); + + + + UInt nb_nodes = residual->getSize(); + + norm = 0; + Real * residual_val = residual->storage(); + Real * mass_val = this->mass->storage(); + bool * blocked_dofs_val = blocked_dofs->storage(); + + for (UInt n = 0; n < nb_nodes; ++n) { + bool is_local_node = mesh.isLocalOrMasterNode(n); + if(is_local_node) { + for (UInt d = 0; d < spatial_dimension; ++d) { + if(!(*blocked_dofs_val)) { + norm += *residual_val * *residual_val/(*mass_val * *mass_val); + } + blocked_dofs_val++; + residual_val++; + mass_val++; + } + } else { + blocked_dofs_val += spatial_dimension; + residual_val += spatial_dimension; + mass_val += spatial_dimension; + } + } + + StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); + + norm = sqrt(norm); + + AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); + + AKANTU_DEBUG_OUT(); + return (norm < tolerance); +} + /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); if(previous_displacement) previous_displacement->copy(*displacement); if(method == _implicit_dynamic) integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); if(method == _implicit_dynamic) { integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); bool * boun_val = blocked_dofs->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){ *incr_val *= (1. - *boun_val); *disp_val += *incr_val; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateIncrement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); Real * prev_disp_val = previous_displacement->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val) *incr_val = (*disp_val - *prev_disp_val); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updatePreviousDisplacement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0.)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::iterator< Vector > eibm = element_index_by_material(*it, ghost_type).begin(2); Array X(0, nb_nodes_per_element*spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::matrix_iterator X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) { elem.element = (*eibm)(1); Real el_h = getFEEngine().getElementInradius(*X_el, *it); Real el_c = mat_val[(*eibm)(0)]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->storage(); Real * mass_val = mass->storage(); for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnQuadraturePoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::vector_iterator vit = vel_on_quad.begin(spatial_dimension); Array::vector_iterator vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[element_index_by_material(type)(index, 0)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = blocked_dofs->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work"){ return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(energy_id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index){ AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } std::vector::iterator mat_it; Vector mat = element_index_by_material(type, _not_ghost).begin(2)[index]; Real energy = materials[mat(0)]->getEnergy(energy_id, type, mat(1)); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if(displacement) displacement->resize(nb_nodes); if(mass ) mass ->resize(nb_nodes); if(velocity ) velocity ->resize(nb_nodes); if(acceleration) acceleration->resize(nb_nodes); if(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(blocked_dofs) blocked_dofs->resize(nb_nodes); if(previous_displacement) previous_displacement->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); if(current_position) current_position->resize(nb_nodes); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } if (method != _explicit_lumped_mass) { delete stiffness_matrix; delete jacobian_matrix; delete solver; SolverOptions solver_options; initImplicit((method == _implicit_dynamic), solver_options); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); /// \todo have rules to choose the correct material UInt mat_id = 0; UInt * mat_id_vect = NULL; try { const NewMaterialElementsEvent & event_mat = dynamic_cast(event); mat_id_vect = event_mat.getMaterialList().storage(); } catch(...) { } for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if(mat_id_vect) mat_id = mat_id_vect[el]; else mat_id = (*material_selector)(elem); Material & mat = *materials[mat_id]; UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type); Vector id(2); id[0] = mat_id; id[1] = mat_index; if(!element_index_by_material.exists(elem.type, elem.ghost_type)) element_index_by_material.alloc(0, 2, elem.type, elem.ghost_type); element_index_by_material(elem.type, elem.ghost_type).push_back(id); } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } if(method != _explicit_lumped_mass) AKANTU_DEBUG_TO_IMPLEMENT(); assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { // MeshUtils::purifyMesh(mesh); getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering); if(mass ) mesh.removeNodesFromArray(*mass , new_numbering); if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if(force ) mesh.removeNodesFromArray(*force , new_numbering); if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); if(blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromArray(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::reassignMaterial() { AKANTU_DEBUG_IN(); std::vector< Array > element_to_add (materials.size()); std::vector< Array > element_to_remove(materials.size()); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; element.ghost_type = ghost_type; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_regular); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type, _ek_regular); for(; it != end; ++it) { ElementType type = *it; element.type = type; element.kind = Mesh::getKind(type); UInt nb_element = mesh.getNbElement(type, ghost_type); Array & el_index_by_mat = element_index_by_material(type, ghost_type); for (UInt el = 0; el < nb_element; ++el) { element.element = el; UInt old_material = el_index_by_mat(el, 0); 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); } } } } std::vector::iterator mat_it; UInt mat_index = 0; for(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::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ internalAddDumpFieldToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "blocked_dofs" ) { ADD_FIELD(blocked_dofs , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField<>(mesh, spatial_dimension, _not_ghost, _ek_regular)); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, spatial_dimension, _not_ghost, _ek_regular)); } else { bool is_internal = false; /// check if at least one material contains field_id as an internal for (UInt m = 0; m < materials.size() && !is_internal; ++m) { is_internal |= materials[m]->isInternal(field_id); } if (is_internal) { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::HomogenizedField (*this, field_id, spatial_dimension, _not_ghost, _ek_regular)); } else { try { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::ElementalField(mesh.getData(field_id), spatial_dimension, _not_ghost, _ek_regular)); } catch (...) {} try { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::ElementalField(mesh.getData(field_id), spatial_dimension, _not_ghost, _ek_regular)); } catch (...) {} } } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ElementTypeMapArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ group.addDumpFieldExternalToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "blocked_dofs") { ADD_FIELD(blocked_dofs , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField(mesh, dim, _not_ghost, _ek_regular, elemental_filter)); } else if(field_id == "element_index_by_material") { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, dim, _not_ghost, _ek_regular, elemental_filter)); } else { typedef DumperIOHelper::HomogenizedField Field; group.addDumpFieldExternalToDumper(dumper_name, field_id, new Field(*this, field_id, dim, _not_ghost, _ek_regular, elemental_filter)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ElementTypeMapArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter); \ f->setPadding(3); \ group.addDumpFieldExternalToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, dim, _not_ghost, _ek_regular, elemental_filter); field->setPadding(3); group.addDumpFieldExternalToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER if(field_id == "stress") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else if(field_id == "grad_u") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #endif } /* ------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name) { EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); synch_registry->synchronize(_gst_for_dump); Dumpable::dump(dumper_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) { EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); synch_registry->synchronize(_gst_for_dump); Dumpable::dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); synch_registry->synchronize(_gst_for_dump); Dumpable::dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump() { Dumpable::dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(UInt step) { Dumpable::dump(step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(Real time, UInt step) { Dumpable::dump(time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeCauchyStresses() { AKANTU_DEBUG_IN(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::saveStressAndStrainBeforeDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::BeginningOfDamageIterationEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateEnergiesAfterDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; element_index_by_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; std::vector::const_iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { const Material & mat = *(*mat_it); mat.printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__