diff --git a/examples/contact_mechanics/CMakeLists.txt b/examples/contact_mechanics/CMakeLists.txt index ae6415a11..be4a64e24 100644 --- a/examples/contact_mechanics/CMakeLists.txt +++ b/examples/contact_mechanics/CMakeLists.txt @@ -1,52 +1,60 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Mon Jan 18 2016 # # @brief configuration for heat transfer example # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_mesh(hertz_2d hertz_2d.geo 2 1) +add_mesh(hertz_implicit_2d hertz_implicit_2d.geo 2 1) add_mesh(hertz_3d hertz_3d.geo 3 1) +register_example(implicit_penalty_2d + SOURCES implicit_penalty_2d.cc + DEPENDS hertz_implicit_2d + FILES_TO_COPY material_implicit.dat + ) + + register_example(explicit_penalty_2d SOURCES explicit_penalty_2d.cc DEPENDS hertz_2d FILES_TO_COPY material.dat ) register_example(explicit_penalty_3d SOURCES explicit_penalty_3d.cc DEPENDS hertz_3d FILES_TO_COPY material.dat ) register_example(solid_contact_explicit SOURCES solid_contact_explicit.cc DEPENDS hertz_2d FILES_TO_COPY material.dat ) diff --git a/examples/contact_mechanics/hertz_implicit_2d.geo b/examples/contact_mechanics/hertz_implicit_2d.geo new file mode 100644 index 000000000..d123ffd31 --- /dev/null +++ b/examples/contact_mechanics/hertz_implicit_2d.geo @@ -0,0 +1,32 @@ +cl1 = 0.01; +cl2 = 0.01; +cl3 = 0.005; +Dy = 0.0; +radius = 0.1; +y = 0.1; +Point(1) = {0, y, 0, cl1}; +Point(2) = {radius, radius + y, 0, cl2}; +Point(3) = {-radius, radius + y, 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, cl2}; +Point(7) = {0.5, -0.25, 0, cl2}; +Point(8) = {0, radius + y, 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 Line("bottom") = {4}; +Physical Surface("top_body") = {9}; +Physical Surface("bot_body") = {11}; \ No newline at end of file diff --git a/examples/contact_mechanics/implicit_penalty_2d.cc b/examples/contact_mechanics/implicit_penalty_2d.cc new file mode 100644 index 000000000..0e77b4729 --- /dev/null +++ b/examples/contact_mechanics/implicit_penalty_2d.cc @@ -0,0 +1,71 @@ +/** + * @file contact_mechanics_penalty.cc + * + * @author Mohit Pundir + * + * @date creation: Mon Jan 21 2019 + * @date last modification: Mon Jan 21 2019 + * + * @brief contact mechanics model with penalty 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 +#include + +/* -------------------------------------------------------------------------- */ +#include "contact_mechanics_model.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + + initialize("material_implicit.dat", argc, argv); + const UInt spatial_dimension = 2; + + Mesh mesh(spatial_dimension); + mesh.read("hertz_implicit_2d.msh"); + + ContactMechanicsModel model(mesh); + model.initFull(_analysis_method = _implicit_contact); + + model.setBaseName("implicit-penalty-2d"); + model.addDumpField("contact_force"); + model.addDumpField("external_force"); + model.addDumpField("blocked_dofs"); + model.addDumpField("gaps"); + model.addDumpField("areas"); + + model.dump(); + + model.search(); + model.dump(); + + model.assembleInternalForces(); + + model.dump(); + + finalize(); + return EXIT_SUCCESS; +} + diff --git a/examples/contact_mechanics/material_implicit.dat b/examples/contact_mechanics/material_implicit.dat new file mode 100644 index 000000000..80c1ccf01 --- /dev/null +++ b/examples/contact_mechanics/material_implicit.dat @@ -0,0 +1,24 @@ +material elastic [ + name = bot_body + rho = 7800 # density + E = 210 # young's modulu + nu = 0.3 # poisson's ratio +] + +material elastic [ + name = top_body + rho = 7800 # density + E = 2.1e11 # young's modulu + nu = 0.3 # poisson's ratio +] + +contact_detector [ + type = implicit + master_surface = rigid + slave_surface = contact_surface +] + +contact_resolution penalty [ + name = contact_surface + epsilon = 10.0 +] \ No newline at end of file diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 4e32c68e3..321b092d6 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,676 +1,536 @@ /** * @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 { /* -------------------------------------------------------------------------- */ /* Constants */ /* -------------------------------------------------------------------------- */ namespace { __attribute__((unused)) constexpr UInt _all_dimensions{ std::numeric_limits::max()}; #ifdef AKANTU_NDEBUG __attribute__((unused)) constexpr Real REAL_INIT_VALUE{0.}; #else __attribute__((unused)) constexpr Real REAL_INIT_VALUE{ std::numeric_limits::quiet_NaN()}; #endif } // namespace /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; using MemoryID = UInt; } /* -------------------------------------------------------------------------- */ #ifndef SWIG #include "aka_enum_macros.hh" #endif /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of /// events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #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_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { _model, _solid_mechanics_model, _solid_mechanics_model_cohesive, _heat_transfer_model, _structural_mechanics_model, _embedded_model, _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 = 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 ContactDetectorType { _explicit, _implicit }; -/// 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)); -} - /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_ask_nodes, _gst_size, //--- SolidMechanicsModel tags --- _gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _gst_smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _gst_smm_uv, ///< synchronization of the nodal velocities and displacement _gst_smm_res, ///< synchronization of the nodal residual _gst_smm_init_mat, ///< synchronization of the data to initialize materials _gst_smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _gst_smmc_facets, ///< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, ///< synchronization of facet global connectivity _gst_smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _gst_smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _gst_giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _gst_ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gst_gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _gst_htm_temperature, ///< synchronization of the nodal temperature _gst_htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _gst_htm_phi, ///< synchronization of the nodal level set value phi _gst_htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _gst_mnl_for_average, ///< synchronization of data to average in non local /// material _gst_mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _gst_nh_criterion, // --- General tags --- _gst_test, ///< Test tag _gst_user_1, ///< tag for user simulations _gst_user_2, ///< tag for user simulations _gst_material_id, ///< synchronization of the material ids _gst_for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _gst_cf_nodal, ///< synchronization of disp, velo, and current position _gst_cf_incr, ///< synchronization of increment // --- Solver tags --- _gst_solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _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/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 841c7e506..6f5ada151 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,565 +1,564 @@ /** * @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 == "implict") { + if (type == "implicit") { this->detection_type = _implicit; } else if (type == "explicit"){ this->detection_type = _explicit; } else { AKANTU_ERROR("Unknown detection type : " << type); } surfaces[Surface::master] = section.getParameterValue("master_surface"); surfaces[Surface::slave ] = section.getParameterValue("slave_surface"); } /* -------------------------------------------------------------------------- */ 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 = 2.0 * max_el_size; + this->max_dd = max_el_size; this->max_bb = max_el_size; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactDetector::search(std::map & contact_map) { this->globalSearch(contact_map); } /* -------------------------------------------------------------------------- */ 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, contact_map); } /* -------------------------------------------------------------------------- */ void ContactDetector::localSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid, 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"); auto natural_projections = std::make_unique>(elements.size(), spatial_dimension - 1, "projections"); this->computeOrthogonalProjection(slave_node, elements, *normals, *gaps, *natural_projections); auto minimum = std::min_element( gaps->begin(), gaps->end()); auto index = std::distance( gaps->begin(), minimum); - Vector master_conn = this->mesh.getConnectivity(elements[index]); + Vector master_conn = + this->mesh.getConnectivity(elements[index]); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave_node; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i-1]; } contact_map[slave_node] = ContactElement(elements[index]); contact_map[slave_node].gap = (*gaps)[index]; contact_map[slave_node].normal = Vector(normals->begin(spatial_dimension)[index], true); contact_map[slave_node].projection = Vector(natural_projections->begin(spatial_dimension - 1)[index], true); contact_map[slave_node].connectivity = elem_conn; } } /* -------------------------------------------------------------------------- */ 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; + 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, Array & natural_projections) { 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), make_view(natural_projections, spatial_dimension - 1))) { const auto & element = std::get<0>(values); auto & gap = std::get<1>(values); auto & normal = std::get<2>(values); auto & natural_projection = std::get<3>(values); this->computeNormalOnElement(element, normal); + Vector real_projection(spatial_dimension); this->computeProjectionOnElement(element, normal, query, natural_projection, real_projection); Vector distance(spatial_dimension); distance = query - real_projection; gap = Math::norm(spatial_dimension, distance.storage()); Vector direction = distance.normalize(); Real cos_angle = direction.dot(normal); Real tolerance = 1e-8; if (std::abs(cos_angle - 1) <= tolerance && detection_type == _explicit) { gap *= -1; } } } /* -------------------------------------------------------------------------- */ 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); } } } /* -------------------------------------------------------------------------- */ void ContactDetector::computeProjectionOnElement(const Element & element, const Vector & normal, const Vector & query, Vector & natural_projection, Vector & real_projection) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); Matrix coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(element, coords); Vector point(coords(0)); Real alpha = (query - point).dot(normal); real_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, real_projection, natural_projection); } /* -------------------------------------------------------------------------- */ bool ContactDetector::isValidProjection(const Element & element, Vector & real_projection, Vector & natural_projection) { const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(this->positions /*mesh.getNodes()*/, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); #define GET_NATURAL_COORDINATE(type) \ ElementClass::inverseMap(real_projection, nodes_coord, natural_projection) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE); #undef GET_NATURAL_COORDINATE /*Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); Real distance = 0; switch (this->spatial_dimension) { case 2: { distance = Math::distance_2d(real_projection.storage(), barycenter.storage()); break; } case 3: { distance = Math::distance_3d(real_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); } } } /* -------------------------------------------------------------------------- */ void ContactDetector::normalProjection(const Element & el, const Vector & slave_coord, Vector & natural_coord, Real & tolerance) { /*Real fmin; auto update_fmin = [&fmin, &slave_coord, &node_coords, &natural_coord]() { Vector physical_guess_v(physical_guess.storage(), spatial_dimension); // interpolate on natual coordinate and get physical guess // compute gradient or jacobian on natural cooordiante // f = slave_coord - physical_guess; return fmin; }; auto closet_point_error = update_fmin(); while (tolerance < closet_point_error) { // compute gradiend on natural coordinate // compute second variation of shape function at natural coord closet_point_error = update_fmin(); }*/ } } // akantu diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 4debf41e4..a1f87b03d 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,589 +1,588 @@ /** * @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 "contact_mechanics_model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #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), current_positions(mesh.getNodes()) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->initDOFManager(); this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, Array & positions, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), current_positions(positions) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->initDOFManager(); this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, positions, 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(); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ 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()); are_resolutions_instantiated = true; } /* -------------------------------------------------------------------------- */ 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)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->contact_force, spatial_dimension, "contact_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->areas, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->contact_force, _dst_nodal); //dof_manager.registerDOFs("displacement", *this->blocked_dofs); } } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", _tsst_dynamic); } case _implicit_contact: { return std::make_tuple("implicit_contact", _tsst_static); } default: return std::make_tuple("unkown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type ) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); // assemble the forces due to local stresses auto assemble = [&] (auto && ghost_type) { for (auto & resolution : resolutions) { resolution->assembleInternalForces(ghost_type); } }; AKANTU_DEBUG_INFO("Assemble residual for local elements"); assemble(_not_ghost); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); assemble(_ghost); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->contact_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { this->detector->search(this->contact_map); auto & blocked_dof = const_cast &>(this->getBlockedDOFs()); auto & gap = const_cast &>(this->getGaps()); for(auto & entry : contact_map) { const auto & element = entry.second; const auto & connectivity = element.connectivity; for (UInt i = 0; i < connectivity.size(); ++i) { UInt n = connectivity(i); blocked_dof[n] = 1.0; gap[n] = element.gap; } } - this->areas->clear(); + this->areas->clear(); this->external_force->clear(); this->applyBC(BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), this->detector->getSurfaceId("slave")); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto areas_it = areas->begin(); UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++areas_it) { const auto & ext_force = *ext_force_it; auto & area = *areas_it; for (UInt i = 0; i < Model::spatial_dimension; ++i) { area += pow(ext_force(i), 2); } area = sqrt(area); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() { } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep() { } /* -------------------------------------------------------------------------- */ 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 == "K") return _symmetric; return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().clearMatrix("K"); for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER - dumper::Field * ContactMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { //std::map *> uint_nodal_fields; //uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = nullptr; //field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * ContactMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["contact_force"] = this->contact_force; real_nodal_fields["external_force"] = this->external_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["areas"] = this->areas; dumper::Field * field = nullptr; if (padding_flag) field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); else field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } #else /* -------------------------------------------------------------------------- */ dumper::Field * ContactMechanicsModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index 1771b32fb..76d286d9d 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,321 +1,359 @@ /** * @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(); auto & internal_force = const_cast &>(model.getInternalForce()); auto & contact_area = const_cast &>(model.getContactArea()); auto & contact_map = model.getContactMap(); const auto slave_nodes = model.getMesh().getElementGroup(name).getNodes(); for (auto & slave: slave_nodes) { if (contact_map.find(slave) == contact_map.end()) continue; auto & master = contact_map[slave].master; auto & gap = contact_map[slave].gap; auto & projection = contact_map[slave].projection; auto & normal = contact_map[slave].normal; const auto & connectivity = contact_map[slave].connectivity; const ElementType & type = master.type; UInt nb_nodes_master = Mesh::getNbNodesPerElement(master.type); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shapes_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Vector elem_force(connectivity.size() * spatial_dimension); Matrix tangents(spatial_dimension - 1, spatial_dimension); Matrix global_coords(nb_nodes_master, spatial_dimension); computeCoordinates(master, global_coords); computeTangents(shapes_derivatives, global_coords, tangents); Matrix surface_matrix(spatial_dimension - 1, spatial_dimension - 1); computeSurfaceMatrix(tangents, surface_matrix); - Array n(connectivity.size() * spatial_dimension, 1, "n_array"); + Vector n(connectivity.size() * spatial_dimension); computeN(n, shapes, normal); computeNormalForce(elem_force, n, gap); - /*Array t_alpha(connectivity.size() * spatial_dimension, "t_alpha_array"); - Array n_alpha(connectivity.size() * spatial_dimension, "n_alpha_array"); - Array d_alpha(connectivity.size() * spatial_dimension, "d_alpha_array"); - + Array t_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); + Array n_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); + Array d_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); + computeTalpha(t_alpha, shapes, tangents); computeNalpha(n_alpha, shapes_derivatives, normal); computeDalpha(d_alpha, n_alpha, t_alpha, surface_matrix, gap); - computeFrictionForce(elem_force, d_alpha, gap);*/ + + //computeFrictionForce(elem_force, d_alpha, gap); UInt nb_degree_of_freedom = internal_force.getNbComponent(); for (UInt i = 0; i < connectivity.size(); ++i) { UInt n = connectivity[i]; for (UInt j = 0; j < nb_degree_of_freedom; ++j) { UInt offset_node = n * nb_degree_of_freedom + j; internal_force[offset_node] += elem_force[i*nb_degree_of_freedom + j]; internal_force[offset_node] *= contact_area[n]; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); - const auto slave_nodes = model.getMesh().getElementGroup(name).getNodes(); + auto & contact_stiffness = + model.getDOFManager().getMatrix("K"); + + const auto slave_nodes = + model.getMesh().getElementGroup(name).getNodes(); + auto & contact_map = model.getContactMap(); for (auto & slave: slave_nodes) { if (contact_map.find(slave) == contact_map.end()) { continue; } auto & master = contact_map[slave].master; auto & gap = contact_map[slave].gap; auto & projection = contact_map[slave].projection; auto & normal = contact_map[slave].normal; const auto & connectivity = contact_map[slave].connectivity; const ElementType & type = master.type; UInt nb_nodes_master = Mesh::getNbNodesPerElement(master.type); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shapes_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Matrix elementary_stiffness(connectivity.size() * spatial_dimension, connectivity.size() * spatial_dimension); Matrix tangents(spatial_dimension - 1, spatial_dimension); Matrix global_coords(nb_nodes_master, spatial_dimension); computeCoordinates(master, global_coords); computeTangents(shapes_derivatives, global_coords, tangents); Matrix surface_matrix(spatial_dimension - 1, spatial_dimension - 1); computeSurfaceMatrix(tangents, surface_matrix); - Array n(connectivity.size() * spatial_dimension, 1, "n_array"); - /*Array t_alpha(connectivity.size() * spatial_dimension, spatial_dimension -1, "t_alpha_array"); - Array n_alpha(connectivity.size() * spatial_dimension, spatial_dimension -1, "n_alpha_array"); - Array d_alpha(connectivity.size() * spatial_dimension, spatial_dimension -1, "d_alpha_array"); + Vector n(connectivity.size() * spatial_dimension); + Array t_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); + Array n_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); + Array d_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); computeN( n, shapes, normal); computeTalpha( t_alpha, shapes, tangents); computeNalpha( n_alpha, shapes_derivatives, normal); - computeDalpha( d_alpha, n_alpha, t_alpha, surface_matrix, gap);*/ + computeDalpha( d_alpha, n_alpha, t_alpha, surface_matrix, gap); //computeTangentModuli(n, n_alpha, t_alpha, d_alpha, gap); + + /*std::vector equations; + UInt nb_degree_of_freedom = Model::spatial_dimension; + for (UInt i : arange(connectivity.size())) { + UInt n = connectivity[i]; + for (UInt j : arange(nb_degree_of_freedom)) + equations.push_back(n * degree_of_freedom + j); + } + + for (UInt i : arange(kc.rows())) { + UInt row = equations[i]; + for (UInt j : arange(kc.cols())) { + UInt col = equations[j]; + contact_stiffness(row, col) += kc(i, j); + } + }*/ } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::computeTangents(Matrix & shapes_derivatives, Matrix & global_coords, Matrix & tangents) { /*UInt index = 0; for (auto && values : zip(make_view(tangents, spatial_dimension))) { auto & tangent = std::get<0>(values); for (UInt n : arange(global_coords.getNbComponent())) { tangent += shapes_derivatives(n, index) * global_coords(n); } ++index; }*/ tangents.mul(shapes_derivatives, global_coords); } /* -------------------------------------------------------------------------- */ void Resolution::computeSurfaceMatrix(Matrix & tangents, Matrix & surface_matrix) { /*Matrix A(surface_matrix); for (UInt i : arange(spatial_dimension - 1)) { for (UInt j : arange(spatial_dimension -1 )) { A(i, j) = tangents(i) * tangents(j); } }*/ surface_matrix.mul(tangents, tangents); surface_matrix = surface_matrix.inverse(); //A.inverse(surface_matrix); } /* -------------------------------------------------------------------------- */ -void Resolution::computeN(Array & n, Vector & shapes, Vector & normal) { +void Resolution::computeN(Vector & n, Vector & shapes, Vector & normal) { UInt dim = normal.size(); for (UInt i = 0; i < dim; ++i) { n[i] = normal[i]; 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) { + Matrix & tangents) { - for (auto && values: - zip(make_view(tangents, spatial_dimension), - make_view(t_alpha, t_alpha.size()))) { + /*for (auto && values: + zip(tangents.transpose(), + make_view(t_alpha, t_alpha.size()))) { auto & tangent = std::get<0>(values); - auto & t = std::get<1>(values); + auto & t_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { - t[i] = -tangent[i]; + t_s[i] = -tangent(i); for (UInt j : arange(shapes.size())) { - t[(1 + j)*spatial_dimension + i] = -shapes[j] * tangent[i]; + t_s[(1 + j)*spatial_dimension + i] = -shapes[j] * tangent(i); } } - } + }*/ } /* -------------------------------------------------------------------------- */ void Resolution::computeNalpha(Array & n_alpha, Matrix & shapes_derivatives, Vector & normal) { - for (auto && values: - zip(make_view(n_alpha, n_alpha.size()))) { - - auto & n = std::get<0>(values); + /*for (auto && values: + zip(shapes_derivatives.transpose(), + make_view(n_alpha, n_alpha.size()))) { + auto & dnds = std::get<0>(values); + auto & n_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { - n[i] = 0; + n_s[i] = 0; for (UInt j : arange(shapes_derivatives.size())) { - n[(1 + j)*spatial_dimension + i] = -shapes_derivatives[j]*normal[i]; + n_s[(1 + j)*spatial_dimension + i] = -shapes_derivatives[j]*normal[i]; } } - } + }*/ } /* -------------------------------------------------------------------------- */ void Resolution::computeDalpha(Array & d_alpha, Array & n_alpha, - Array & t_alpha, Matrix & surface_matrix, Real & gap) { - - + Array & t_alpha, Matrix & surface_matrix, + Real & gap) { + + /*for (auto && entry : zip(surface_matrix.transpose(), + make_view(d_alpha, d_alpha.size()))) { + auto & a_s = std::get<0>(entry); + auto & d_s = std::get<1>(entry); + for (auto && values : + enumerate(make_view(t_alpha, t_alpha.size()), + make_view(n_alpha, n_alpha.size()))) { + auto & index = std::get<0>(values); + auto & t_s = std::get<1>(values); + auto & n_s = std::get<2>(values); + + d_s += (t_s + gap * n_s); + d_s *= a_s(index); + } + }*/ } /* -------------------------------------------------------------------------- */ void Resolution::computeCoordinates(const Element & el, Matrix & 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]; // change this to current position auto & positions = model.getMesh().getNodes(); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; for (UInt s: arange(spatial_dimension)) { coords(n, s) = positions(node, s); } } } } // akantu diff --git a/src/model/contact_mechanics/resolution.hh b/src/model/contact_mechanics/resolution.hh index 0545ae441..5830b84ed 100644 --- a/src/model/contact_mechanics/resolution.hh +++ b/src/model/contact_mechanics/resolution.hh @@ -1,227 +1,229 @@ /** * @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(); /// computes coordinates of a given element void computeCoordinates(const Element & , Matrix &); - /// - void computeTalpha(Array & , Vector & , Array & ); - /// computes tangents void computeTangents(Matrix & /* shapes_derivatives */, Matrix & /* global_coords */, Matrix & /* tangents */); /// computes surface metric matrix void computeSurfaceMatrix(Matrix & /* tangents */, Matrix & /* surface_matrix */); /// computes N array - void computeN(Array & /* n */, + void computeN(Vector & /* n */, Vector & /* shapes */, Vector & /* normal */); /// computes N_{\alpha} where \alpha is number of surface dimensions void computeNalpha(Array & /* n_alpha */, Matrix & /* shapes_derivatives */, Vector & /* normal */); + /// computes T_{\alpha} where \alpha is surface dimenion + void computeTalpha(Array & , + Vector & , + Matrix & ); + /// computes D_{\alpha} where \alpha is number of surface dimensions void computeDalpha(Array & /* d_alpha */, Array & /* n_alpha */, Array & /* t_alpha */, Matrix & /* surface_matrix */, Real & /* gap */); /* ------------------------------------------------------------------------ */ /* Functions that resolutions can/should reimplement */ /* ------------------------------------------------------------------------ */ protected: /// computes the normal force virtual void computeNormalForce(__attribute__((unused)) Vector & /* force */, - __attribute__((unused)) Array & /* n */, + __attribute__((unused)) Vector & /* 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)) Array & /* d_alpha */, __attribute__((unused)) Real & /* gap */) { AKANTU_TO_IMPLEMENT(); } /// compute the tangent moduli - virtual void computeTangentModuli(__attribute__((unused)) Array & /* n */, + virtual void computeTangentModuli(__attribute__((unused)) Vector & /* n */, __attribute__((unused)) Array & /* n_alpha */, __attribute__((unused)) Array & /* t_alpha */, __attribute__((unused)) Array & /* d_alpha */, __attribute__((unused)) Matrix & /* A */, __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 index b7eb11c3a..be8932500 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,105 +1,159 @@ /** * @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, +void ResolutionPenalty::computeNormalForce(Vector & force, Vector & 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 - Real tt = 0; + Vector tractions(d_alpha.getNbComponent()); + computeFrictionalTraction(tractions); + 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; - } + zip(tractions, + make_view(d_alpha, d_alpha.size()))) { + auto & t_s = std::get<0>(values); + auto & d_s = std::get<1>(values); + force += d_s * t_s; } - } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentModuli(Array & n, Array & n_alpha, +void ResolutionPenalty::computeFrictionalTraction(Vector & tractions) { + +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeTangentModuli(Vector & n, Array & n_alpha, Array & t_alpha, Array & d_alpha, Matrix & surface_matrix, Real & gap) { - computeNormalStiffness(n, n_alpha, d_alpha, gap); + computeNormalStiffness(n, n_alpha, d_alpha, surface_matrix, gap); computeFrictionalStiffness(n, n_alpha, d_alpha, gap); } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalStiffness(Array & n, - Array & n_alpha, Array & d_alpha, - Real & gap) { - +void ResolutionPenalty::computeNormalStiffness(Vector & n, + Array & n_alpha, Array & d_alpha, + Matrix & surface_matrix, Real & gap) { + + Real tn = gap * epsilon; + tn = macaulay(tn); + + Matrix ke(n_alpha.size(), n_alpha.size()); + Matrix n_mat(n.storage(), n.size(), 1); + + ke.mul(n_mat, n_mat); + ke *= epsilon * heaviside(gap); + + for (auto && entry1 : + enumerate(make_view(n_alpha, n_alpha.size()))) { + + auto & i = std::get<0>(entry1); + auto & ni = std::get<1>(entry1); + Matrix ni_mat(ni.storage(), ni.size(), 1); + + for (auto && entry2: + enumerate(make_view(n_alpha, n_alpha.size()))) { + + auto & j = std::get<0>(entry2); + auto & nj = std::get<1>(entry2); + Matrix nj_mat(nj.storage(), nj.size(), 1); + + Matrix tmp(nj.size(), nj.size()); + tmp.mul(ni_mat, nj_mat); + tmp *= surface_matrix[i, j]; + + ke += tmp; + } + } + + for (auto && values: + zip(make_view(n_alpha, n_alpha.size()), + make_view(d_alpha, d_alpha.size()))) { + auto & n_s = std::get<0>(values); + auto & d_s = std::get<1>(values); + + Matrix ns_mat(n_s.storage(), n_s.size(), 1); + Matrix ds_mat(d_s.storage(), d_s.size(), 1); + + Matrix tmp1(n_s.size(), n_s.size()); + tmp1.mul(ns_mat, ds_mat); + + Matrix tmp2(n_s.size(), n_s.size()); + tmp1.mul(ds_mat, ns_mat); + + ke -= tmp1 + tmp2; + } + } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeFrictionalStiffness(Array & n, +void ResolutionPenalty::computeFrictionalStiffness(Vector & n, Array & n_alpha, Array & d_alpha, Real & gap) { } INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty); } // akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.hh b/src/model/contact_mechanics/resolutions/resolution_penalty.hh index 128168e8c..bfb255ca2 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh @@ -1,98 +1,104 @@ /** * @file contact_resolution_penalty.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Jan 29 2018 * * @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_RESOLUTION_PENALTY_HH__ #define __AKANTU_RESOLUTION_PENALTY_HH__ namespace akantu { class ResolutionPenalty : public Resolution { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ResolutionPenalty(ContactMechanicsModel & model, const ID & id = ""); ~ResolutionPenalty() override = default; 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, + void computeNormalStiffness(Vector & n, + Array & n_alpha, + Array & d_alpha, + Matrix & surface_matrix, Real & gap); /// local computation of stiffness matrix due to frictional stress - void computeFrictionalStiffness(Array & n, - Array & n_alpha, Array & d_alpha, + void computeFrictionalStiffness(Vector & n, + Array & n_alpha, + Array & d_alpha, Real & gap); + + /// computation of frictional tractions + void computeFrictionalTraction(Vector & tractions); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// local computation of tangent moduli - void computeTangentModuli(Array & /* N */, + void computeTangentModuli(Vector & /* N */, Array & /* N_alpha */, Array & /* T_alpha */, Array & /* D_alpha */, Matrix & /* A matrix */, Real & /* gap */ ) override; /// local computation of normal force void computeNormalForce(Vector & /* force vector */, - Array & /* n */, + Vector & /* n */, Real & /* gap */) override; /// local computation of friction force void computeFrictionForce(Vector & /* force vector */, Array & /* D_alpha */, Real & /* gap */) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// penalty parameter Real epsilon; }; } // akantu #endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */