diff --git a/packages/contact_mechanics.cmake b/packages/contact_mechanics.cmake index 57fac5659..7b437511f 100644 --- a/packages/contact_mechanics.cmake +++ b/packages/contact_mechanics.cmake @@ -1,47 +1,54 @@ #=============================================================================== # @file contact_mechanics.cmake # # @author Mohit Pundir # # @date creation: Sun Oct 21 2018 # @date last modification: Sun Oct 21 2018 # # @brief package description for contact mechanics # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(contact_mechanics DESCRIPTION "Use Contact Mechanics package of Akantu") package_declare_sources(contact_mechanics model/contact_mechanics/contact_mechanics_model.hh model/contact_mechanics/contact_mechanics_model.cc model/contact_mechanics/contact_detector.hh model/contact_mechanics/contact_detector.cc - model/contact_mechanics/contact_element.hh) + model/contact_mechanics/contact_element.hh + + model/contact_mechanics/resolution.hh + model/contact_mechanics/resolution.cc + model/contact_mechanics/resolutions/resolution_penalty.hh + model/contact_mechanics/resolutions/resolution_penalty.cc) package_declare_documentation_files(contact_mechanics - manual-contact_detector.tex + manual-contactmechanicsmodel.tex + manual-contact-detector.tex + manual-contact-resolutions.tex ) package_declare_documentation(contact_mechanics - "This package activates the contact ") + "This package activates the contact mechanics model") diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index f5faebfc3..ad815b370 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,580 +1,581 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ using MemoryID = UInt; // using Surface = std::string; // using SurfacePair= std::pair; // using SurfacePairList = std::list; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; #define AKANTU_PP_ENUM(s, data, i, elem) \ BOOST_PP_TUPLE_REM() \ elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data))) } // namespace akantu #if (defined(__GNUC__) || defined(__GNUG__)) #define AKA_GCC_VERSION \ (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) #if AKA_GCC_VERSION < 60000 #define AKANTU_ENUM_HASH(type_name) \ namespace std { \ template <> struct hash<::akantu::type_name> { \ using argument_type = ::akantu::type_name; \ size_t operator()(const argument_type & e) const noexcept { \ auto ue = underlying_type_t(e); \ return uh(ue); \ } \ \ private: \ const hash> uh{}; \ }; \ } #else #define AKANTU_ENUM_HASH(type_name) #endif // AKA_GCC_VERSION #endif // GNU #include "aka_element_classes_info.hh" namespace akantu { #define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) #define AKANTU_PP_TYPE_TO_STR(s, data, elem) \ ({BOOST_PP_CAT(data::_, elem), BOOST_PP_STRINGIZE(elem)}) #define AKANTU_PP_STR_TO_TYPE(s, data, elem) \ ({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data::_, elem)}) #define AKANTU_ENUM_DECLARE(type_name, list) \ enum class type_name { \ BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \ }; #define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ } \ AKANTU_ENUM_HASH(type_name) \ namespace aka { \ inline std::string to_string(const ::akantu::type_name & type) { \ static std::unordered_map<::akantu::type_name, std::string> convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, \ ::akantu::type_name, list))}; \ return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << aka::to_string(type); \ return stream; \ } #define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ inline std::istream & operator>>(std::istream & stream, type_name & type) { \ std::string str; \ stream >> str; \ static std::unordered_map convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, type_name, list))}; \ type = convert.at(str); \ return stream; \ } /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #ifndef SWIG // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) \ (contact_mechanics_model) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { _model, _solid_mechanics_model, _solid_mechanics_model_cohesive, _heat_transfer_model, _structural_mechanics_model, _embedded_model, _contact_mechanics_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_penalty = 4 + _explicit_contact = 4, + _implicit_contact = 5 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent _nls_lumped, ///< Case of lumped mass or equivalent matrix _nls_auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_ask_nodes, _gst_size, //--- SolidMechanicsModel tags --- _gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _gst_smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _gst_smm_uv, ///< synchronization of the nodal velocities and displacement _gst_smm_res, ///< synchronization of the nodal residual _gst_smm_init_mat, ///< synchronization of the data to initialize materials _gst_smm_stress, ///< synchronization of the stresses to compute the internal /// forces _gst_smmc_facets, ///< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, ///< synchronization of facet global connectivity _gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet /// synch _gst_smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _gst_giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _gst_ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gst_gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _gst_htm_temperature, ///< synchronization of the nodal temperature _gst_htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _gst_htm_phi, ///< synchronization of the nodal level set value phi _gst_htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _gst_mnl_for_average, ///< synchronization of data to average in non local /// material _gst_mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _gst_nh_criterion, // --- General tags --- _gst_test, ///< Test tag _gst_user_1, ///< tag for user simulations _gst_user_2, ///< tag for user simulations _gst_material_id, ///< synchronization of the material ids _gst_for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _gst_cf_nodal, ///< synchronization of disp, velo, and current position _gst_cf_incr, ///< synchronization of increment // --- Solver tags --- _gst_solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// @enum ContactElementType type of contact element enum class ContactElementType { _node_to_node, _node_to_segment }; /// @enum ContactDetectorType type of contact detection enum class ContactDetectorType { extrinsic, intrinsic }; /// 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)); } inline std::ostream & operator<<(std::ostream & stream, const NodeFlag & flag) { using under = std::underlying_type_t; stream << under(flag); return stream; } } // namespace akantu #ifndef SWIG AKANTU_ENUM_HASH(GhostType) #endif namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ struct TensorTrait {}; /* -------------------------------------------------------------------------- */ template using is_tensor = std::is_base_of; /* -------------------------------------------------------------------------- */ template using is_scalar = std::is_arithmetic; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ // #if defined(__INTEL_COMPILER) // /// remark #981: operands are evaluated in unspecified order // #pragma warning(disable : 981) // /// remark #383: value copied to temporary, reference to temporary used // #pragma warning(disable : 383) // #endif // defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); 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); /* -------------------------------------------------------------------------- */ } // namespace akantu #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic number * is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/io/parser/parser.hh b/src/io/parser/parser.hh index be92244f2..750d6c7b8 100644 --- a/src/io/parser/parser.hh +++ b/src/io/parser/parser.hh @@ -1,510 +1,511 @@ /** * @file parser.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Fri Dec 08 2017 * * @brief File parser interface * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARSER_HH__ #define __AKANTU_PARSER_HH__ namespace akantu { #ifndef SWIG // clang-format off #define AKANTU_SECTION_TYPES \ (cohesive_inserter) \ (contact) \ (embedded_interface) \ (friction) \ (global) \ (heat) \ (integration_scheme) \ (material) \ (mesh) \ (model) \ (model_solver) \ (neighborhood) \ (neighborhoods) \ (non_linear_solver) \ (non_local) \ (rules) \ (solver) \ (time_step_solver) \ (user) \ (weight_function) \ (contact_detector) \ + (contact_resolution) \ (not_defined) // clang-format on /// Defines the possible section types AKANTU_ENUM_DECLARE(ParserType, AKANTU_SECTION_TYPES) AKANTU_ENUM_OUTPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) AKANTU_ENUM_INPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) /// Defines the possible search contexts/scopes (for parameter search) enum ParserParameterSearchCxt { _ppsc_current_scope = 0x1, _ppsc_parent_scope = 0x2, _ppsc_current_and_parent_scope = 0x3 }; #endif /* ------------------------------------------------------------------------ */ /* Parameters Class */ /* ------------------------------------------------------------------------ */ class ParserSection; /// @brief The ParserParameter objects represent the end of tree branches as /// they /// are the different informations contained in the input file. class ParserParameter { public: ParserParameter() : name(std::string()), value(std::string()), dbg_filename(std::string()) { } ParserParameter(const std::string & name, const std::string & value, const ParserSection & parent_section) : parent_section(&parent_section), name(name), value(value), dbg_filename(std::string()) {} ParserParameter(const ParserParameter & param) = default; virtual ~ParserParameter() = default; /// Get parameter name const std::string & getName() const { return name; } /// Get parameter value const std::string & getValue() const { return value; } /// Set info for debug output void setDebugInfo(const std::string & filename, UInt line, UInt column) { dbg_filename = filename; dbg_line = line; dbg_column = column; } template inline operator T() const; // template inline operator Vector() const; // template inline operator Matrix() const; /// Print parameter info in stream void printself(std::ostream & stream, __attribute__((unused)) unsigned int indent = 0) const { stream << name << ": " << value << " (" << dbg_filename << ":" << dbg_line << ":" << dbg_column << ")"; } private: void setParent(const ParserSection & sect) { parent_section = § } friend class ParserSection; private: /// Pointer to the parent section const ParserSection * parent_section{nullptr}; /// Name of the parameter std::string name; /// Value of the parameter std::string value; /// File for debug output std::string dbg_filename; /// Position of parameter in parsed file UInt dbg_line, dbg_column; }; /* ------------------------------------------------------------------------ */ /* Sections Class */ /* ------------------------------------------------------------------------ */ /// ParserSection represents a branch of the parsing tree. class ParserSection { public: using SubSections = std::multimap; using Parameters = std::map; private: using const_section_iterator_ = SubSections::const_iterator; public: /* ------------------------------------------------------------------------ */ /* SubSection iterator */ /* ------------------------------------------------------------------------ */ /// Iterator on sections class const_section_iterator { public: using iterator_category = std::forward_iterator_tag; using value_type = ParserSection; using pointer = ParserSection *; using reference = ParserSection &; const_section_iterator() = default; const_section_iterator(const const_section_iterator_ & it) : it(it) {} const_section_iterator(const const_section_iterator & other) = default; const_section_iterator & operator=(const const_section_iterator & other) = default; const ParserSection & operator*() const { return it->second; } const ParserSection * operator->() const { return &(it->second); } bool operator==(const const_section_iterator & other) const { return it == other.it; } bool operator!=(const const_section_iterator & other) const { return it != other.it; } const_section_iterator & operator++() { ++it; return *this; } const_section_iterator operator++(int) { const_section_iterator tmp = *this; operator++(); return tmp; } private: const_section_iterator_ it; }; /* ------------------------------------------------------------------------ */ /* Parameters iterator */ /* ------------------------------------------------------------------------ */ /// Iterator on parameters class const_parameter_iterator { public: const_parameter_iterator(const const_parameter_iterator & other) = default; const_parameter_iterator(const Parameters::const_iterator & it) : it(it) {} const_parameter_iterator & operator=(const const_parameter_iterator & other) { if (this != &other) { it = other.it; } return *this; } const ParserParameter & operator*() const { return it->second; } const ParserParameter * operator->() { return &(it->second); }; bool operator==(const const_parameter_iterator & other) const { return it == other.it; } bool operator!=(const const_parameter_iterator & other) const { return it != other.it; } const_parameter_iterator & operator++() { ++it; return *this; } const_parameter_iterator operator++(int) { const_parameter_iterator tmp = *this; operator++(); return tmp; } private: Parameters::const_iterator it; }; /* ---------------------------------------------------------------------- */ ParserSection() : name(std::string()) {} ParserSection(const std::string & name, ParserType type) : parent_section(nullptr), name(name), type(type) {} ParserSection(const std::string & name, ParserType type, const std::string & option, const ParserSection & parent_section) : parent_section(&parent_section), name(name), type(type), option(option) {} ParserSection(const ParserSection & section) : parent_section(section.parent_section), name(section.name), type(section.type), option(section.option), parameters(section.parameters), sub_sections_by_type(section.sub_sections_by_type) { setChldrenPointers(); } ParserSection & operator=(const ParserSection & other) { if (&other != this) { parent_section = other.parent_section; name = other.name; type = other.type; option = other.option; parameters = other.parameters; sub_sections_by_type = other.sub_sections_by_type; setChldrenPointers(); } return *this; } virtual ~ParserSection(); virtual void printself(std::ostream & stream, unsigned int indent = 0) const; /* ---------------------------------------------------------------------- */ /* Creation functions */ /* ---------------------------------------------------------------------- */ public: ParserParameter & addParameter(const ParserParameter & param); ParserSection & addSubSection(const ParserSection & section); protected: /// Clean ParserSection content void clean() { parameters.clear(); sub_sections_by_type.clear(); } private: void setChldrenPointers() { for (auto && param_pair : this->parameters) param_pair.second.setParent(*this); for (auto && sub_sect_pair : this->sub_sections_by_type) sub_sect_pair.second.setParent(*this); } /* ---------------------------------------------------------------------- */ /* Accessors */ /* ---------------------------------------------------------------------- */ public: #ifndef SWIG class SubSectionsRange : public std::pair { public: SubSectionsRange(const const_section_iterator & first, const const_section_iterator & second) : std::pair(first, second) {} auto begin() { return this->first; } auto end() { return this->second; } }; /// Get begin and end iterators on subsections of certain type auto getSubSections(ParserType type = ParserType::_not_defined) const { if (type != ParserType::_not_defined) { auto range = sub_sections_by_type.equal_range(type); return SubSectionsRange(range.first, range.second); } else { return SubSectionsRange(sub_sections_by_type.begin(), sub_sections_by_type.end()); } } /// Get number of subsections of certain type UInt getNbSubSections(ParserType type = ParserType::_not_defined) const { if (type != ParserType::_not_defined) { return this->sub_sections_by_type.count(type); } else { return this->sub_sections_by_type.size(); } } /// Get begin and end iterators on parameters auto getParameters() const { return std::pair( parameters.begin(), parameters.end()); } #endif /* ---------------------------------------------------------------------- */ /// Get parameter within specified context const ParserParameter & getParameter( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if (search_ctx & _ppsc_current_scope) it = parameters.find(name); if (it == parameters.end()) { if ((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->getParameter(name, search_ctx); else { AKANTU_SILENT_EXCEPTION( "The parameter " << name << " has not been found in the specified context"); } } return it->second; } /* ------------------------------------------------------------------------ */ /// Get parameter within specified context, with a default value in case the /// parameter does not exists template const T getParameter( const std::string & name, const T & default_value, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { try { T tmp = this->getParameter(name, search_ctx); return tmp; } catch (debug::Exception &) { return default_value; } } /* ------------------------------------------------------------------------ */ /// Check if parameter exists within specified context bool hasParameter( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if (search_ctx & _ppsc_current_scope) it = parameters.find(name); if (it == parameters.end()) { if ((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->hasParameter(name, search_ctx); else { return false; } } return true; } /* -------------------------------------------------------------------------- */ /// Get value of given parameter in context template T getParameterValue( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { const ParserParameter & tmp_param = getParameter(name, search_ctx); T t = tmp_param; return t; } /* -------------------------------------------------------------------------- */ /// Get section name const std::string getName() const { return name; } /// Get section type ParserType getType() const { return type; } /// Get section option const std::string getOption(const std::string & def = "") const { return option != "" ? option : def; } protected: void setParent(const ParserSection & sect) { parent_section = § } /* ---------------------------------------------------------------------- */ /* Members */ /* ---------------------------------------------------------------------- */ private: /// Pointer to the parent section const ParserSection * parent_section{nullptr}; /// Name of section std::string name; /// Type of section, see AKANTU_SECTION_TYPES ParserType type{ParserType::_not_defined}; /// Section option std::string option; /// Map of parameters in section Parameters parameters; /// Multi-map of subsections SubSections sub_sections_by_type; }; /* ------------------------------------------------------------------------ */ /* Parser Class */ /* ------------------------------------------------------------------------ */ /// Root of parsing tree, represents the global ParserSection class Parser : public ParserSection { public: Parser() : ParserSection("global", ParserType::_global) {} void parse(const std::string & filename); std::string getLastParsedFile() const; static bool isPermissive() { return permissive_parser; } public: /// Parse real scalar static Real parseReal(const std::string & value, const ParserSection & section); /// Parse real vector static Vector parseVector(const std::string & value, const ParserSection & section); /// Parse real matrix static Matrix parseMatrix(const std::string & value, const ParserSection & section); #ifndef SWIG /// Parse real random parameter static RandomParameter parseRandomParameter(const std::string & value, const ParserSection & section); #endif protected: /// General parse function template static T parseType(const std::string & value, Grammar & grammar); protected: // friend class Parsable; static bool permissive_parser; std::string last_parsed_file; }; inline std::ostream & operator<<(std::ostream & stream, const ParserParameter & _this) { _this.printself(stream); return stream; } inline std::ostream & operator<<(std::ostream & stream, const ParserSection & section) { section.printself(stream); return stream; } } // akantu namespace std { template <> struct iterator_traits<::akantu::Parser::const_section_iterator> { using iterator_category = input_iterator_tag; using value_type = ::akantu::ParserParameter; using difference_type = ptrdiff_t; using pointer = const ::akantu::ParserParameter *; using reference = const ::akantu::ParserParameter &; }; } #include "parser_tmpl.hh" #endif /* __AKANTU_PARSER_HH__ */ diff --git a/src/model/boundary_condition_functor_inline_impl.cc b/src/model/boundary_condition_functor_inline_impl.cc index 41a7d2b87..3a23bf6c4 100644 --- a/src/model/boundary_condition_functor_inline_impl.cc +++ b/src/model/boundary_condition_functor_inline_impl.cc @@ -1,158 +1,150 @@ /** * @file boundary_condition_functor_inline_impl.cc * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Mon Feb 19 2018 * * @brief implementation of the BC::Functors * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition_functor.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_BOUNDARY_CONDITION_FUNCTOR_INLINE_IMPL_CC__ #define __AKANTU_BOUNDARY_CONDITION_FUNCTOR_INLINE_IMPL_CC__ /* -------------------------------------------------------------------------- */ #define DIRICHLET_SANITY_CHECK \ - AKANTU_DEBUG_ASSERT( \ - coord.size() <= flags.size(), \ - "The coordinates and flags vectors given to the boundary" \ - << " condition functor have different sizes!"); \ - AKANTU_DEBUG_ASSERT( \ - primal.size() <= coord.size(), \ - "The primal vector and coordinates vector given" \ + AKANTU_DEBUG_ASSERT( \ + primal.size() <= flags.size(), \ + "The primal vector and flags vectors given" \ << " to the boundary condition functor have different sizes!"); #define NEUMANN_SANITY_CHECK \ AKANTU_DEBUG_ASSERT( \ coord.size() <= normals.size(), \ "The coordinates and normals vectors given to the" \ << " boundary condition functor have different sizes!"); \ - AKANTU_DEBUG_ASSERT( \ - dual.size() <= coord.size(), \ - "The dual vector and coordinates vector given to" \ - << " the boundary condition functor have different sizes!"); - + namespace akantu { namespace BC { namespace Dirichlet { /* ---------------------------------------------------------------------- */ inline void FlagOnly:: operator()(__attribute__((unused)) UInt node, Vector & flags, __attribute__((unused)) Vector & primal, __attribute__((unused)) const Vector & coord) const { DIRICHLET_SANITY_CHECK; flags(this->axis) = true; } /* ---------------------------------------------------------------------- */ inline void FreeBoundary:: operator()(__attribute__((unused)) UInt node, Vector & flags, __attribute__((unused)) Vector & primal, __attribute__((unused)) const Vector & coord) const { DIRICHLET_SANITY_CHECK; flags(this->axis) = false; } /* ---------------------------------------------------------------------- */ inline void FixedValue::operator()(__attribute__((unused)) UInt node, Vector & flags, Vector & primal, __attribute__((unused)) const Vector & coord) const { DIRICHLET_SANITY_CHECK; flags(this->axis) = true; primal(this->axis) = value; } /* ---------------------------------------------------------------------- */ inline void IncrementValue::operator()(__attribute__((unused)) UInt node, Vector & flags, Vector & primal, __attribute__((unused)) const Vector & coord) const { DIRICHLET_SANITY_CHECK; flags(this->axis) = true; primal(this->axis) += value; } /* ---------------------------------------------------------------------- */ inline void Increment::operator()(__attribute__((unused)) UInt node, Vector & flags, Vector & primal, __attribute__((unused)) const Vector & coord) const { DIRICHLET_SANITY_CHECK; flags.set(true); primal += value; } } // end namespace Dirichlet /* ------------------------------------------------------------------------ */ /* Neumann */ /* ------------------------------------------------------------------------ */ namespace Neumann { /* ---------------------------------------------------------------------- */ inline void FreeBoundary:: operator()(__attribute__((unused)) const IntegrationPoint & quad_point, Vector & dual, __attribute__((unused)) const Vector & coord, __attribute__((unused)) const Vector & normals) const { for (UInt i(0); i < dual.size(); ++i) { dual(i) = 0.0; } } /* ---------------------------------------------------------------------- */ inline void FromHigherDim::operator()(__attribute__((unused)) const IntegrationPoint & quad_point, Vector & dual, __attribute__((unused)) const Vector & coord, const Vector & normals) const { dual.mul(this->bc_data, normals); } /* ---------------------------------------------------------------------- */ inline void FromSameDim:: operator()(__attribute__((unused)) const IntegrationPoint & quad_point, Vector & dual, __attribute__((unused)) const Vector & coord, __attribute__((unused)) const Vector & normals) const { dual = this->bc_data; } } // namespace Neumann } // namespace BC } // namespace akantu #endif /* __AKANTU_BOUNDARY_CONDITION_FUNCTOR_INLINE_IMPL_CC__ */ diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 4a1a2c7f8..a87f56863 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,495 +1,506 @@ /** * @file contact_detector.cc * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Fri Sep 21 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "contact_detector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, const ID & id, UInt memory_id) : ContactDetector(mesh, mesh.getNodes(), id, memory_id) { } /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, Array & positions, const ID & id, UInt memory_id) : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), mesh(mesh), positions(positions) { AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); this->mesh.fillNodesToElements(this->spatial_dimension - 1); this->parseSection(); this->getMaximalDetectionDistance(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactDetector::parseSection() { const Parser & parser = getStaticParser(); const ParserSection & section = *(parser.getSubSections(ParserType::_contact_detector).first); auto type = section.getParameterValue("type"); if (type == "intrinsic") { this->detection_type = ContactDetectorType::intrinsic; } else if (type == "extrinsic"){ this->detection_type = ContactDetectorType::extrinsic; } else { AKANTU_ERROR("Unknown detection type : " << type); } surfaces[Surface::master] = section.getParameterValue("master_surface"); surfaces[Surface::slave ] = section.getParameterValue("slave_surface"); //this->max_bb = section.getParameterValue("max_bb"); } /* -------------------------------------------------------------------------- */ void ContactDetector::getMaximalDetectionDistance() { AKANTU_DEBUG_IN(); Real el_size; Real max_el_size = std::numeric_limits::min(); auto & master_group = mesh.getElementGroup(surfaces[Surface::master]); for (auto & type: master_group.elementTypes(spatial_dimension - 1, _not_ghost)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array coord(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); UInt nb_element = mesh.getNbElement(type); for (UInt el =0; el < nb_element; ++el, ++el_coord) { el_size = FEEngine::getElementInradius(*el_coord, type); max_el_size = std::max(max_el_size, el_size); } AKANTU_DEBUG_INFO("The maximum element size : " << max_el_size ); } this->max_dd = max_el_size; this->max_bb = max_el_size; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void ContactDetector::search(std::vector & elements) { - this->globalSearch(elements); + void ContactDetector::search(std::map & contact_map) { + this->globalSearch(contact_map); } /* -------------------------------------------------------------------------- */ -void ContactDetector::globalSearch(std::vector & elements) { +void ContactDetector::globalSearch(std::map & contact_map) { auto & master_list = mesh.getElementGroup(surfaces[Surface::master]).getNodeGroup().getNodes(); auto & slave_list = mesh.getElementGroup(surfaces[Surface::slave]).getNodeGroup().getNodes(); BBox bbox_master(spatial_dimension); this->constructBoundingBox(bbox_master, master_list); BBox bbox_slave(spatial_dimension); this->constructBoundingBox(bbox_slave, slave_list); auto && bbox_intersection = bbox_master.intersection(bbox_slave); AKANTU_DEBUG_INFO( "Intersection BBox " << bbox_intersection ); Vector center(spatial_dimension); bbox_intersection.getCenter(center); Vector spacing(spatial_dimension); this->computeCellSpacing(spacing); auto & master_surface_list = mesh.getElementGroup(surfaces[Surface::master]).getNodeGroup().getNodes(); auto & slave_surface_list = mesh.getElementGroup(surfaces[Surface::slave]).getNodeGroup().getNodes(); SpatialGrid master_grid(spatial_dimension, spacing, center); this->constructGrid(master_grid, bbox_intersection, master_surface_list); SpatialGrid slave_grid(spatial_dimension, spacing, center); this->constructGrid(slave_grid, bbox_intersection, slave_surface_list); if (AKANTU_DEBUG_TEST(dblDump)) { Mesh mesh(spatial_dimension, "save"); master_grid.saveAsMesh(mesh); mesh.write("master_grid.msh"); } if (AKANTU_DEBUG_TEST(dblDump)) { Mesh mesh2(spatial_dimension, "save"); slave_grid.saveAsMesh(mesh2); mesh2.write("slave_grid.msh"); } AKANTU_DEBUG_INFO( "Grid Details " << master_grid ); // search slave grid nodes in contactelement array and if they exits // and still have orthogonal projection on its associated master // facetremove it from the spatial grid or do not consider it for // local search, maybe better option will be to have spatial grid of // type node info and one of the variable of node info should be // facet already exits // so contact eleemnts will be updated based on the above // consideration , this means only those contact elements will be // keep whose slave node is still in intersection bbox and still has // projection in its master facet // also if slave node is already exists in contact element and // orthogonal projection does not exits then search the associated // master facets with the current master facets within a given // radius , this is subjected to computational cost as searching // neighbbor cells can be more effective. - this->localSearch(slave_grid, master_grid, elements); + this->localSearch(slave_grid, master_grid, contact_map); } /* -------------------------------------------------------------------------- */ void ContactDetector::localSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid, - std::vector & contact_elements) { + std::map & contact_map) { // local search // out of these array check each cell for closet node in that cell // and neighbouring cells find the actual orthogonally closet // check the projection of slave node on master facets connected to // the closet master node, if yes update the contact element with // slave node and master node and master surfaces connected to the // master node // these master surfaces will be needed later to update contact // elements Array slave_nodes; Array master_nodes; BBox bbox_master_grid(spatial_dimension); BBox bbox_slave_grid(spatial_dimension); bbox_master_grid += master_grid.getUpperBounds(); bbox_master_grid += master_grid.getLowerBounds(); bbox_slave_grid += slave_grid.getUpperBounds(); bbox_slave_grid += slave_grid.getLowerBounds(); auto && bbox_intersection = bbox_master_grid.intersection(bbox_slave_grid); // find the closet master node for each slave node for (auto && cell_id : slave_grid) { AKANTU_DEBUG_INFO("Looping on next cell"); for (auto && q1: slave_grid.getCell(cell_id)) { Vector pos(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos(s) = this->positions(q1, s); } if (!bbox_intersection.contains(pos)) { continue; } Real closet_distance = std::numeric_limits::max(); UInt closet_master_node; // loop over all the neighboring cells of the current cells for (auto && neighbor_cell : cell_id.neighbors()) { // loop over the data of neighboring cells from master grid for (auto && q2 : master_grid.getCell(neighbor_cell)) { AKANTU_DEBUG_INFO("Looping on neighbor cell in master"); Vector pos2(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos2(s) = this->positions(q2, s); } Real distance = pos.distance(pos2); if (distance <= closet_distance) { closet_master_node = q2; closet_distance = distance; } } } slave_nodes.push_back(q1); master_nodes.push_back(closet_master_node); } } for (auto && values : zip(slave_nodes, master_nodes)) { const auto & slave_node = std::get<0>(values); const auto & master_node = std::get<1>(values); Array elements; this->mesh.getAssociatedElements(master_node, elements); auto normals = std::make_unique>(elements.size(), spatial_dimension, "normals"); auto gaps = std::make_unique>(elements.size(), 1, "gaps"); this->computeOrthogonalProjection(slave_node, elements, *normals, *gaps); auto minimum = std::min_element( gaps->begin(), gaps->end()); auto index = std::distance( gaps->begin(), minimum); auto normal = Vector(normals->begin(spatial_dimension)[index]); - auto contact_element = ContactElement(slave_node, elements[index]); - contact_element.setGap( (*gaps)[index]); - contact_element.setNormal( normal); + //auto contact_element = ContactElement(slave_node, elements[index]); + //contact_element.setGap( (*gaps)[index]); + //contact_element.setNormal( normal); + + //const Array & connectivity = this->mesh.getConnectivity(elements[index].type, ghost_type); + //contact_map[slave].setConnectivity(connectivity(elements[index].element)); + + contact_map[slave_node] = ContactElement(slave_node, elements[index]); + contact_map[slave_node].gap = (*gaps)[index]; + contact_map[slave_node].normal = normal; - contact_elements.push_back(contact_element); } } /* -------------------------------------------------------------------------- */ void ContactDetector::constructGrid(SpatialGrid & grid, BBox & bbox, const Array & nodes_list) { auto to_grid = [&](UInt node) { Vector pos(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos(s) = this->positions(node, s); } if (bbox.contains(pos)) { grid.insert(node, pos); } }; std::for_each(nodes_list.begin(), nodes_list.end(), to_grid); } /* -------------------------------------------------------------------------- */ void ContactDetector::constructBoundingBox(BBox & bbox, const Array & nodes_list) { auto to_bbox = [&](UInt node) { Vector pos(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos(s) = this->positions(node, s); } bbox += pos; }; std::for_each(nodes_list.begin(), nodes_list.end(), to_bbox); AKANTU_DEBUG_INFO("BBox" << bbox); auto & lower_bound = bbox.getLowerBounds(); auto & upper_bound = bbox.getUpperBounds(); for (UInt s: arange(spatial_dimension)) { lower_bound(s) -= this->max_bb; upper_bound(s) += this->max_bb; } AKANTU_DEBUG_INFO("BBox" << bbox); } /* -------------------------------------------------------------------------- */ void ContactDetector::computeCellSpacing(Vector & spacing) { for (UInt s: arange(spatial_dimension)) spacing(s) = std::sqrt(2.0) * max_dd; } /* -------------------------------------------------------------------------- */ void ContactDetector::computeOrthogonalProjection(const UInt & node, const Array & elements, Array & normals, Array & gaps) { Vector query(spatial_dimension); for (UInt s: arange(spatial_dimension)) { query(s) = this->positions(node, s); } for (auto && values : zip( elements, gaps, make_view(normals , spatial_dimension))) { const auto & element = std::get<0>(values); auto & gap = std::get<1>(values); auto & normal = std::get<2>(values); Vector projection(spatial_dimension); this->computeNormalOnElement(element, normal); this->computeProjectionOnElement(element, normal, query, projection); projection -= query; gap = Math::norm(spatial_dimension, projection.storage()); } } /* -------------------------------------------------------------------------- */ void ContactDetector::computeNormalOnElement(const Element & element, Vector & normal) { Matrix vectors(spatial_dimension, spatial_dimension - 1); this->vectorsAlongElement(element, vectors); switch (this->spatial_dimension) { case 2: { Math::normal2(vectors.storage(), normal.storage()); break; } case 3: { Math::normal3(vectors(0).storage(), vectors(1).storage(), normal.storage()); break; } default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } } switch (this->detection_type) { case ContactDetectorType::extrinsic: { normal *= -1.0; break; } default: break; } } /* -------------------------------------------------------------------------- */ void ContactDetector::computeProjectionOnElement(const Element & element, const Vector & normal, const Vector & query, Vector & projection) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); Real alpha = (query - barycenter).dot(normal); projection = query - alpha * normal; + // use contains function to check whether projection lies inside + // the element, if yes it is a valid projection otherwise no + // still have to think about what to do if normal exists but + // projection doesnot lie inside the element + bool validity = this->isValidProjection(element, projection); if (!validity) { projection *= std::numeric_limits::max(); } } /* -------------------------------------------------------------------------- */ bool ContactDetector::isValidProjection(const Element & element, Vector & projection) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); Real distance = 0; switch (this->spatial_dimension) { case 2: { distance = Math::distance_2d(projection.storage(), barycenter.storage()); break; } case 3: { distance = Math::distance_3d(projection.storage(), barycenter.storage()); break; } default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } } if (distance <= max_dd) { return true; } return false; } /* -------------------------------------------------------------------------- */ void ContactDetector::coordinatesOfElement(const Element & el, Matrix & coords) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Vector connect = mesh.getConnectivity(el.type, _not_ghost) .begin(nb_nodes_per_element)[el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; for (UInt s: arange(spatial_dimension)) { coords(s, n) = this->positions(node, s); } } } /* -------------------------------------------------------------------------- */ void ContactDetector::vectorsAlongElement(const Element & el, Matrix & vectors) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Matrix coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(el, coords); switch (spatial_dimension) { case 2: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); break; } case 3: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); vectors(1) = Vector(coords(2)) - Vector(coords(0)); break; } default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } } } } // akantu diff --git a/src/model/contact_mechanics/contact_detector.hh b/src/model/contact_mechanics/contact_detector.hh index c19ecdb12..7e08bdbb4 100644 --- a/src/model/contact_mechanics/contact_detector.hh +++ b/src/model/contact_mechanics/contact_detector.hh @@ -1,165 +1,165 @@ /** * @file contact_detection.hh * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Tue Oct 23 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_grid_dynamic.hh" #include "aka_bbox.hh" #include "mesh.hh" #include "mesh_io.hh" #include "fe_engine.hh" #include "parsable.hh" #include "element_group.hh" #include "contact_element.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_HH__ #define __AKANTU_CONTACT_DETECTOR_HH__ namespace akantu { enum class Surface { master, slave }; /* -------------------------------------------------------------------------- */ class ContactDetector : private Memory, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: ContactDetector(Mesh &, const ID & id = "contact_detector", UInt memory_id = 0); ContactDetector(Mesh &, Array & positions, const ID & id = "contact_detector", UInt memory_id = 0); ~ContactDetector() = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: /// - void search(std::vector &); + void search(std::map &); /// computes orthogonal projection on master elements void computeOrthogonalProjection(const UInt & /* slave node */, const Array & /* master elements */, Array & /* normals */, Array & /* projections */); private: /// reads the input file to get contact detection options void parseSection(); /// performs global spatial search - void globalSearch(std::vector &); + void globalSearch(std::map &); /// performs local search to create contact element /// TODO: templated function typename void localSearch(SpatialGrid &, SpatialGrid &, - std::vector &); + std::map &); /// constructs a grid containing nodes lying within bounding box /// TODO : templated fucntion to created template Spatial Grid void constructGrid(SpatialGrid &, BBox &, const Array &); /// constructs the bounding box based on nodes list void constructBoundingBox(BBox &, const Array &); /// computes the optimal cell size for grid void computeCellSpacing(Vector &); /// computes the maximum in radius for a given mesh void getMaximalDetectionDistance(); /// extracts the coordinates of an element void coordinatesOfElement(const Element & /* element id */, Matrix & /* coordinates */); /// extracts vectors which forms the plane of element void vectorsAlongElement(const Element & /* element id */, Matrix & /* vectors matrix */); /// computes normal on an element void computeNormalOnElement(const Element & /* element id */, Vector & /* normal vector */); /// computes projection of a query point on an element void computeProjectionOnElement(const Element & /* element */, const Vector & /* normal */, const Vector & /* query */, Vector & /* projection */); /// checks for the validity of a projection bool isValidProjection(const Element & /* element */, Vector & /* projection */); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// maximal detection distance for grid spacing Real max_dd; /// maximal bounding box extension Real max_bb; /// Mesh Mesh & mesh; /// dimension of the model UInt spatial_dimension{0}; /// map to contain ids for surfaces std::map surfaces; /// contains the updated positions of the nodes Array & positions; /// type of detection extrinisic/intrinsic ContactDetectorType detection_type; }; } // namespace akantu #endif /* __AKANTU_CONTACT_DETECTOR_HH__ */ diff --git a/src/model/contact_mechanics/contact_element.hh b/src/model/contact_mechanics/contact_element.hh index 2b1b34a55..9a2802167 100644 --- a/src/model/contact_mechanics/contact_element.hh +++ b/src/model/contact_mechanics/contact_element.hh @@ -1,113 +1,117 @@ /** * @file contact_element.hh * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Tue Oct 23 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" - /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_ELEMENT_HH__ #define __AKANTU_CONTACT_ELEMENT_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { using SlaveType = UInt; using MasterType = Element; class ContactElement { /* ------------------------------------------------------------------------ */ /* Constructor/ Destructors */ /* ------------------------------------------------------------------------ */ public: ContactElement(SlaveType slave, MasterType master) : node(slave), master(master) { + + } ~ContactElement() = default; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ + public: /// sets the value of normal vector AKANTU_SET_MACRO(Normal, normal, Vector); /// sets the value of tangent vector AKANTU_SET_MACRO(Tangent, tangent, Vector); /// sets the value of gap AKANTU_SET_MACRO(Gap, gap, Real); /// gets the value of gap AKANTU_GET_MACRO(Gap, gap, Real); /// gets the normal to the master element AKANTU_GET_MACRO(Normal, normal, Vector); // gets the node of slave AKANTU_GET_MACRO(Slave, node, SlaveType); - // sets the value of normal vector AKANTU_SET_MACRO(Patch, patch, Array); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: /// slave node SlaveType node; /// master element/node MasterType master; /// normalized normal direction Vector normal; /// normalized tangent direction Vector tangent; + /// + Vector connectivity; + /// penetration gap between slave and master Real gap; /// an array of master nodes/elements around slave node Array patch; }; } // akantu #endif /* __AKANTU_CONTACT_ELEMENT_HH__ */ diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 29aaa564d..802edf6af 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,68 +1,265 @@ /** * @file coontact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Contact mechanics model * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ -#include "fe_engine.hh" #include "contact_mechanics_model.hh" +#include "integrator_gauss.hh" +#include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id) { AKANTU_DEBUG_IN(); + this->registerFEEngineObject("ContactMechanicsModel", mesh, + Model::spatial_dimension); + this->initDOFManager(); - //this->detector = std::make_unique( - // this->mesh, id + ":contact_detector"); + this->detector = std::make_unique(this->mesh, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } - +/* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { + + Model::initFullImpl(options); + + // initalize the resolutions + if (this->parser.getLastParsedFile() != "") { + this->instantiateResolutions(); + } + + this->initResolutions(); + +} + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::instantiateResolutions() { + ParserSection model_section; + bool is_empty; + std::tie(model_section, is_empty) = this->getParserSection(); + + if (not is_empty) { + auto model_resolutions = model_section.getSubSections(ParserType::_contact_resolution); + for (const auto & section : model_resolutions) { + this->registerNewResolution(section); + } + } + + auto sub_sections = this->parser.getSubSections(ParserType::_contact_resolution); + for (const auto & section : sub_sections) { + this->registerNewResolution(section); + } + + if (resolutions.empty()) + AKANTU_EXCEPTION("No contact resolutions where instantiated for the model" + << getID()); +} + +/* -------------------------------------------------------------------------- */ +Resolution & +ContactMechanicsModel::registerNewResolution(const ParserSection & section) { + std::string res_name; + std::string res_type = section.getName(); + std::string opt_param = section.getOption(); + + try { + std::string tmp = section.getParameter("name"); + res_name = tmp; /** this can seem weird, but there is an ambiguous operator + * overload that i couldn't solve. @todo remove the + * weirdness of this code + */ + } catch (debug::Exception &) { + AKANTU_ERROR("A contact resolution of type \'" + << res_type + << "\' in the input file has been defined without a name!"); + } + Resolution & res = this->registerNewResolution(res_name, res_type, opt_param); + + res.parseSection(section); + + return res; +} + +/* -------------------------------------------------------------------------- */ +Resolution & ContactMechanicsModel::registerNewResolution(const ID & res_name, + const ID & res_type, + const ID & opt_param) { + AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) == + resolutions_names_to_id.end(), + "A resolution with this name '" + << res_name << "' has already been registered. " + << "Please use unique names for resolutions"); + + UInt res_count = resolutions.size(); + resolutions_names_to_id[res_name] = res_count; + + std::stringstream sstr_res; + sstr_res << this->id << ":" << res_count << ":" << res_type; + ID res_id = sstr_res.str(); + + std::unique_ptr resolution = ResolutionFactory::getInstance().allocate( + res_type, spatial_dimension, opt_param, *this, res_id); + + resolutions.push_back(std::move(resolution)); + + return *(resolutions.back()); +} + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::initResolutions() { + AKANTU_DEBUG_ASSERT(resolutions.size() != 0, "No resolutions to initialize !"); + + if (!are_resolutions_instantiated) + instantiateResolutions(); + + // \TODO check if each resolution needs a initResolution() method +} + +/* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initModel() { + getFEEngine().initShapeFunctions(_not_ghost); + getFEEngine().initShapeFunctions(_ghost); +} + + +/* -------------------------------------------------------------------------- */ +FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { + return dynamic_cast( + getFEEngineClassBoundary(name)); +} + +/* -------------------------------------------------------------------------- */ +std::tuple +ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { + + /* @todo contact mechanics model doesnt really needs a solver have + to find a way to fix this absurd part + */ + return std::make_tuple("contact", _tsst_static); } + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::assembleResidual() { + AKANTU_DEBUG_IN(); + + /* ------------------------------------------------------------------------ */ + // computes the internal forces + this->assembleInternalForces(); + + /* ------------------------------------------------------------------------ */ + this->getDOFManager().assembleToResidual("displacement", + *this->contact_force, 1); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::assembleInternalForces() { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_INFO("Assemble the contact forces"); + + // assemble the forces due to local stresses + AKANTU_DEBUG_INFO("Assemble residual for local elements"); + for (auto & resolution : resolutions) { + resolution->assembleInternalForces(_not_ghost); + } + // assemble the stresses due to ghost elements + AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); + for (auto & resolution : resolutions) { + resolution->assembleInternalForces(_ghost); + } + AKANTU_DEBUG_OUT(); + +} + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::printself(std::ostream & stream, int indent) const { + std::string space; + for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) + ; + + stream << space << "Contact Mechanics Model [" << std::endl; + stream << space << " + id : " << id << std::endl; + stream << space << " + spatial dimension : " << Model::spatial_dimension + << std::endl; + stream << space << " + fem [" << std::endl; + getFEEngine().printself(stream, indent + 2); + stream << space << AKANTU_INDENT << "]" << std::endl; + + stream << space << " + resolutions [" << std::endl; + for (auto & resolution : resolutions) { + resolution->printself(stream, indent + 1); + } + stream << space << AKANTU_INDENT << "]" << std::endl; + + stream << space << "]" << std::endl; +} + + +/* -------------------------------------------------------------------------- */ +MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) { + // \TODO correct it for contact mechanics model, only one type of matrix + if (matrix_id == "C") + return _mt_not_defined; + + return _symmetric; +} + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) { + AKANTU_TO_IMPLEMENT(); +} + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { + AKANTU_TO_IMPLEMENT(); } diff --git a/src/model/contact_mechanics/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh index 33a9d97ef..d116b4f7b 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,127 +1,221 @@ /** * @file contact_mechanics_model.hh * * @author Mohit Pundir * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Contact Mechanics * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" +#include "fe_engine.hh" #include "contact_detector.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__ #define __AKANTU_CONTACT_MECHANICS_MODEL_HH__ namespace akantu { +class Resolution; +template +class IntegratorGauss; +template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class ContactMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ - + + using MyFEEngineType = FEEngineTemplate; public: ContactMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_contact_mechanics_model); ~ContactMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: + /// initialize completely the model + void initFullImpl(const ModelOptions & options) override; + /// initialize all internal arrays for resolutions + void initResolutions(); + /// initialize the modelType void initModel() override; /// computes the force residual void assembleResidual() override; - /// callback to assemble a Matrix - void assembleMatrix(const ID &) override; + /// get the type of matrix needed + MatrixType getMatrixType(const ID & matrix_id) override; + + /// callback for the solver, this assembles different matrices + void assembleMatrix(const ID & matrix_id) override; + + /// callback for the solver, this assembles the stiffness matrix + void assembleLumpedMatrix(const ID & matrix_id) override; + + /// get some default values for derived classes + std::tuple + getDefaultSolverID(const AnalysisMethod & method) override; + /// function to print the containt of the class + void printself(std::ostream & stream, int indent = 0) const override; + + /* ------------------------------------------------------------------------ */ + /* Contact Detection */ + /* ------------------------------------------------------------------------ */ +public: + void search(); + /* ------------------------------------------------------------------------ */ + /* Contact Resolution */ + /* ------------------------------------------------------------------------ */ +public: + /// register an empty contact resolution of a given type + Resolution & registerNewResolution(const ID & res_name, const ID & res_type, + const ID & opt_param); + +protected: + /// register a resolution in the dynamic database + Resolution & registerNewResolution(const ParserSection & res_section); + + /// read the resolution files to instantiate all the resolutions + void instantiateResolutions(); + + /* ------------------------------------------------------------------------ */ + /* Solver Interface */ + /* ------------------------------------------------------------------------ */ +public: + /// assembles the contact stiffness matrix + //virtual void assembleStiffnessMatrix(); + + /// assembles the contant internal forces + virtual void assembleInternalForces(); + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +protected: + + FEEngine & getFEEngineBoundary(const ID & name = "") override; + /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: - dumper::Field * createNodalFieldReal(const std::string & field_name, + /*dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); - virtual void dump(Real time, UInt step); + virtual void dump(Real time, UInt step);*/ - +protected: + /// contact detection class friend class ContactDetector; - -private: + + /// contact resolution class + friend class Resolution; + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + /// get the ConatctMechanics::contact_force vector (internal forces) + AKANTU_GET_MACRO(InternalForce, *contact_force, Array &); + + /// get the contact element for a given slave node + inline ContactElement & getContactElement(const UInt slave) { return contact_map[slave]; } + /// get the contat map + inline std::map & getContactMap() { return contact_map; } + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: + /// tells if the resolutions are instantiated + bool are_resolutions_instantiated; + + /// contact forces array + Array * contact_force{nullptr}; + /// contact detection std::unique_ptr detector; -}; - + /// list of contact resolutions + std::vector> resolutions; + /// mapping between resolution name and resolution internal id + std::map resolutions_names_to_id; + /// mapping between slave node its respective contact element + std::map contact_map; +}; + } // namespace akantu +/* ------------------------------------------------------------------------ */ +/* inline functions */ +/* ------------------------------------------------------------------------ */ +#include "resolution.hh" +#include "parser.hh" +/* ------------------------------------------------------------------------ */ + -#endif +#endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc new file mode 100644 index 000000000..a4ad14d3d --- /dev/null +++ b/src/model/contact_mechanics/resolution.cc @@ -0,0 +1,281 @@ +/** + * @file resolution.cc + * + * @author Mohit Pundir + * + * @date creation: Mon Jan 7 2019 + * @date last modification: Mon Jan 7 2019 + * + * @brief Implementation of common part of the contact resolution class + * + * @section LICENSE + * + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "resolution.hh" +#include "contact_mechanics_model.hh" +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +Resolution::Resolution(ContactMechanicsModel & model, const ID & id) + : Memory(id, model.getMemoryID()), Parsable(ParserType::_contact_resolution, id), + fem(model.getFEEngine()), + name(""), model(model), + spatial_dimension(model.getMesh().getSpatialDimension()){ + + AKANTU_DEBUG_IN(); + + this->initialize(); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +Resolution::~Resolution() = default; + +/* -------------------------------------------------------------------------- */ +void Resolution::initialize() { + registerParam("name", name, std::string(), _pat_parsable | _pat_readable); + registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable, + "Friciton Coefficient"); +} + +/* -------------------------------------------------------------------------- */ +void Resolution::printself(std::ostream & stream, int indent) const { + std::string space; + for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) + ; + + std::string type = getID().substr(getID().find_last_of(':') + 1); + + stream << space << "Contact Resolution " << type << " [" << std::endl; + Parsable::printself(stream, indent); + stream << space << "]" << std::endl; +} + +/* -------------------------------------------------------------------------- */ +void Resolution::assembleInternalForces(GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + const Array & equation_array = + model.getDOFManager().getEquationsNumbers(); + + auto & internal_force = + const_cast &>(model.getInternalForce()); + + auto & contact_map = model.getContactMap(); + + const auto slave_nodes = model.getMesh().getNodeGroup(name); + + for (auto & slave: slave_nodes) { + + auto & master = contact_map[slave].master; + auto & gap = contact_map[slave].gap; + auto & projection = contact_map[slave].projection; + + UInt nb_nodes_master = Mesh::getNbNodesPerElement(master.type) + + Vector shapes(nb_nodes_master); + fem.computeShapes(projection, master, master.type, shapes, ghost_type); + + Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); + fem.computeShapeDerivatives(projection, master, master.type, shapes_derivatives, ghost_type); + + const auto & connectivity = contact_map[slave].connectivity; + Vector elementary_force(connectivity.size() * spatial_dimension); + + Array * tangents = + new Array(spatial_dimension - 1, spatial_dimension, "surface_tangents"); + + Array * global_coords = + new Array(nb_nodes_master, spatial_dimension); + + computeCoordinates(master.type, *global_coords); + computeTangents(shapes_derivatives, *global_coords, *tangents); + + Matrix surface_matrix(spatial_dimension - 1, spatial_dimension - 1); + computeSurfaceMatrix(*tangents, surface_matrix); + + computeN(*n, shapes, normal); + computeNormalForce(elementary_force, *n, gap); + + computeTalpha(*t_alpha, shapes, *tangents); + computeNalpha(*n_alpha, *shapes_derivatives, normal); + computeDalpha(*d_alpha, *n_alpha, *t_alpha, surface_matrix); + computeFrictionForce(elementary_force, *d_alpha, gap); + + for (UInt i = 0; i < connectivity.size(); ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) { + UInt offset_node = connectivity(i) * spatial_dimension + j; + auto & equation_num = equation_array(offset_node); + internal_force(equation_num) += elementary_force(i + j); + } + } + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void Resolution::assembleStiffnessMatrix(GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + const auto slave_nodes = model.getMesh().getNodeGroup(name); + auto & contact_map = model.getContactMap(); + + for (auto & slave: slave_nodes) { + + auto & master = contact_map[slave].master; + auto & gap = contact_map[slave].gap; + auto & projection = contact_map[slave].projection; + + Vector shapes(master.nb_nodes); + fem.computeShapes(projection, master, master.type, shapes, ghost_type); + + Vector shapes_derivatives(master.nb_nodes * spatial_dimension); + fem.computeShapeDerivatives(projection, master, master.type, shapes_derivatives, ghost_type); + + const auto & connectivity = contact_map[slave].connectivity; + Matrix elementary_stiffness(connectivity.size() * spatial_dimension, + connectivity.size() * spatial_dimension); + + Array * tangents = + new Array(spatial_dimension - 1, spatial_dimension, "surface_tangents"); + + Array * global_coords = + new Array(nb_nodes_master, spatial_dimension); + + computeCoordinates(master.type, *global_coords); + computeTangents(shapes_derivatives, *global_coords, *tangents); + + Matrix surface_matrix(spatial_dimension - 1, spatial_dimension - 1); + computeSurfaceMatrix(*tangents, surface_matrix); + + computeN( *n, shapes, normal); + computeTalpha(*t_alpha, shapes, *tangents); + computeNalpha(*n_alpha, *shapes_derivatives, normal); + computeDalpha(*d_alpha, *n_alpha, *t_alpha, surface_matrix); + + computeTangentModuli(*n, *n_alpha, *t_alpha, *d_alpha, gap); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computeTangents(Matrix & shapes_derivatives, Array & global_coords, + Array & tangents) { + + UInt i = 0; + for (auto && values : zip(make_view(tangents, spatial_dimension))) { + auto & tangent = std::get<0>(values); + for (UInt n : arange(global_coords.nb_components)) { + tangent += shapes_derivaties(n, i) * global_coords(n); + } + ++i; + } + +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computeSurfaceMatrix(Array & tangents, Matrix & surface_matrix) { + + for (UInt i : arange(spatial_dimension - 1)) { + for (UInt j : arange(spatial_dimension -1 )) { + surface_matrix(i, j) = tangents(i) * tangents(j); + } + } + + inverse(surface_matrix); +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computeN(Array & n, Vector & shapes, Vector & normal) { + + UInt dim = normal.size(); + for (UInt i = 0; i < dim; ++i) { + n[i] = normal[i] * tn; + for (UInt j = 0; j < shapes.size(); ++j) { + n[(1 + j) * dim + i] = -normal[i] * shapes[j]; + } + } +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computeTalpha(Array & t_alpha, Vector & shapes, + Array & tangents) { + + for (auto && values: + zip(make_view(tangents, spatial_dimension), + make_view(t_alpha, t_alpha.size()))) { + + auto & tangent = std::get<0>(values); + auto & t = std::get<1>(values); + for (UInt i : arange(spatial_dimension)) { + t[i] = -tangent[i]; + for (UInt j : arange(shapes.size())) { + t[(1 + j)*spatial_dimension + i] = -shapes[j] * tangent[i]; + } + } + } +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computeNalpha(Array & n_alpha, Array & shapes_derivatives, + Vector & normal) { + + for (auto && values: + zip(make_view(shapes_derivatives, shapes_derivatives.size(), + n_alpha, n_alpha.size()))) { + + auto & shape_derivative = std::get<0>(values); + auto & n = std::get<1>(values); + for (UInt i : arange(spatial_dimension)) { + n[i] = 0; + for (UInt j : arange(shapes.size())) { + n[(1 + j)*spatial_dimension + i] = -shape_derivative[j]*normal[i]; + } + } + } +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computeDalpha(Array & d_alpha, Array & n_alpha, + Array & t_alpha, Matrix & surface_matrix, Real gap) { + + +} + + +/* -------------------------------------------------------------------------- */ +void Resolution::computeCoordiantes(const Element & el, Array & coords) { + UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); + Vector connect = model.getMesh().getConnectivity(el.type, _not_ghost) + .begin(nb_nodes_per_element)[el.element]; + + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt node = connect[n]; + for (UInt s: arange(spatial_dimension)) { + coords(s, n) = this->positions(node, s); + } + } +} + diff --git a/src/model/contact_mechanics/resolution.hh b/src/model/contact_mechanics/resolution.hh new file mode 100644 index 000000000..5792a1e6e --- /dev/null +++ b/src/model/contact_mechanics/resolution.hh @@ -0,0 +1,195 @@ +/** + * @file resolution.hh + * + * @author Mohit Pundir + * + * @date creation: Mon Jan 7 2019 + * @date last modification: Mon Jan 7 2019 + * + * @brief Mother class for all contact resolutions + * + * @section LICENSE + * + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "aka_factory.hh" +#include "aka_memory.hh" +#include "parsable.hh" +#include "parser.hh" +#include "fe_engine.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_RESOLUTION_HH__ +#define __AKANTU_RESOLUTION_HH__ + +/* -------------------------------------------------------------------------- */ +namespace akantu { + class Model; + class ContactMechanicsModel; +} // namespace akantu + + +namespace akantu { + +/** + * Interface of all contact resolutions + * Prerequisites for a new resolution + * - inherit from this class + * - implement the following methods: + * \code + * + * virtual void computeNormalForce(); + * virtual void computeFrictionForce(); + * + * \endcode + * + */ +class Resolution : public Memory, + public Parsable { + + /* ------------------------------------------------------------------------ */ + /* Constructor/Destructor */ + /* ------------------------------------------------------------------------ */ +public: + /// instantiate contact resolution with defaults + Resolution(ContactMechanicsModel & model, const ID & id = ""); + + /// Destructor + ~Resolution() override; + +protected: + void initialize(); + + /* ------------------------------------------------------------------------ */ + /* Functions that resolutions can/should reimplement */ + /* ------------------------------------------------------------------------ */ +protected: + /// computes the normal force + virtual void computeNormalForce(__attribute__((unused)) Vector & force, + __attribute__((unused)) Array & n, + __attribute__((unused)) Real & gap) { + AKANTU_TO_IMPLEMENT(); + } + + /// computes the friction force + virtual void computeFrictionForce(__attribute__((unused)) Vector & force, + __attribute__((unused)) Array & d_alpha, + __attribute__((unused)) Real & gap) { + AKANTU_TO_IMPLEMENT(); + } + + /// compute the tangent moduli + virtual void computeTangentModuli(__attribute__((unused)) Array & n, + __attribute__((unused)) Array & n_alpha, + __attribute__((unused)) Array & t_alpha, + __attribute__((unused)) Array & d_alpha, + __attribute__((unused)) Real & gap + ) { + AKANTU_TO_IMPLEMENT(); + } + + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + /// assemble the residual for this resolution + void assembleInternalForces(GhostType ghost_type); + + /// assemble the stiffness matrix for this resolution + void assembleStiffnessMatrix(GhostType ghost_type); + +public: + /// function to print the contain of the class + void printself(std::ostream & stream, int indent = 0) const override; + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +protected: + /// Link to the fem object in the model + FEEngine & fem; + + /// resolution name + std::string name; + + /// model to which the resolution belong + ContactMechanicsModel & model; + + /// friciton coefficient : mu + Real mu; + + /// spatial dimension + UInt spatial_dimension; + +}; + +/// standard output stream operator +inline std::ostream & operator<<(std::ostream & stream, + const Resolution & _this) { + _this.printself(stream); + return stream; +} + + +} // namespace akantu + +/* -------------------------------------------------------------------------- */ +namespace akantu { +using ResolutionFactory = + Factory; + +template +T macaulay(T var) {return var < 0 ? 0 : var; } + +template +T heaviside(T var) {return var < 0 ? 0 : 1; } +} // namespace akantu + +#define INSTANTIATE_RESOLUTION_ONLY(res_name) \ + class res_name + +#define RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name) \ + [](UInt dim, const ID &, ContactMechanicsModel & model, \ + const ID & id) -> std::unique_ptr { \ + switch (dim) { \ + case 1: \ + return std::make_unique(model, id); \ + case 2: \ + return std::make_unique(model, id); \ + case 3: \ + return std::make_unique(model, id); \ + default: \ + AKANTU_EXCEPTION("The dimension " \ + << dim << "is not a valid dimension for the contact resolution " \ + << #id); \ + } \ + } + + +#define INSTANTIATE_RESOLUTION(id, res_name) \ + INSTANTIATE_RESOLUTION_ONLY(res_name); \ + static bool resolution_is_alocated_##id[[gnu::unused]] = \ + ResolutionFactory::getInstance().registerAllocator( \ + #id, RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name)) + +#endif /* __AKANTU_RESOLUTION_HH__ */ + + diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc new file mode 100644 index 000000000..72f7fe0fb --- /dev/null +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -0,0 +1,106 @@ +/** + * @file resolution_penalty.cc + * + * @author Mohit Pundir + * + * @date creation: Mon Jan 7 2019 + * @date last modification: Mon Jan 7 2019 + * + * @brief Specialization of the resolution class for the penalty method + * + * @section LICENSE + * + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + + +/* -------------------------------------------------------------------------- */ +#include "resolution_penalty.hh" + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model, + const ID & id) + : Resolution(model, id) { + AKANTU_DEBUG_IN(); + this->initialize(); + AKANTU_DEBUG_OUT(); +} + + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::initialize() { + this->registerParam("epsilon", epsilon, Real(0.), _pat_parsable | _pat_modifiable, + "Penalty parameter"); +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeNormalForce(Vector & force, Array & n, + Real & gap) { + Real tn = gap * epsilon; + tn = macaulay(tn); + + for (UInt i : arange(force.size())) { + force[i] += tn * n[i]; + } +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeFrictionForce(Vector & force, Array & d_alpha, + Real & gap) { + ///needs to calculate tt + + UInt alpha = d_alpha.getNbCompoment(); + for (auto && values: + zip(make_view(d_alpha, d_alpha.size()))) { + auto & d = std::get<0>(values); + for (UInt i : arange(d.size())) { + force[i] += d[i] * tt; + } + } + +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeTangentModuli(Array & n, Array & n_alpha, + Array & t_alpha, Array & d_alpha, + Matrix & surface_matrix, Real & gap) { + computeNormalStiffness(n, n_alpha, d_alpha, gap); + computeFrictionalStiffness(n, n_alpha, d_alpha, gap); +} + + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeNormalStiffness(Array & n, + Array & n_alpha, Array & d_alpha, + Real & gap) { + +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeFrictionalStiffness(Array & n, + Array & n_alpha, Array & d_alpha, + Real & gap) { + +} + + +INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty); + +} // akantu diff --git a/src/model/contact_mechanics/contact_element.hh b/src/model/contact_mechanics/resolutions/resolution_penalty.hh similarity index 50% copy from src/model/contact_mechanics/contact_element.hh copy to src/model/contact_mechanics/resolutions/resolution_penalty.hh index 2b1b34a55..b246ba19e 100644 --- a/src/model/contact_mechanics/contact_element.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh @@ -1,113 +1,97 @@ /** - * @file contact_element.hh + * @file contact_resolution_penalty.hh * * @author Mohit Pundir * - * @date creation: Wed Sep 12 2018 - * @date last modification: Tue Oct 23 2018 + * @date creation: Fri Jun 18 2010 + * @date last modification: Mon Jan 29 2018 * - * @brief Mother class for all detection algorithms + * @brief Material isotropic thermo-elastic * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" - - +#include "resolution.hh" /* -------------------------------------------------------------------------- */ -#ifndef __AKANTU_CONTACT_ELEMENT_HH__ -#define __AKANTU_CONTACT_ELEMENT_HH__ - -/* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_RESOLUTION_PENALTY_HH__ +#define __AKANTU_RESOLUTION_PENALTY_HH__ namespace akantu { - -using SlaveType = UInt; -using MasterType = Element; - -class ContactElement { - +class ResolutionPenalty : public Resolution { /* ------------------------------------------------------------------------ */ - /* Constructor/ Destructors */ + /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: + ResolutionPenalty(ContactMechanicsModel & model, const ID & id = ""); + ~ResolutionPenalty() override = default; - ContactElement(SlaveType slave, MasterType master) - : node(slave), master(master) { - - } +protected: + /// initialize the resolution + void initialize(); + + /// local computation of stifnness matrix due to normal stress + void computeNormalStiffness(Array & n, + Array & n_alpha, Array & d_alpha, + Real & gap); + + /// local computation of stiffness matrix due to frictional stress + void computeFricitonalStiffness(Array & n, + Array & n_alpha, Array & d_alpha, + Real & gap); - ~ContactElement() = default; - /* ------------------------------------------------------------------------ */ - /* Accessors */ + /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// sets the value of normal vector - AKANTU_SET_MACRO(Normal, normal, Vector); - - /// sets the value of tangent vector - AKANTU_SET_MACRO(Tangent, tangent, Vector); - - /// sets the value of gap - AKANTU_SET_MACRO(Gap, gap, Real); - - /// gets the value of gap - AKANTU_GET_MACRO(Gap, gap, Real); - - /// gets the normal to the master element - AKANTU_GET_MACRO(Normal, normal, Vector); - - // gets the node of slave - AKANTU_GET_MACRO(Slave, node, SlaveType); - + /// local computation of tangent moduli + void computeTangentModuli(Array & /* N */, + Array & /* N_alpha */, + Array & /* T_alpha */, + Array & /* D_alpha */, + Real & /* gap */ + ) override; + + /// local computation of normal force + void computeNormalForce(Vector & /* force vector */, + Array & /* n */, + Real & /* gap */) override; + + /// local computation of friction force + void computeFrictionForce(Vector & /* force vector */, + Array & /* D_alpha */, + Real & /* gap */) override; - // sets the value of normal vector - AKANTU_SET_MACRO(Patch, patch, Array); - /* ------------------------------------------------------------------------ */ /* Class Members */ - /* ------------------------------------------------------------------------ */ -public: - - /// slave node - SlaveType node; + /* ------------------------------------------------------------------------ */ +protected: + /// penalty parameter + Real epsilon; - /// master element/node - MasterType master; - - /// normalized normal direction - Vector normal; - - /// normalized tangent direction - Vector tangent; - - /// penetration gap between slave and master - Real gap; - - /// an array of master nodes/elements around slave node - Array patch; }; + } // akantu -#endif /* __AKANTU_CONTACT_ELEMENT_HH__ */ + +#endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */ diff --git a/src/model/model_options.hh b/src/model/model_options.hh index a50c85055..6bd681587 100644 --- a/src/model/model_options.hh +++ b/src/model/model_options.hh @@ -1,162 +1,162 @@ /** * @file model_options.hh * * @author Nicolas Richart * * @date creation: Mon Dec 04 2017 * @date last modification: Wed Jan 31 2018 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_named_argument.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_OPTIONS_HH__ #define __AKANTU_MODEL_OPTIONS_HH__ namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(analysis_method); } struct ModelOptions { explicit ModelOptions(AnalysisMethod analysis_method = _static) : analysis_method(analysis_method) {} template ModelOptions(use_named_args_t, pack &&... _pack) : ModelOptions(OPTIONAL_NAMED_ARG(analysis_method, _static)) {} virtual ~ModelOptions() = default; AnalysisMethod analysis_method; }; #ifdef AKANTU_SOLID_MECHANICS /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelOptions : public ModelOptions { explicit SolidMechanicsModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass) : ModelOptions(analysis_method) {} template SolidMechanicsModelOptions(use_named_args_t, pack &&... _pack) : SolidMechanicsModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} }; #endif /* -------------------------------------------------------------------------- */ #ifdef AKANTU_COHESIVE_ELEMENT namespace { DECLARE_NAMED_ARGUMENT(is_extrinsic); } /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelCohesiveOptions : public SolidMechanicsModelOptions { SolidMechanicsModelCohesiveOptions( AnalysisMethod analysis_method = _explicit_lumped_mass, bool extrinsic = false) : SolidMechanicsModelOptions(analysis_method), is_extrinsic(extrinsic) {} template SolidMechanicsModelCohesiveOptions(use_named_args_t, pack &&... _pack) : SolidMechanicsModelCohesiveOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass), OPTIONAL_NAMED_ARG(is_extrinsic, false)) {} bool is_extrinsic{false}; }; #endif #ifdef AKANTU_HEAT_TRANSFER /* -------------------------------------------------------------------------- */ struct HeatTransferModelOptions : public ModelOptions { explicit HeatTransferModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass) : ModelOptions(analysis_method) {} template HeatTransferModelOptions(use_named_args_t, pack &&... _pack) : HeatTransferModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} }; #endif #ifdef AKANTU_EMBEDDED namespace { DECLARE_NAMED_ARGUMENT(init_intersections); } /* -------------------------------------------------------------------------- */ struct EmbeddedInterfaceModelOptions : SolidMechanicsModelOptions { /** * @brief Constructor for EmbeddedInterfaceModelOptions * @param analysis_method see SolidMechanicsModelOptions * @param init_intersections compute intersections */ EmbeddedInterfaceModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass, bool init_intersections = true) : SolidMechanicsModelOptions(analysis_method), has_intersections(init_intersections) {} template EmbeddedInterfaceModelOptions(use_named_args_t, pack &&... _pack) : EmbeddedInterfaceModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass), OPTIONAL_NAMED_ARG(init_intersections, true)) {} /// Should consider reinforcements bool has_intersections; }; #endif #ifdef AKANTU_CONTACT_MECHANICS namespace { DECLARE_NAMED_ARGUMENT(is_explicit); } /* -------------------------------------------------------------------------- */ struct ContactMechanicsModelOptions : public ModelOptions { explicit ContactMechanicsModelOptions( - AnalysisMethod analysis_method = _explicit_contact_penalty, + AnalysisMethod analysis_method = _explicit_contact, bool explicit = true) : ModelOptions(analysis_method), is_explicit(explicit) {} template ContactMechanicsModelOptions(use_named_args_t, pack &&... _pack) : ContactMechanicsModelOptions( - OPTIONAL_NAMED_ARG(analysis_method, _explicit_contact_penalty), + OPTIONAL_NAMED_ARG(analysis_method, _explicit_contact), OPTIONAL_NAMED_ARG(is_explicit, true)) {} bool is_explicit{true}; }; #endif } // akantu #endif /* __AKANTU_MODEL_OPTIONS_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 0edab13c7..4915be9bf 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,568 +1,568 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); - /// apply a constant eigen_grad_u on all quadrature points of a given material + /// applya constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the SolidMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); AKANTU_SET_MACRO(MaterialSelector, material_selector, std::shared_ptr); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// lumped mass array Array * mass{nullptr}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force{nullptr}; /// internal forces array Array * internal_force{nullptr}; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs{nullptr}; /// array of current position used during update residual Array * current_position{nullptr}; UInt current_position_release{0}; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; /// class defining of to choose a material std::shared_ptr material_selector; /// flag defining if the increment must be computed or not bool increment_flag; /// tells if the material are instantiated bool are_materials_instantiated; using flatten_internal_map = std::map, ElementTypeMapArray *>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/test/test_model/test_contact_mechanics_model/CMakeLists.txt b/test/test_model/test_contact_mechanics_model/CMakeLists.txt index 5500401e5..f03aa5a6f 100644 --- a/test/test_model/test_contact_mechanics_model/CMakeLists.txt +++ b/test/test_model/test_contact_mechanics_model/CMakeLists.txt @@ -1,25 +1,26 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Tue Dec 18 2018 # @date last modification: Tue Dec 18 2018 # # @brief configuration for ContactMechanicsModel tests # # @section LICENSE # # Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_akantu_test(test_detector "detection_test") +add_akantu_test(test_resolution "resolution_test") diff --git a/test/test_model/test_contact_mechanics_model/CMakeLists.txt b/test/test_model/test_contact_mechanics_model/test_resolution/CMakeLists.txt similarity index 76% copy from test/test_model/test_contact_mechanics_model/CMakeLists.txt copy to test/test_model/test_contact_mechanics_model/test_resolution/CMakeLists.txt index 5500401e5..3334ff98e 100644 --- a/test/test_model/test_contact_mechanics_model/CMakeLists.txt +++ b/test/test_model/test_contact_mechanics_model/test_resolution/CMakeLists.txt @@ -1,25 +1,32 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # -# @date creation: Tue Dec 18 2018 -# @date last modification: Tue Dec 18 2018 +# @date creation: Mon Jan 7 2019 +# @date last modification: Mon Jan 7 2019 # -# @brief configuration for ContactMechanicsModel tests +# @brief configuration for contact resolution tests # # @section LICENSE # # Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== -add_akantu_test(test_detector "detection_test") +add_mesh(explicit_2d explicit_2d.geo 2 1) + +register_test(test_explicit_resolution + SOURCES test_explicit_resolution.cc + DEPENDS explicit_2d + FILES_TO_COPY options.dat + PACKAGE contact_mechanics + ) diff --git a/test/test_model/test_contact_mechanics_model/test_resolution/explicit_2d.geo b/test/test_model/test_contact_mechanics_model/test_resolution/explicit_2d.geo new file mode 100644 index 000000000..4bd8b36d5 --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_resolution/explicit_2d.geo @@ -0,0 +1,29 @@ +cl1 = 0.05; +cl2 = 0.05; +cl3 = 0.05; +Dy = 0.05; +Point(1) = {0, 0.1-Dy, 0, cl1}; +Point(2) = {0.5, 0.6-Dy, 0, cl2}; +Point(3) = {-0.5, 0.6-Dy, 0, cl2}; +Point(4) = {0.5, 0.1, 0, cl3}; +Point(5) = {-0.5, 0.1, 0, cl3}; +Point(6) = {-0.5, -0.25, 0, cl3}; +Point(7) = {0.5, -0.25, 0, cl3}; +Point(8) = {0, 0.6-Dy, 0, cl2}; +Circle(1) = {3, 8, 1}; +Circle(2) = {1, 8, 2}; +Line(3) = {2, 8}; +Line(13) = {8, 3}; +Line(4) = {6, 7}; +Line(5) = {7, 4}; +Line(6) = {4, 5}; +Line(7) = {5, 6}; +Line Loop(9) = {2, 3, 13, 1}; +Plane Surface(9) = {9}; +Line Loop(11) = {6, 7, 4, 5}; +Plane Surface(11) = {11}; +Physical Line("contact_surface") = {6}; +Physical Line("rigid") = {1, 2}; +Physical Line("top") = {3, 13}; +Physical Surface("top_body") = {9}; +Physical Surface("bot_body") = {11}; \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_resolution/options.dat b/test/test_model/test_contact_mechanics_model/test_resolution/options.dat new file mode 100644 index 000000000..93c6fa956 --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_resolution/options.dat @@ -0,0 +1,4 @@ +contact_resolution_penalty [ + name = dummy + epsilon = 10.0 +] \ No newline at end of file diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/test/test_model/test_contact_mechanics_model/test_resolution/test_explicit_resolution.cc similarity index 58% copy from src/model/contact_mechanics/contact_mechanics_model.cc copy to test/test_model/test_contact_mechanics_model/test_resolution/test_explicit_resolution.cc index 29aaa564d..26c9c647c 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/test/test_model/test_contact_mechanics_model/test_resolution/test_explicit_resolution.cc @@ -1,68 +1,53 @@ /** - * @file coontact_mechanics_model.cc + * @file test_resolution_explicit.cc * * @author Mohit Pundir * - * @date creation: Tue May 08 2012 - * @date last modification: Wed Feb 21 2018 + * @date creation: Mon Jan 7 2019 + * @date last modification: Mon Jan 7 2019 * - * @brief Contact mechanics model + * @brief Test for explicit contact resolution * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ -#include "fe_engine.hh" #include "contact_mechanics_model.hh" /* -------------------------------------------------------------------------- */ -#include -/* -------------------------------------------------------------------------- */ - - -namespace akantu { +using namespace akantu; -ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, - const MemoryID & memory_id, - const ModelType model_type) - : Model(mesh, model_type, dim, id, memory_id) { +const UInt spatial_dimension = 2; - AKANTU_DEBUG_IN(); +int main(int argc, char *argv[]) +{ + initialize("options.dat", argc, argv); - - //this->detector = std::make_unique( - // this->mesh, id + ":contact_detector"); - - AKANTU_DEBUG_OUT(); + Mesh mesh(spatial_dimension); + //mesh.read("explicit_2d.msh"); -} + ContactMechanicsModel model(mesh); + model.initFull(_analysis_method = _static); - -ContactMechanicsModel::~ContactMechanicsModel() { - AKANTU_DEBUG_IN(); - - AKANTU_DEBUG_OUT(); -} - -void ContactMechanicsModel::initModel() { + std::cout << model; + finalize(); + return EXIT_SUCCESS; } - -}