diff --git a/examples/contact_mechanics/test_node_selector.cc b/examples/contact_mechanics/test_node_selector.cc new file mode 100644 index 000000000..3d641aa23 --- /dev/null +++ b/examples/contact_mechanics/test_node_selector.cc @@ -0,0 +1,32 @@ +#include "aka_common.hh" +#include "node_selector.hh" +#include "contact_mechanics_model.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char * argv[]) { + + const UInt spatial_dimension = 3; + + initialize("material.dat", argc, argv); + + Mesh mesh(spatial_dimension); + mesh.read("contact_hertz_2d.msh"); + + ContactMechanicsModel model(mesh); + + PhysicalSurfaceNodeSelector selector(model); + auto & slave = selector.getSlaveList(); + auto & master = selector.getMasterList(); + + for (auto & s : slave) { + std::cerr << s << std::endl; + } + + for (auto m : master) { + std::cerr << m << std::endl; + } + + return 0; +} diff --git a/packages/contact_mechanics.cmake b/packages/contact_mechanics.cmake index db13ad023..964f6f572 100644 --- a/packages/contact_mechanics.cmake +++ b/packages/contact_mechanics.cmake @@ -1,57 +1,60 @@ #=============================================================================== # @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_detector_inline_impl.cc model/contact_mechanics/contact_element.hh model/contact_mechanics/resolution.hh model/contact_mechanics/resolution.cc model/contact_mechanics/resolution_utils.hh model/contact_mechanics/resolution_utils.cc model/contact_mechanics/resolutions/resolution_penalty.hh - model/contact_mechanics/resolutions/resolution_penalty.cc) + model/contact_mechanics/resolutions/resolution_penalty.cc + + model/contact_mechanics/node_selector.hh + model/contact_mechanics/node_selector.cc) package_declare_documentation_files(contact_mechanics manual-contactmechanicsmodel.tex manual-contact-detector.tex manual-contact-resolutions.tex ) package_declare_documentation(contact_mechanics "This package activates the contact mechanics model") diff --git a/src/common/aka_bbox.hh b/src/common/aka_bbox.hh index cd0fe08b1..1df2cb61d 100644 --- a/src/common/aka_bbox.hh +++ b/src/common/aka_bbox.hh @@ -1,272 +1,273 @@ /** * @file aka_bbox.hh * * @author Nicolas Richart * * @date creation Mon Feb 12 2018 * * @brief A simple bounding box class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ +#include "aka_common.hh" #include "aka_iterators.hh" #include "aka_types.hh" #include "communicator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_BBOX_HH__ #define __AKANTU_AKA_BBOX_HH__ namespace akantu { class BBox { public: BBox() = default; BBox(UInt spatial_dimension) : dim(spatial_dimension), lower_bounds(spatial_dimension, std::numeric_limits::max()), upper_bounds(spatial_dimension, std::numeric_limits::lowest()) {} BBox(const BBox & other) : dim(other.dim), empty{false}, lower_bounds(other.lower_bounds), upper_bounds(other.upper_bounds) {} BBox & operator=(const BBox & other) { if (this != &other) { this->dim = other.dim; this->lower_bounds = other.lower_bounds; this->upper_bounds = other.upper_bounds; this->empty = other.empty; } return *this; } inline BBox & operator+=(const Vector & position) { AKANTU_DEBUG_ASSERT( this->dim == position.size(), "You are adding a point of a wrong dimension to the bounding box"); this->empty = false; for (auto s : arange(dim)) { lower_bounds(s) = std::min(lower_bounds(s), position(s)); upper_bounds(s) = std::max(upper_bounds(s), position(s)); } return *this; } /* ------------------------------------------------------------------------ */ inline bool intersects(const BBox & other, const SpatialDirection & direction) const { AKANTU_DEBUG_ASSERT( this->dim == other.dim, "You are intersecting bounding boxes of different dimensions"); return Math::intersects(lower_bounds(direction), upper_bounds(direction), other.lower_bounds(direction), other.upper_bounds(direction)); } inline bool intersects(const BBox & other) const { if (this->empty or other.empty) return false; bool intersects_ = true; for (auto s : arange(this->dim)) { intersects_ &= this->intersects(other, SpatialDirection(s)); } return intersects_; } /* ------------------------------------------------------------------------ */ inline BBox intersection(const BBox & other) const { AKANTU_DEBUG_ASSERT( this->dim == other.dim, "You are intersecting bounding boxes of different dimensions"); BBox intersection_(this->dim); intersection_.empty = not this->intersects(other); if (intersection_.empty) return intersection_; for (auto s : arange(this->dim)) { // is lower point in range ? bool point1 = Math::is_in_range(other.lower_bounds(s), lower_bounds(s), upper_bounds(s)); // is upper point in range ? bool point2 = Math::is_in_range(other.upper_bounds(s), lower_bounds(s), upper_bounds(s)); if (point1 and not point2) { // |-----------| this (i) // |-----------| other(i) // 1 2 intersection_.lower_bounds(s) = other.lower_bounds(s); intersection_.upper_bounds(s) = upper_bounds(s); } else if (point1 && point2) { // |-----------------| this (i) // |-----------| other(i) // 1 2 intersection_.lower_bounds(s) = other.lower_bounds(s); intersection_.upper_bounds(s) = other.upper_bounds(s); } else if (!point1 && point2) { // |-----------| this (i) // |-----------| other(i) // 1 2 intersection_.lower_bounds(s) = this->lower_bounds(s); intersection_.upper_bounds(s) = other.upper_bounds(s); } else { // |-----------| this (i) // |-----------------| other(i) // 1 2 intersection_.lower_bounds(s) = this->lower_bounds(s); intersection_.upper_bounds(s) = this->upper_bounds(s); } } return intersection_; } /* ------------------------------------------------------------------------ */ inline bool contains(const Vector & point) const { return (point >= lower_bounds) and (point <= upper_bounds); } /* ------------------------------------------------------------------------ */ inline void reset() { lower_bounds.set(std::numeric_limits::max()); upper_bounds.set(std::numeric_limits::lowest()); } /* -------------------------------------------------------------------------- */ inline void getCenter(Vector & center) { center = upper_bounds; center += lower_bounds; center /= 2.; } /* ------------------------------------------------------------------------ */ const Vector & getLowerBounds() const { return lower_bounds; } const Vector & getUpperBounds() const { return upper_bounds; } Vector & getLowerBounds() { return lower_bounds; } Vector & getUpperBounds() { return upper_bounds; } /* ------------------------------------------------------------------------ */ inline Real size(const SpatialDirection & direction) const { return upper_bounds(direction) - lower_bounds(direction); } Vector size() const { Vector size_(dim); for (auto s : arange(this->dim)) { size_(s) = this->size(SpatialDirection(s)); } return size_; } inline operator bool() const { return not empty; } /* ------------------------------------------------------------------------ */ BBox allSum(const Communicator & communicator) const { Matrix reduce_bounds(dim, 2); Vector(reduce_bounds(0)) = lower_bounds; Vector(reduce_bounds(1)) = Real(-1.) * upper_bounds; communicator.allReduce(reduce_bounds, SynchronizerOperation::_min); BBox global(dim); global.lower_bounds = Vector(reduce_bounds(0)); global.upper_bounds = Real(-1.) * Vector(reduce_bounds(1)); global.empty = false; return global; } std::vector allGather(const Communicator & communicator) const { auto prank = communicator.whoAmI(); auto nb_proc = communicator.getNbProc(); Array bboxes_data(nb_proc, dim * 2 + 1); auto * base = bboxes_data.storage() + prank * (2 * dim + 1); Vector(base + dim * 0, dim) = lower_bounds; Vector(base + dim * 1, dim) = upper_bounds; base[dim * 2] = empty ? 1. : 0.; // ugly trick communicator.allGather(bboxes_data); std::vector bboxes; bboxes.reserve(nb_proc); for (auto p : arange(nb_proc)) { bboxes.emplace_back(dim); auto & bbox = bboxes.back(); auto * base = bboxes_data.storage() + p * (2 * dim + 1); bbox.lower_bounds = Vector(base + dim * 0, dim); bbox.upper_bounds = Vector(base + dim * 1, dim); bbox.empty = base[dim * 2] == 1. ? true : false; } return bboxes; } std::map intersection(const BBox & other, const Communicator & communicator) { // todo: change for a custom reduction algorithm auto other_bboxes = other.allGather(communicator); std::map intersections; for (const auto & bbox : enumerate(other_bboxes)) { auto && tmp = this->intersection(std::get<1>(bbox)); if (tmp) { intersections[std::get<0>(bbox)] = tmp; } } return intersections; } void printself(std::ostream & stream) const { stream << "BBox["; if (not empty) { stream << lower_bounds << " - " << upper_bounds; } stream << "]"; } protected: UInt dim{0}; bool empty{true}; Vector lower_bounds; Vector upper_bounds; }; inline std::ostream & operator<<(std::ostream & stream, const BBox & bbox) { bbox.printself(stream); return stream; } } // akantu #endif /* __AKANTU_AKA_BBOX_HH__ */ diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 3f5407f52..bb7edc6b1 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,634 +1,637 @@ /** * @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 #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; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "aka_enum_macros.hh" /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of /// events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) \ (contact_mechanics_model) \ (coupler_solid_contact) // 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) /// 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, _explicit_dynamic_contact = 6 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_NON_LINEAR_SOLVER_TYPES \ (linear) \ (newton_raphson) \ (newton_raphson_modified) \ (lumped) \ (gmres) \ (bfgs) \ (cg) \ (newton_raphson_contact) \ (auto) // clang-format on AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) #else /// Type of non linear resolution available in akantu enum class NonLinearSolverType { _linear, ///< No non linear convergence loop _newton_raphson, ///< Regular Newton-Raphson _newton_raphson_modified, ///< Newton-Raphson with initial tangent _lumped, ///< Case of lumped mass or equivalent matrix _gmres, _bfgs, _cg, _newton_raphson_contact, ///< Regular Newton-Raphson modified /// for contact problem _auto, ///< This will take a default value that make sense in case of /// model::getNewSolver }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_TIME_STEP_SOLVER_TYPE \ (static) \ (dynamic) \ (dynamic_lumped) \ (not_defined) // clang-format on AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) #else /// Type of time stepping solver enum class TimeStepSolverType { _static, ///< Static solution _dynamic, ///< Dynamic solver _dynamic_lumped, ///< Dynamic solver with lumped mass _not_defined, ///< For not defined cases }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_INTEGRATION_SCHEME_TYPE \ (pseudo_time) \ (forward_euler) \ (trapezoidal_rule_1) \ (backward_euler) \ (central_difference) \ (fox_goodwin) \ (trapezoidal_rule_2) \ (linear_acceleration) \ (newmark_beta) \ (generalized_trapezoidal) // clang-format on AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) #else /// Type of integration scheme enum class IntegrationSchemeType { _pseudo_time, ///< Pseudo Time _forward_euler, ///< GeneralizedTrapezoidal(0) _trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _backward_euler, ///< GeneralizedTrapezoidal(1) _central_difference, ///< NewmarkBeta(0, 1/2) _fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ (residual) \ (solution) \ (residual_mass_wgh) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) #else /// enum SolveConvergenceCriteria different convergence criteria enum class SolveConvergenceCriteria { _residual, ///< Use residual to test the convergence _solution, ///< Use solution to test the convergence _residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; #endif /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /// @enum Type of contact detection enum DetectionType { _explicit, _implicit}; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_SYNCHRONIZATION_TAG \ (whatever) \ (update) \ (ask_nodes) \ (size) \ (smm_mass) \ (smm_for_gradu) \ (smm_boundary) \ (smm_uv) \ (smm_res) \ (smm_init_mat) \ (smm_stress) \ (smmc_facets) \ (smmc_facets_conn) \ (smmc_facets_stress) \ (smmc_damage) \ (giu_global_conn) \ (ce_groups) \ (gm_clusters) \ (htm_temperature) \ (htm_gradient_temperature) \ (htm_phi) \ (htm_gradient_phi) \ (mnl_for_average) \ (mnl_weight) \ (nh_criterion) \ (test) \ (user_1) \ (user_2) \ (material_id) \ (for_dump) \ (cf_nodal) \ (cf_incr) \ (solver_solution) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) #else /// @enum SynchronizationTag type of synchronizations enum class SynchronizationTag { //--- Generic tags --- _whatever, _update, _ask_nodes, _size, //--- SolidMechanicsModel tags --- _smm_mass, ///< synchronization of the SolidMechanicsModel.mass _smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _smm_uv, ///< synchronization of the nodal velocities and displacement _smm_res, ///< synchronization of the nodal residual _smm_init_mat, ///< synchronization of the data to initialize materials _smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _smmc_facets, ///< synchronization of facet data to setup facet synch _smmc_facets_conn, ///< synchronization of facet global connectivity _smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _htm_temperature, ///< synchronization of the nodal temperature _htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _htm_phi, ///< synchronization of the nodal level set value phi _htm_gradient_phi, ///< synchronization of the element gradient phi + _pfm_damage, + _pfm_gradient_damage, + //--- Material non local --- _mnl_for_average, ///< synchronization of data to average in non local /// material _mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _nh_criterion, // --- General tags --- _test, ///< Test tag _user_1, ///< tag for user simulations _user_2, ///< tag for user simulations _material_id, ///< synchronization of the material ids _for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _cf_nodal, ///< synchronization of disp, velo, and current position _cf_incr, ///< synchronization of increment // --- Solver tags --- _solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; #endif /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } std::ostream & operator<<(std::ostream & stream, NodeFlag flag); } // namespace akantu AKANTU_ENUM_HASH(GhostType) namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; 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 /* -------------------------------------------------------------------------- */ #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); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ struct TensorTrait {}; } // namespace akantu /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ namespace aka { /* ------------------------------------------------------------------------ */ template using is_tensor = std::is_base_of; /* ------------------------------------------------------------------------ */ template using is_scalar = std::is_arithmetic; /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> bool is_of_type(T && t) { return ( dynamic_cast>::value, std::add_const_t, R>>>(&t) != nullptr); } /* -------------------------------------------------------------------------- */ template bool is_of_type(std::unique_ptr & t) { return ( dynamic_cast::value, std::add_const_t, R>>>( t.get()) != nullptr); } /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { static_assert( disjunction< std::is_base_of, std::decay_t>, // down-cast std::is_base_of, std::decay_t> // up-cast >::value, "Type T and R are not valid for a as_type conversion"); return dynamic_cast>::value, std::add_const_t, R>>>(t); } /* -------------------------------------------------------------------------- */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { return &as_type(*t); } /* -------------------------------------------------------------------------- */ template decltype(auto) as_type(const std::shared_ptr & t) { return std::dynamic_pointer_cast(t); } } // namespace aka #include "aka_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/fe_engine/shape_lagrange_base.cc b/src/fe_engine/shape_lagrange_base.cc index 71d0f3942..f8484b988 100644 --- a/src/fe_engine/shape_lagrange_base.cc +++ b/src/fe_engine/shape_lagrange_base.cc @@ -1,160 +1,163 @@ /** * @file shape_lagrange_base.cc * * @author Nicolas Richart * * @date creation: Wed Aug 09 2017 * @date last modification: Tue Feb 20 2018 * * @brief common par for the shape lagrange * * @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 "shape_lagrange_base.hh" #include "mesh_iterators.hh" /* -------------------------------------------------------------------------- */ namespace akantu { ShapeLagrangeBase::ShapeLagrangeBase(const Mesh & mesh, const ElementKind & kind, const ID & id, const MemoryID & memory_id) : ShapeFunctions(mesh, id, memory_id), _kind(kind) {} /* -------------------------------------------------------------------------- */ ShapeLagrangeBase::~ShapeLagrangeBase() = default; /* -------------------------------------------------------------------------- */ #define AKANTU_COMPUTE_SHAPES(type) \ _this.template computeShapesOnIntegrationPoints( \ nodes, integration_points, shapes, ghost_type, filter_elements) namespace shape_lagrange { namespace details { template struct Helper { template static void call(const S &, const Array &, const Matrix &, Array &, const ElementType &, const GhostType &, const Array &) { AKANTU_TO_IMPLEMENT(); } }; #define AKANTU_COMPUTE_SHAPES_KIND(kind) \ template <> struct Helper { \ template \ static void call(const S & _this, const Array & nodes, \ const Matrix & integration_points, \ Array & shapes, const ElementType & type, \ const GhostType & ghost_type, \ const Array & filter_elements) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES, kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST(AKANTU_COMPUTE_SHAPES_KIND, AKANTU_FE_ENGINE_LIST_LAGRANGE_BASE) } // namespace details } // namespace shape_lagrange /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::computeShapesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shapes, const ElementType & type, const GhostType & ghost_type, const Array & filter_elements) const { auto kind = Mesh::getKind(type); #define AKANTU_COMPUTE_SHAPES_KIND_SWITCH(kind) \ shape_lagrange::details::Helper::call( \ *this, nodes, integration_points, shapes, type, ghost_type, \ filter_elements); AKANTU_BOOST_LIST_SWITCH( AKANTU_COMPUTE_SHAPES_KIND_SWITCH, BOOST_PP_LIST_TO_SEQ(AKANTU_FE_ENGINE_LIST_LAGRANGE_BASE), kind); #undef AKANTU_COMPUTE_SHAPES #undef AKANTU_COMPUTE_SHAPES_KIND #undef AKANTU_COMPUTE_SHAPES_KIND_SWITCH } /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::onElementsAdded(const Array & new_elements) { AKANTU_DEBUG_IN(); const auto & nodes = mesh.getNodes(); for (auto elements_range : MeshElementsByTypes(new_elements)) { auto type = elements_range.getType(); auto ghost_type = elements_range.getGhostType(); if (mesh.getKind(type) != _kind) continue; auto & elements = elements_range.getElements(); auto itp_type = FEEngine::getInterpolationType(type); if (not this->shapes_derivatives.exists(itp_type, ghost_type)) { auto size_of_shapesd = this->getShapeDerivativesSize(type); this->shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type); } if (not shapes.exists(itp_type, ghost_type)) { auto size_of_shapes = this->getShapeSize(type); this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type); } const auto & natural_coords = integration_points(type, ghost_type); computeShapesOnIntegrationPoints(nodes, natural_coords, shapes(itp_type, ghost_type), type, ghost_type, elements); - computeShapeDerivativesOnIntegrationPoints( - nodes, natural_coords, shapes_derivatives(itp_type, ghost_type), type, - ghost_type, elements); + if (mesh.getSpatialDimension() == mesh.getNaturalSpaceDimension(type) || + _kind != _ek_regular) { + computeShapeDerivativesOnIntegrationPoints( + nodes, natural_coords, shapes_derivatives(itp_type, ghost_type), type, + ghost_type, elements); + } } #undef INIT_SHAPE_FUNCTIONS AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::onElementsRemoved( const Array &, const ElementTypeMapArray & new_numbering) { this->shapes.onElementsRemoved(new_numbering); this->shapes_derivatives.onElementsRemoved(new_numbering); } /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Shapes Lagrange [" << std::endl; ShapeFunctions::printself(stream, indent + 1); shapes.printself(stream, indent + 1); shapes_derivatives.printself(stream, indent + 1); stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index 002b01a60..f04ac6468 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,201 +1,201 @@ /** * @file element_group.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Wed Nov 08 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @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_memory.hh" #include "dumpable.hh" #include "element_type_map.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_GROUP_HH__ #define __AKANTU_ELEMENT_GROUP_HH__ namespace akantu { class Mesh; class Element; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class ElementGroup : private Memory, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ElementGroup(const std::string & name, const Mesh & mesh, NodeGroup & node_group, UInt dimension = _all_dimensions, const std::string & id = "element_group", const MemoryID & memory_id = 0); ElementGroup(const ElementGroup &); - + /* ------------------------------------------------------------------------ */ /* Type definitions */ /* ------------------------------------------------------------------------ */ public: using ElementList = ElementTypeMapArray; using NodeList = Array; -/* ------------------------------------------------------------------------ */ -/* Element iterator */ -/* ------------------------------------------------------------------------ */ + /* ------------------------------------------------------------------------ */ + /* Element iterator */ + /* ------------------------------------------------------------------------ */ using type_iterator = ElementList::type_iterator; [[deprecated("Use elementTypes instead")]] inline type_iterator firstType(UInt dim = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_regular) const; [[deprecated("Use elementTypes instead")]] inline type_iterator lastType(UInt dim = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_regular) const; template inline decltype(auto) elementTypes(pack &&... _pack) const { return elements.elementTypes(_pack...); } using const_element_iterator = Array::const_iterator; inline const_element_iterator begin(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline const_element_iterator end(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// empty the element group void empty(); /// append another group to this group /// BE CAREFUL: it doesn't conserve the element order void append(const ElementGroup & other_group); /// add an element to the group. By default the it does not add the nodes to /// the group inline void add(const Element & el, bool add_nodes = false, bool check_for_duplicate = true); /// \todo fix the default for add_nodes : make it coherent with the other /// method inline void add(const ElementType & type, UInt element, const GhostType & ghost_type = _not_ghost, bool add_nodes = true, bool check_for_duplicate = true); inline void addNode(UInt node_id, bool check_for_duplicate = true); inline void removeNode(UInt node_id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// fill the elements based on the underlying node group. virtual void fillFromNodeGroup(); // sort and remove duplicated values void optimize(); /// change the dimension if needed void addDimension(UInt dimension); private: inline void addElement(const ElementType & elem_type, UInt elem_id, const GhostType & ghost_type); friend class GroupManager; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: const Array & getElements(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; AKANTU_GET_MACRO(Elements, elements, const ElementTypeMapArray &); AKANTU_GET_MACRO_NOT_CONST(Elements, elements, ElementTypeMapArray &); - -// AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array &); + + // AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array &); AKANTU_GET_MACRO(NodeGroup, node_group, const NodeGroup &); AKANTU_GET_MACRO_NOT_CONST(NodeGroup, node_group, NodeGroup &); AKANTU_GET_MACRO(Dimension, dimension, UInt); AKANTU_GET_MACRO(Name, name, std::string); inline UInt getNbNodes() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Mesh to which this group belongs const Mesh & mesh; /// name of the group std::string name; /// list of elements composing the group ElementList elements; /// sub list of nodes which are composing the elements NodeGroup & node_group; /// group dimension UInt dimension; /// empty arry for the iterator to work when an element type not present Array empty_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const ElementGroup & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "element.hh" #include "element_group_inline_impl.cc" #endif /* __AKANTU_ELEMENT_GROUP_HH__ */ diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh index 0ca606684..522337572 100644 --- a/src/mesh/mesh.hh +++ b/src/mesh/mesh.hh @@ -1,701 +1,704 @@ /** * @file mesh.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Feb 19 2018 * * @brief the class representing the meshes * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_bbox.hh" #include "aka_event_handler_manager.hh" #include "aka_memory.hh" #include "communicator.hh" #include "dumpable.hh" #include "element.hh" #include "element_class.hh" #include "element_type_map.hh" #include "group_manager.hh" #include "mesh_data.hh" #include "mesh_events.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementSynchronizer; class NodeSynchronizer; class PeriodicNodeSynchronizer; class MeshGlobalDataUpdater; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(communicator); DECLARE_NAMED_ARGUMENT(edge_weight_function); DECLARE_NAMED_ARGUMENT(vertex_weight_function); } // namespace /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Array, and the connectivity. The connectivity are stored in by element * types. * * In order to loop on all element you have to loop on all types like this : * @code{.cpp} for(auto & type : mesh.elementTypes()) { UInt nb_element = mesh.getNbElement(type); const Array & conn = mesh.getConnectivity(type); for(UInt e = 0; e < nb_element; ++e) { ... } } or for_each_element(mesh, [](Element & element) { std::cout << element << std::endl }); @endcode */ class Mesh : protected Memory, public EventHandlerManager, public GroupManager, public MeshData, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: /// default constructor used for chaining, the last parameter is just to /// differentiate constructors Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator); public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// mesh not distributed and not using the default communicator Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); ~Mesh() override; /// read the mesh from a file void read(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); /// write the mesh to a file void write(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); protected: void makeReady(); private: /// initialize the connectivity to NULL and other stuff void init(); /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /* ------------------------------------------------------------------------ */ /* Distributed memory methods and accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// patitionate the mesh among the processors involved in their computation virtual void distributeImpl( Communicator & communicator, std::function edge_weight_function, std::function vertex_weight_function); public: /// with the arguments to pass to the partitionner template std::enable_if_t::value> distribute(pack &&... _pack) { distributeImpl( OPTIONAL_NAMED_ARG(communicator, Communicator::getStaticCommunicator()), OPTIONAL_NAMED_ARG(edge_weight_function, [](auto &&, auto &&) { return 1; }), OPTIONAL_NAMED_ARG(vertex_weight_function, [](auto &&) { return 1; })); } /// defines is the mesh is distributed or not inline bool isDistributed() const { return this->is_distributed; } /* ------------------------------------------------------------------------ */ /* Periodicity methods and accessors */ /* ------------------------------------------------------------------------ */ public: /// set the periodicity in a given direction void makePeriodic(const SpatialDirection & direction); void makePeriodic(const SpatialDirection & direction, const ID & list_1, const ID & list_2); protected: void makePeriodic(const SpatialDirection & direction, const Array & list_1, const Array & list_2); /// Removes the face that the mesh is periodic void wipePeriodicInfo(); inline void addPeriodicSlave(UInt slave, UInt master); template void synchronizePeriodicSlaveDataWithMaster(Array & data); // update the periodic synchronizer (creates it if it does not exists) void updatePeriodicSynchronizer(); public: /// defines if the mesh is periodic or not inline bool isPeriodic() const { return (this->is_periodic != 0); } inline bool isPeriodic(const SpatialDirection & direction) const { return ((this->is_periodic & (1 << direction)) != 0); } class PeriodicSlaves; /// get the master node for a given slave nodes, except if node not a slave inline UInt getPeriodicMaster(UInt slave) const; /// get an iterable list of slaves for a given master node inline decltype(auto) getPeriodicSlaves(UInt master) const; /* ------------------------------------------------------------------------ */ /* General Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Array & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; // /// extract coordinates of nodes from a reversed element // inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, // UInt * connectivity, // UInt n_nodes); /// add a Array of connectivity for the type . inline void addConnectivityType(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ template inline void sendEvent(Event & event) { // if(event.getList().size() != 0) EventHandlerManager::sendEvent(event); } /// prepare the event to remove the elements listed void eraseElements(const Array & elements); /* ------------------------------------------------------------------------ */ template inline void removeNodesFromArray(Array & vect, const Array & new_numbering); /// initialize normals void initNormals(); /// init facets' mesh Mesh & initMeshFacets(const ID & id = "mesh_facets"); /// define parent mesh void defineMeshParent(const Mesh & mesh); /// get global connectivity array void getGlobalConnectivity(ElementTypeMapArray & global_connectivity); public: void getAssociatedElements(const Array & node_list, Array & elements); void getAssociatedElements(const UInt & node, Array & elements); public: /// fills the nodes_to_elements for given dimension elements void fillNodesToElements(UInt dimension = _all_dimensions); private: /// update the global ids, nodes type, ... std::tuple updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event); void registerGlobalDataUpdater( std::unique_ptr && global_data_updater); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the id of the mesh AKANTU_GET_MACRO(ID, Memory::id, const ID &); /// get the id of the mesh AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &); /// get the spatial dimension of the mesh = number of component of the /// coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Array aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Array &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array &); /// get the normals for the elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt); /// get the Array of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array &); // AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array // &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global id of a node inline UInt getNodeLocalId(UInt global_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Array AKANTU_GET_MACRO(NodesFlags, *nodes_flags, const Array &); protected: AKANTU_GET_MACRO_NOT_CONST(NodesFlags, *nodes_flags, Array &); public: inline NodeFlag getNodeFlag(UInt local_id) const; inline Int getNodePrank(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; inline bool isPeriodicSlave(UInt n) const; inline bool isPeriodicMaster(UInt n) const; const Vector & getLowerBounds() const { return bbox.getLowerBounds(); } const Vector & getUpperBounds() const { return bbox.getUpperBounds(); } AKANTU_GET_MACRO(BBox, bbox, const BBox &); const Vector & getLocalLowerBounds() const { return bbox_local.getLowerBounds(); } const Vector & getLocalUpperBounds() const { return bbox_local.getUpperBounds(); } AKANTU_GET_MACRO(LocalBBox, bbox_local, const BBox &); /// get the connectivity Array for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); AKANTU_GET_MACRO(Connectivities, connectivities, const ElementTypeMapArray &); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get the number of element for a given ghost_type and a given dimension inline UInt getNbElement(const UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_not_defined) const; /// compute the barycenter of a given element inline void getBarycenter(const Element & element, Vector & barycenter) const; void getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const; /// get the element connected to a subelement (element of lower dimension) const auto & getElementToSubelement() const; /// get the element connected to a subelement const auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the element connected to a subelement auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the elements connected to a subelement const auto & getElementToSubelement(const Element & element) const; /// get the subelement (element of lower dimension) connected to a element const auto & getSubelementToElement() const; /// get the subelement connected to an element const auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the subelement connected to an element auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the subelement (element of lower dimension) connected to a element VectorProxy getSubelementToElement(const Element & element) const; /// get connectivity of a given element inline VectorProxy getConnectivity(const Element & element) const; inline Vector getConnectivityWithPeriodicity(const Element & element) const; protected: inline auto & getElementToSubelement(const Element & element); inline VectorProxy getSubelementToElement(const Element & element); inline VectorProxy getConnectivity(const Element & element); public: /// get a name field associated to the mesh template inline const Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get a name field associated to the mesh template inline Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get a name field associated to the mesh template inline const ElementTypeMapArray & getData(const ID & data_name) const; /// get a name field associated to the mesh template inline ElementTypeMapArray & getData(const ID & data_name); template ElementTypeMap getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); template std::shared_ptr createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); /// templated getter returning the pointer to data in MeshData (modifiable) template inline Array & getDataPointer(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost, UInt nb_component = 1, bool size_to_nb_element = true, bool resize_with_parent = false); template inline Array & getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_); /// Facets mesh accessor inline const Mesh & getMeshFacets() const; inline Mesh & getMeshFacets(); /// Parent mesh accessor inline const Mesh & getMeshParent() const; inline bool isMeshFacets() const { return this->is_mesh_facets; } /// return the dumper from a group and and a dumper name DumperIOHelper & getGroupDumper(const std::string & dumper_name, const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); + /// get the natural space dimension of a type of element + static inline UInt getNaturalSpaceDimension(const ElementType & type); + /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type, UInt t); /// get local connectivity of a facet for a given facet type static inline auto getFacetLocalConnectivity(const ElementType & type, UInt t = 0); /// get connectivity of facets for a given element inline auto getFacetConnectivity(const Element & element, UInt t = 0) const; /// get the number of type of the surface element associated to a given /// element type static inline UInt getNbFacetTypes(const ElementType & type, UInt t = 0); /// get the type of the surface element associated to a given element static inline constexpr auto getFacetType(const ElementType & type, UInt t = 0); /// get all the type of the surface element associated to a given element static inline constexpr auto getAllFacetTypes(const ElementType & type); /// get the number of nodes in the given element list static inline UInt getNbNodesPerElementList(const Array & elements); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ using type_iterator [[deprecated]] = ElementTypeMapArray::type_iterator; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; template ElementTypesIteratorHelper elementTypes(pack &&... _pack) const; [[deprecated("Use elementTypes instead")]] inline decltype(auto) firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).begin(); } [[deprecated("Use elementTypes instead")]] inline decltype(auto) lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).end(); } AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer, const ElementSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer, ElementSynchronizer &); AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer, const NodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer, NodeSynchronizer &); AKANTU_GET_MACRO(PeriodicNodeSynchronizer, *periodic_node_synchronizer, const PeriodicNodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(PeriodicNodeSynchronizer, *periodic_node_synchronizer, PeriodicNodeSynchronizer &); // AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator // &); AKANTU_GET_MACRO(Communicator, *communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &); AKANTU_GET_MACRO(PeriodicMasterSlaves, periodic_master_slave, const auto &); /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshAccessor; friend class MeshUtils; AKANTU_GET_MACRO(NodesPointer, *nodes, Array &); /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline Array & getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Array and create it if necessary inline Array & getNodesFlagsPointer(); /// get a pointer to the connectivity Array for the given type and create it /// if necessary inline Array & getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the ghost element counter inline Array & getGhostsCounters(const ElementType & type, const GhostType & ghost_type = _ghost) { AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost, "No ghost counter for _not_ghost elements"); return ghosts_counters(type, ghost_type); } /// get a pointer to the element_to_subelement Array for the given type and /// create it if necessary inline Array> & getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Array for the given type and /// create it if necessary inline Array & getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of the nodes coordinates std::shared_ptr> nodes; /// global node ids std::shared_ptr> nodes_global_ids; /// node flags (shared/periodic/...) std::shared_ptr> nodes_flags; /// processor handling the node when not local or master std::unordered_map nodes_prank; /// global number of nodes; UInt nb_global_nodes{0}; /// all class of elements present in this mesh (for heterogenous meshes) ElementTypeMapArray connectivities; /// count the references on ghost elements ElementTypeMapArray ghosts_counters; /// map to normals for all class of elements present in this mesh ElementTypeMapArray normals; /// the spatial dimension of this mesh UInt spatial_dimension{0}; /// size covered by the mesh on each direction Vector size; /// global bounding box BBox bbox; /// local bounding box BBox bbox_local; /// Extra data loaded from the mesh file // MeshData mesh_data; /// facets' mesh std::unique_ptr mesh_facets; /// parent mesh (this is set for mesh_facets meshes) const Mesh * mesh_parent{nullptr}; /// defines if current mesh is mesh_facets or not bool is_mesh_facets{false}; /// defines if the mesh is centralized or distributed bool is_distributed{false}; /// defines if the mesh is periodic bool is_periodic{false}; /// Communicator on which mesh is distributed Communicator * communicator; /// Element synchronizer std::unique_ptr element_synchronizer; /// Node synchronizer std::unique_ptr node_synchronizer; /// Node synchronizer for periodic nodes std::unique_ptr periodic_node_synchronizer; using NodesToElements = std::vector>>; /// class to update global data using external knowledge std::unique_ptr global_data_updater; /// This info is stored to simplify the dynamic changes NodesToElements nodes_to_elements; /// periodicity local info std::unordered_map periodic_slave_master; std::unordered_multimap periodic_master_slave; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ #include "element_type_map_tmpl.hh" #include "mesh_inline_impl.cc" #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/mesh/mesh_inline_impl.cc b/src/mesh/mesh_inline_impl.cc index af4c9e1ef..e7f395620 100644 --- a/src/mesh/mesh_inline_impl.cc +++ b/src/mesh/mesh_inline_impl.cc @@ -1,767 +1,775 @@ /** * @file mesh_inline_impl.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Thu Jul 15 2010 * @date last modification: Mon Dec 18 2017 * * @brief Implementation of the inline functions of the mesh 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 "aka_iterators.hh" #include "element_class.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_INLINE_IMPL_CC__ #define __AKANTU_MESH_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline ElementKind Element::kind() const { return Mesh::getKind(type); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template Mesh::ElementTypesIteratorHelper Mesh::elementTypes(pack &&... _pack) const { return connectivities.elementTypes(_pack...); } /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh) : new_numbering(mesh.getNbNodes(), 1, "new_numbering") {} /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh, const ID & new_numbering_id) : new_numbering(new_numbering_id, mesh.getID(), mesh.getMemoryID()) {} /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewElementsEvent & event) { - this->nodes_to_elements.resize(nodes->size()); + this->fillNodesToElements(); + /*this->nodes_to_elements.resize(nodes->size()); for (const auto & elem : event.getList()) { const Array & conn = connectivities(elem.type, elem.ghost_type); UInt nb_nodes_per_elem = this->getNbNodesPerElement(elem.type); - for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = conn(elem.element, n); if (not nodes_to_elements[node]) nodes_to_elements[node] = std::make_unique>(); nodes_to_elements[node]->insert(elem); } - } + }*/ EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewNodesEvent & event) { this->computeBoundingBox(); this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedElementsEvent & event) { this->connectivities.onElementsRemoved(event.getNewNumbering()); this->fillNodesToElements(); this->computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedNodesEvent & event) { const auto & new_numbering = event.getNewNumbering(); this->removeNodesFromArray(*nodes, new_numbering); if (nodes_global_ids and not is_mesh_facets) this->removeNodesFromArray(*nodes_global_ids, new_numbering); if (not is_mesh_facets) this->removeNodesFromArray(*nodes_flags, new_numbering); if (not nodes_to_elements.empty()) { std::vector>> tmp( nodes_to_elements.size()); auto it = nodes_to_elements.begin(); UInt new_nb_nodes = 0; for (auto new_i : new_numbering) { if (new_i != UInt(-1)) { tmp[new_i] = std::move(*it); ++new_nb_nodes; } ++it; } tmp.resize(new_nb_nodes); std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin()); } computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template inline void Mesh::removeNodesFromArray(Array & vect, const Array & new_numbering) { Array tmp(vect.size(), vect.getNbComponent()); UInt nb_component = vect.getNbComponent(); UInt new_nb_nodes = 0; for (UInt i = 0; i < new_numbering.size(); ++i) { UInt new_i = new_numbering(i); if (new_i != UInt(-1)) { T * to_copy = vect.storage() + i * nb_component; std::uninitialized_copy(to_copy, to_copy + nb_component, tmp.storage() + new_i * nb_component); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if (not nodes_global_ids) { nodes_global_ids = std::make_shared>( nodes->size(), 1, getID() + ":nodes_global_ids"); for (auto && global_ids : enumerate(*nodes_global_ids)) { std::get<1>(global_ids) = std::get<0>(global_ids); } } AKANTU_DEBUG_OUT(); return *nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getConnectivityPointer(const ElementType & type, const GhostType & ghost_type) { if (connectivities.exists(type, ghost_type)) return connectivities(type, ghost_type); if (ghost_type != _not_ghost) { ghosts_counters.alloc(0, 1, type, ghost_type, 1); } AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); return connectivities.alloc(0, nb_nodes_per_element, type, ghost_type); } /* -------------------------------------------------------------------------- */ inline Array> & Mesh::getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type) { return getDataPointer>("element_to_subelement", type, ghost_type, 1, true); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type) { auto & array = getDataPointer( "subelement_to_element", type, ghost_type, getNbFacetsPerElement(type), true, is_mesh_facets, ElementNull); return array; } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement() const { return getData>("element_to_subelement"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) const { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(const Element & element) const { return getData>("element_to_subelement")(element); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelement(const Element & element) { return getData>("element_to_subelement")(element); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement() const { return getData("subelement_to_element"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) const { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElement(const Element & element) const { const auto & sub_to_element = this->getSubelementToElement(element.type, element.ghost_type); auto it = sub_to_element.begin(sub_to_element.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElement(const Element & element) { auto & sub_to_element = this->getSubelementToElement(element.type, element.ghost_type); auto it = sub_to_element.begin(sub_to_element.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent) { Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) tmp.resize(mesh_parent->getNbElement(el_type, ghost_type)); else tmp.resize(this->getNbElement(el_type, ghost_type)); } else { tmp.resize(0); } return tmp; } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_) { Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) tmp.resize(mesh_parent->getNbElement(el_type, ghost_type), defaul_); else tmp.resize(this->getNbElement(el_type, ghost_type), defaul_); } else { tmp.resize(0); } return tmp; } /* -------------------------------------------------------------------------- */ template inline const Array & Mesh::getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type) const { return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type) { return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline const ElementTypeMapArray & Mesh::getData(const ID & data_name) const { return this->getElementalData(data_name); } /* -------------------------------------------------------------------------- */ template inline ElementTypeMapArray & Mesh::getData(const ID & data_name) { return this->getElementalData(data_name); } - /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const ElementType & type, const GhostType & ghost_type) const { try { const Array & conn = connectivities(type, ghost_type); return conn.size(); } catch (...) { return 0; } } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & kind) const { AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == UInt(-1), "spatial_dimension is " << spatial_dimension << " and is greater than 3 !"); UInt nb_element = 0; for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) nb_element += getNbElement(type, ghost_type); return nb_element; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(const Element & element, Vector & barycenter) const { Vector conn = getConnectivity(element); Matrix local_coord(spatial_dimension, conn.size()); auto node_begin = make_view(*nodes, spatial_dimension).begin(); for (auto && node : enumerate(conn)) { local_coord(std::get<0>(node)) = Vector(node_begin[std::get<1>(node)]); } Math::barycenter(local_coord.storage(), conn.size(), spatial_dimension, barycenter.storage()); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElement(const ElementType & type) { UInt nb_nodes_per_element = 0; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass::getNbNodesPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT); #undef GET_NB_NODES_PER_ELEMENT return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getP1ElementType(const ElementType & type) { ElementType p1_type = _not_defined; #define GET_P1_TYPE(type) p1_type = ElementClass::getP1ElementType() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_P1_TYPE); #undef GET_P1_TYPE return p1_type; } /* -------------------------------------------------------------------------- */ inline ElementKind Mesh::getKind(const ElementType & type) { ElementKind kind = _ek_not_defined; #define GET_KIND(type) kind = ElementClass::getKind() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_KIND); #undef GET_KIND return kind; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getSpatialDimension(const ElementType & type) { UInt spatial_dimension = 0; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass::getSpatialDimension() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION); #undef GET_SPATIAL_DIMENSION return spatial_dimension; } +/* -------------------------------------------------------------------------- */ +inline UInt Mesh::getNaturalSpaceDimension(const ElementType & type) { + UInt natural_dimension = 0; +#define GET_NATURAL_DIMENSION(type) \ + natural_dimension = ElementClass::getNaturalSpaceDimension() + AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_DIMENSION); +#undef GET_NATURAL_DIMENSION + + return natural_dimension; +} + /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetTypes(const ElementType & type, __attribute__((unused)) UInt t) { UInt nb = 0; #define GET_NB_FACET_TYPE(type) nb = ElementClass::getNbFacetTypes() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET_TYPE); #undef GET_NB_FACET_TYPE return nb; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getFacetType(const ElementType & type, UInt t) { #define GET_FACET_TYPE(type) return ElementClass::getFacetType(t); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return _not_defined; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getAllFacetTypes(const ElementType & type) { #define GET_FACET_TYPE(type) return ElementClass::getFacetTypes(); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return ElementClass<_not_defined>::getFacetTypes(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) n_facet = ElementClass::getNbFacetsPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetLocalConnectivity(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); #define GET_FACET_CON(type) \ AKANTU_DEBUG_OUT(); \ return ElementClass::getFacetLocalConnectivityPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON); #undef GET_FACET_CON AKANTU_DEBUG_OUT(); return ElementClass<_not_defined>::getFacetLocalConnectivityPerElement(0); // This avoid a compilation warning but will certainly // also cause a segfault if reached } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetConnectivity(const Element & element, UInt t) const { AKANTU_DEBUG_IN(); Matrix local_facets(getFacetLocalConnectivity(element.type, t)); Matrix facets(local_facets.rows(), local_facets.cols()); const Array & conn = connectivities(element.type, element.ghost_type); for (UInt f = 0; f < facets.rows(); ++f) { for (UInt n = 0; n < facets.cols(); ++n) { facets(f, n) = conn(element.element, local_facets(f, n)); } } AKANTU_DEBUG_OUT(); return facets; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivity(const Element & element) const { const auto & conn = connectivities(element.type, element.ghost_type); auto it = conn.begin(conn.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivity(const Element & element) { auto & conn = connectivities(element.type, element.ghost_type); auto it = conn.begin(conn.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ template inline void Mesh::extractNodalValuesFromElement( const Array & nodal_values, T * local_coord, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const { for (UInt n = 0; n < n_nodes; ++n) { memcpy(local_coord + n * nb_degree_of_freedom, nodal_values.storage() + connectivity[n] * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(T)); } } /* -------------------------------------------------------------------------- */ inline void Mesh::addConnectivityType(const ElementType & type, const GhostType & ghost_type) { getConnectivityPointer(type, ghost_type); } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPureGhostNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_pure_ghost; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_local_master_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_master; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_slave; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPeriodicSlave(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) == NodeFlag::_periodic_slave; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPeriodicMaster(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) == NodeFlag::_periodic_master; } /* -------------------------------------------------------------------------- */ inline NodeFlag Mesh::getNodeFlag(UInt local_id) const { return (*nodes_flags)(local_id); } /* -------------------------------------------------------------------------- */ inline Int Mesh::getNodePrank(UInt local_id) const { auto it = nodes_prank.find(local_id); return it == nodes_prank.end() ? -1 : it->second; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeLocalId(UInt global_id) const { if (nodes_global_ids == nullptr) return global_id; return nodes_global_ids->find(global_id); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->size(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElementList(const Array & elements) { UInt nb_nodes_per_element = 0; UInt nb_nodes = 0; ElementType current_element_type = _not_defined; for (const auto & el : elements) { if (el.type != current_element_type) { current_element_type = el.type; nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type); } nb_nodes += nb_nodes_per_element; } return nb_nodes; } /* -------------------------------------------------------------------------- */ inline Mesh & Mesh::getMeshFacets() { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshFacets() const { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshParent() const { if (!this->mesh_parent) AKANTU_SILENT_EXCEPTION( "No parent mesh is defined! This is only valid in a mesh_facets"); return *this->mesh_parent; } /* -------------------------------------------------------------------------- */ void Mesh::addPeriodicSlave(UInt slave, UInt master) { if (master == slave) return; // if pair already registered auto master_slaves = periodic_master_slave.equal_range(master); auto slave_it = std::find_if(master_slaves.first, master_slaves.second, [&](auto & pair) { return pair.second == slave; }); if (slave_it == master_slaves.second) { // no duplicates periodic_master_slave.insert(std::make_pair(master, slave)); AKANTU_DEBUG_INFO("adding periodic slave, slave gid:" << getNodeGlobalId(slave) << " [lid: " << slave << "]" << ", master gid:" << getNodeGlobalId(master) << " [lid: " << master << "]"); - // std::cout << "adding periodic slave, slave gid:" << getNodeGlobalId(slave) + // std::cout << "adding periodic slave, slave gid:" << + // getNodeGlobalId(slave) // << " [lid: " << slave << "]" // << ", master gid:" << getNodeGlobalId(master) // << " [lid: " << master << "]" << std::endl; } periodic_slave_master[slave] = master; auto set_flag = [&](auto node, auto flag) { (*nodes_flags)[node] &= ~NodeFlag::_periodic_mask; // clean periodic flags (*nodes_flags)[node] |= flag; }; set_flag(slave, NodeFlag::_periodic_slave); set_flag(master, NodeFlag::_periodic_master); } /* -------------------------------------------------------------------------- */ UInt Mesh::getPeriodicMaster(UInt slave) const { return periodic_slave_master.at(slave); } /* -------------------------------------------------------------------------- */ class Mesh::PeriodicSlaves { using internal_iterator = std::unordered_multimap::const_iterator; std::pair pair; public: PeriodicSlaves(const Mesh & mesh, UInt master) : pair(mesh.getPeriodicMasterSlaves().equal_range(master)) {} PeriodicSlaves(const PeriodicSlaves & other) = default; PeriodicSlaves(PeriodicSlaves && other) = default; PeriodicSlaves & operator=(const PeriodicSlaves & other) = default; class const_iterator { internal_iterator it; public: const_iterator(internal_iterator it) : it(std::move(it)) {} const_iterator operator++() { ++it; return *this; } - bool - operator!=(const const_iterator & other) { - return other.it != it; - } + bool operator!=(const const_iterator & other) { return other.it != it; } auto operator*() { return it->second; } }; auto begin() { return const_iterator(pair.first); } auto end() { return const_iterator(pair.second); } }; /* -------------------------------------------------------------------------- */ inline decltype(auto) Mesh::getPeriodicSlaves(UInt master) const { return PeriodicSlaves(*this, master); } /* -------------------------------------------------------------------------- */ inline Vector Mesh::getConnectivityWithPeriodicity(const Element & element) const { Vector conn = getConnectivity(element); if (not isPeriodic()) { return conn; } for (auto && node : conn) { if (isPeriodicSlave(node)) { node = getPeriodicMaster(node); } } return conn; } } // namespace akantu #endif /* __AKANTU_MESH_INLINE_IMPL_CC__ */ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index ed51811d6..f30b07edb 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,1817 +1,1847 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Leonardo Snozzi * @author Marco Vocialta * * @date creation: Fri Aug 20 2010 * @date last modification: Wed Feb 21 2018 * * @brief All mesh utils necessary for various tasks * * @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 "mesh_utils.hh" #include "element_synchronizer.hh" #include "fe_engine.hh" #include "mesh_accessor.hh" #include "mesh_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); /// count number of occurrence of each node UInt nb_nodes = mesh.getNbNodes(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); for_each_element(mesh, [&](auto && element) { Vector conn = mesh.getConnectivity(element); for (auto && node : conn) { ++node_to_elem.rowOffset(node); } }, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list // Element e; node_to_elem.beginInsertions(); for_each_element(mesh, [&](auto && element) { Vector conn = mesh.getConnectivity(element); for (auto && node : conn) { node_to_elem.insertInRow(node, element); } }, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh, CSR & node_to_elem, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, ghost_type).size(); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) ++node_to_elem.rowOffset(conn_val[el_offset + n]); } /// convert the occurrence array in a csr one node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(spatial_dimension - 1, ghost_type)) { mesh.getConnectivity(type, ghost_type).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt to_dimension) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt from_dimension, UInt to_dimension) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT( mesh_facets.isMeshFacets(), "The mesh_facets should be initialized with initMeshFacets"); /// generate facets buildFacetsDimension(mesh, mesh_facets, false, from_dimension); /// sort facets and generate sub-facets for (UInt i = from_dimension - 1; i > to_dimension; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension) { AKANTU_DEBUG_IN(); // save the current parent of mesh_facets and set it temporarly to mesh since // mesh is the one containing the elements for which mesh_facets has the // sub-elements // example: if the function is called with mesh = mesh_facets const Mesh * mesh_facets_parent = nullptr; try { mesh_facets_parent = &mesh_facets.getMeshParent(); } catch (...) { } mesh_facets.defineMeshParent(mesh); MeshAccessor mesh_accessor(mesh_facets); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const auto mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; std::vector connected_elements; // init the SubelementToElement data to improve performance for (auto && ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(dimension, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type); auto facet_types = mesh.getAllFacetTypes(type); for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); mesh_accessor.getElementToSubelement(facet_type, ghost_type); mesh_accessor.getConnectivity(facet_type, ghost_type); } } } const ElementSynchronizer * synchronizer = nullptr; if (mesh.isDistributed()) { synchronizer = &(mesh.getElementSynchronizer()); } Element current_element; for (auto && ghost_type : ghost_types) { GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; for (auto && type : mesh.elementTypes(dimension, ghost_type)) { auto facet_types = mesh.getAllFacetTypes(type); current_element.type = type; for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); auto nb_element = mesh.getNbElement(type, ghost_type); auto element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, ghost_type); auto connectivity_facets = &mesh_facets.getConnectivity(facet_type, ghost_type); auto nb_facet_per_element = mesh.getNbFacetsPerElement(type, ft); const auto & element_connectivity = mesh.getConnectivity(type, ghost_type); Matrix facet_local_connectivity( mesh.getFacetLocalConnectivity(type, ft)); auto nb_nodes_per_facet = connectivity_facets->getNbComponent(); Vector facet(nb_nodes_per_facet); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; for (UInt f = 0; f < nb_facet_per_element; ++f) { for (UInt n = 0; n < nb_nodes_per_facet; ++n) facet(n) = element_connectivity(el, facet_local_connectivity(f, n)); UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0)); counter.resize(first_node_nb_elements); counter.clear(); // loop over the other nodes to search intersecting elements, // which are the elements that share another node with the // starting element after first_node UInt local_el = 0; auto first_node_elements = node_to_elem.begin(facet(0)); auto first_node_elements_end = node_to_elem.end(facet(0)); for (; first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (UInt n = 1; n < nb_nodes_per_facet; ++n) { auto node_elements_begin = node_to_elem.begin(facet(n)); auto node_elements_end = node_to_elem.end(facet(n)); counter(local_el) += std::count(node_elements_begin, node_elements_end, *first_node_elements); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; connected_elements.clear(); for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) { Element real_el = node_to_elem(facet(0), el_f); if (not(counter(el_f) == nb_nodes_per_facet - 1)) continue; ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } if (minimum_el != current_element) continue; bool full_ghost_facet = false; UInt n = 0; while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n))) ++n; if (n == nb_nodes_per_facet) full_ghost_facet = true; if (full_ghost_facet) continue; if (boundary_only and nb_element_connected_to_facet != 1) continue; std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(current_element); if (nb_element_connected_to_facet == 1) { /// boundary facet elements.push_back(ElementNull); } else if (nb_element_connected_to_facet == 2) { /// internal facet elements.push_back(connected_elements[1]); /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2] = {_not_ghost, _not_ghost}; for (UInt el = 0; el < connected_elements.size(); ++el) gt[el] = connected_elements[el].ghost_type; if ((gt[0] == _not_ghost) xor (gt[1] == _not_ghost)) { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { prank[el] = synchronizer->getRank(connected_elements[el]); } // ugly trick from Marco detected :P bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) facet_ghost_type = _not_ghost; else facet_ghost_type = _ghost; connectivity_facets = &mesh_facets.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_facets.getElementToSubelement( facet_type, facet_ghost_type); } } else { /// facet of facet for (UInt i = 1; i < nb_element_connected_to_facet; ++i) { elements.push_back(connected_elements[i]); } } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->size() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.size(); ++elem) { Element loc_el = elements[elem]; if (loc_el.type == _not_defined) continue; Array & subelement_to_element = mesh_facets.getSubelementToElement(loc_el.type, loc_el.ghost_type); UInt nb_facet_per_loc_element = subelement_to_element.getNbComponent(); for (UInt f_in = 0; f_in < nb_facet_per_loc_element; ++f_in) { auto & el = subelement_to_element(loc_el.element, f_in); if (el.type != _not_defined) continue; el.type = facet_type; el.element = current_facet; el.ghost_type = facet_ghost_type; break; } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = &mesh_accessor.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_accessor.getElementToSubelement( facet_type, facet_ghost_type); } } } } } } // restore the parent of mesh_facet if (mesh_facets_parent) mesh_facets.defineMeshParent(*mesh_facets_parent); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, Array & local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes_numbers) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; for (UInt i = 0; i < old_nodes_numbers.size(); ++i) { renumbering_map[old_nodes_numbers(i)] = i; } /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element) * nb_nodes_per_element, renumbering_map); old_nodes_numbers.resize(renumbering_map.size()); for (auto & renumber_pair : renumbering_map) { old_nodes_numbers(renumber_pair.second) = renumber_pair.first; } renumbering_map.clear(); MeshAccessor mesh_accessor(mesh); /// copy the renumbered connectivity to the right place auto & local_conn = mesh_accessor.getConnectivity(type); local_conn.resize(nb_local_element); - - if(nb_local_element > 0) { + + if (nb_local_element > 0) { memcpy(local_conn.storage(), local_connectivities.storage(), nb_local_element * nb_nodes_per_element * sizeof(UInt)); } auto & ghost_conn = mesh_accessor.getConnectivity(type, _ghost); ghost_conn.resize(nb_ghost_element); - - if(nb_ghost_element > 0) { + + if (nb_ghost_element > 0) { std::memcpy(ghost_conn.storage(), local_connectivities.storage() + - nb_local_element * nb_nodes_per_element, + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); } - + auto & ghost_counter = mesh_accessor.getGhostsCounters(type, _ghost); ghost_counter.resize(nb_ghost_element, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity( Array & list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes.storage(); UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; auto it = renumbering_map.find(node); if (it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Array & nodes_removed = remove_nodes.getList(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(_all_dimensions, ghost_type, _ek_not_defined)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); Array & connectivity = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity.size()); renumberNodesInConnectivity( connectivity, nb_element * nb_nodes_per_element, renumbering_map); } } Array & new_numbering = remove_nodes.getNewNumbering(); std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); for (auto && pair : renumbering_map) { new_numbering(std::get<0>(pair)) = std::get<1>(pair); } for (UInt i = 0; i < new_numbering.size(); ++i) { if (new_numbering(i) == UInt(-1)) nodes_removed.push_back(i); } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_COHESIVE_ELEMENT) /* -------------------------------------------------------------------------- */ UInt MeshUtils::insertCohesiveElements( Mesh & mesh, Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion, Array & doubled_nodes, Array & new_elements, bool only_double_facets) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt elements_to_insert = updateFacetToDouble(mesh_facets, facet_insertion); if (elements_to_insert > 0) { if (spatial_dimension == 1) { doublePointFacet(mesh, mesh_facets, doubled_nodes); } else { doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes, true); findSubfacetToDouble(mesh_facets); if (spatial_dimension == 2) { doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes); } else if (spatial_dimension == 3) { doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false); findSubfacetToDouble(mesh_facets); doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes); } } if (!only_double_facets) updateCohesiveData(mesh, mesh_facets, new_elements); } return elements_to_insert; } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doubleNodes(Mesh & mesh, const std::vector & old_nodes, Array & doubled_nodes) { AKANTU_DEBUG_IN(); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt old_nb_nodes = position.size(); UInt new_nb_nodes = old_nb_nodes + old_nodes.size(); UInt old_nb_doubled_nodes = doubled_nodes.size(); UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.size(); position.resize(new_nb_nodes); doubled_nodes.resize(new_nb_doubled_nodes); Array::iterator> position_begin = position.begin(spatial_dimension); for (UInt n = 0; n < old_nodes.size(); ++n) { UInt new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes[n]; doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node; /// update position std::copy(position_begin + old_nodes[n], position_begin + old_nodes[n] + 1, position_begin + new_node); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets, UInt facet_dimension, Array & doubled_nodes, bool facet_mode) { AKANTU_DEBUG_IN(); + NewElementsEvent event; + for (auto gt_facet : ghost_types) { for (auto && type_facet : mesh_facets.elementTypes(facet_dimension, gt_facet)) { auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); if (nb_facet_to_double == 0) continue; // this while fail if multiple facet types // \TODO handle multiple sub-facet types auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(type_facet); auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto nb_nodes_per_facet = conn_facet.getNbComponent(); auto old_nb_facet = conn_facet.size(); auto new_nb_facet = old_nb_facet + nb_facet_to_double; #ifndef AKANTU_NDEBUG // memory initialization are slow but help debug conn_facet.resize(new_nb_facet, UInt(-1)); #else conn_facet.resize(new_nb_facet); #endif + Element facet{type_facet, 0, gt_facet}; + for (auto el : arange(old_nb_facet, new_nb_facet)) { + facet.element = el; + event.getList().push_back(facet); + } + auto conn_facet_begin = conn_facet.begin(nb_nodes_per_facet); auto & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); #ifndef AKANTU_NDEBUG subfacet_to_facet.resize(new_nb_facet, ElementNull); #else subfacet_to_facet.resize(new_nb_facet); #endif auto subfacet_to_facet_begin = subfacet_to_facet.begin(nb_subfacet_per_facet); Element new_facet{type_facet, old_nb_facet, gt_facet}; auto conn_facet_new_it = conn_facet_begin + new_facet.element; auto subfacet_to_facet_new_it = subfacet_to_facet_begin + new_facet.element; for (UInt facet = 0; facet < nb_facet_to_double; ++facet, ++new_facet.element, ++conn_facet_new_it, ++subfacet_to_facet_new_it) { UInt old_facet = facets_to_double(facet); /// adding a new facet by copying original one /// copy connectivity in new facet *conn_facet_new_it = conn_facet_begin[old_facet]; /// update subfacet_to_facet *subfacet_to_facet_new_it = subfacet_to_facet_begin[old_facet]; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; /// update facet_to_subfacet array mesh_facets.getElementToSubelement(subfacet).push_back(new_facet); } } /// update facet_to_subfacet and _segment_3 facets if any if (not facet_mode) { updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true); updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true); updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } else updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } } + mesh_facets.sendEvent(event); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt MeshUtils::updateFacetToDouble( Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); UInt nb_facets_to_double = 0.; for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { const auto & f_insertion = facet_insertion(type_facet, gt_facet); auto & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); Element old_facet_el{type_facet, 0, gt_facet}; UInt nb_facets = mesh_facets.getNbElement(type_facet, gt_facet); for (UInt f = 0; f < f_insertion.size(); ++f) { if (f_insertion(f) == false) continue; ++nb_facets_to_double; if (element_to_facet(f)[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || element_to_facet(f)[1].kind() == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } f_to_double.push_back(f); UInt new_facet = nb_facets + f_to_double.size() - 1; old_facet_el.element = f; /// update facet_to_element vector auto & elem_to_update = element_to_facet(f)[1]; UInt el = elem_to_update.element; auto & facet_to_element = mesh_facets.getSubelementToElement( elem_to_update.type, elem_to_update.ghost_type); auto el_facets = Vector( make_view(facet_to_element, facet_to_element.getNbComponent()) .begin()[el]); auto f_update = std::find(el_facets.begin(), el_facets.end(), old_facet_el); AKANTU_DEBUG_ASSERT(f_update != el_facets.end(), "Facet not found"); f_update->element = new_facet; /// update elements connected to facet const auto & first_facet_list = element_to_facet(f); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1]; element_to_facet(new_facet)[1] = ElementNull; element_to_facet(f)[1] = ElementNull; } } } AKANTU_DEBUG_OUT(); return nb_facets_to_double; } /* -------------------------------------------------------------------------- */ void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { for (auto type : mesh_facets.elementTypes(_all_dimensions, gt)) { mesh_facets.getDataPointer("facet_to_double", type, gt, 1, false); mesh_facets.getDataPointer>( "facets_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "elements_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "subfacets_to_subsubfacet_double", type, gt, 1, false); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::findSubfacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); if (spatial_dimension == 1) return; for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); GhostType gt_subfacet = _casper; ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet); GhostType gt_subsubfacet = _casper; Array * conn_subfacet = nullptr; Array * sf_to_double = nullptr; Array> * sf_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; Array> * el_to_subfacet_double = nullptr; UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet); UInt nb_subsubfacet; UInt nb_nodes_per_sf_el; if (subsubfacet_mode) { nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet); nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet); } else nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array * subsubfacet_to_subfacet = nullptr; UInt old_nb_facet = subfacet_to_facet.size() - nb_facet_to_double; Element current_facet{type_facet, 0, gt_facet}; std::vector element_list; std::vector facet_list; std::vector * subfacet_list; if (subsubfacet_mode) subfacet_list = new std::vector; /// map to filter subfacets Array> * facet_to_subfacet = nullptr; /// this is used to make sure that both new and old facets are /// checked UInt facets[2]; /// loop on every facet for (UInt f_index = 0; f_index < 2; ++f_index) { for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { facets[bool(f_index)] = facets_to_double(facet); facets[!bool(f_index)] = old_nb_facet + facet; UInt old_facet = facets[0]; UInt new_facet = facets[1]; Element & starting_element = element_to_facet(new_facet)[0]; current_facet.element = old_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; if (subsubfacet_mode) { subsubfacet_to_subfacet = &mesh_facets.getSubelementToElement( type_subfacet, gt_subfacet); } else { conn_subfacet = &mesh_facets.getConnectivity(type_subfacet, gt_subfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subfacet, gt_subfacet); } } if (subsubfacet_mode) { /// loop on every subsubfacet for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) { Element & subsubfacet = (*subsubfacet_to_subfacet)(subfacet.element, ssf); if (subsubfacet == ElementNull) continue; if (gt_subsubfacet != subsubfacet.ghost_type) { gt_subsubfacet = subsubfacet.ghost_type; conn_subfacet = &mesh_facets.getConnectivity(type_subsubfacet, gt_subsubfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subsubfacet, gt_subsubfacet); sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subsubfacet, gt_subsubfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subsubfacet, gt_subsubfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subsubfacet, gt_subsubfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subsubfacet, gt_subsubfacet); } UInt global_ssf = subsubfacet.element; Vector subsubfacet_connectivity( conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet( mesh_facets, starting_element, current_facet, subsubfacet_connectivity, element_list, facet_list, subfacet_list) == false && removeElementsInVector(*subfacet_list, (*facet_to_subfacet)(global_ssf)) == false) { sf_to_double->push_back(global_ssf); sf_to_subfacet_double->push_back(*subfacet_list); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } else { const UInt global_sf = subfacet.element; Vector subfacet_connectivity( conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subfacet is to be doubled if (findElementsAroundSubfacet( mesh_facets, starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false && removeElementsInVector( facet_list, (*facet_to_subfacet)(global_sf)) == false) { sf_to_double->push_back(global_sf); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } } } if (subsubfacet_mode) delete subfacet_list; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets, Array & new_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; MeshAccessor mesh_facets_accessor(mesh_facets); for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); auto & facet_to_coh_element = mesh_facets_accessor.getSubelementToElement(type_cohesive, gt_facet); auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto & conn_cohesive = mesh.getConnectivity(type_cohesive, gt_facet); UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt old_nb_cohesive_elements = conn_cohesive.size(); UInt new_nb_cohesive_elements = conn_cohesive.size() + nb_facet_to_double; UInt old_nb_facet = element_to_facet.size() - nb_facet_to_double; facet_to_coh_element.resize(new_nb_cohesive_elements); conn_cohesive.resize(new_nb_cohesive_elements); UInt new_elements_old_size = new_elements.size(); new_elements.resize(new_elements_old_size + nb_facet_to_double); Element c_element{type_cohesive, 0, gt_facet}; Element f_element{type_facet, 0, gt_facet}; UInt facets[2]; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension ? 1 : 0] = f_to_double(facet); facets[third_dimension ? 0 : 1] = old_nb_facet + facet; UInt cohesive_element = old_nb_cohesive_elements + facet; /// store doubled facets f_element.element = facets[0]; facet_to_coh_element(cohesive_element, 0) = f_element; f_element.element = facets[1]; facet_to_coh_element(cohesive_element, 1) = f_element; /// modify cohesive elements' connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) { conn_cohesive(cohesive_element, n) = conn_facet(facets[0], n); conn_cohesive(cohesive_element, n + nb_nodes_per_facet) = conn_facet(facets[1], n); } /// update element_to_facet vectors c_element.element = cohesive_element; element_to_facet(facets[0])[1] = c_element; element_to_facet(facets[1])[1] = c_element; /// add cohesive element to the element event list new_elements(new_elements_old_size + facet) = c_element; } } } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); if (spatial_dimension != 1) return; + NewElementsEvent event; + auto & position = mesh.getNodes(); for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); const auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); auto new_nb_facet = element_to_facet.size(); auto old_nb_facet = element_to_facet.size() - nb_facet_to_double; auto old_nb_nodes = position.size(); auto new_nb_nodes = old_nb_nodes + nb_facet_to_double; position.resize(new_nb_nodes); conn_facet.resize(new_nb_facet); + Element facet{type_facet, 0, gt_facet}; + for (auto el : arange(old_nb_facet, new_nb_facet)) { + facet.element = el; + event.getList().push_back(facet); + } + auto old_nb_doubled_nodes = doubled_nodes.size(); doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double); for (auto && data_facet : enumerate(facets_to_double)) { const auto & old_facet = std::get<1>(data_facet); auto facet = std::get<0>(data_facet); auto new_facet = old_nb_facet + facet; auto el = element_to_facet(new_facet)[0]; auto old_node = conn_facet(old_facet); auto new_node = old_nb_nodes + facet; /// update position position(new_node) = position(old_node); conn_facet(new_facet) = new_node; Vector conn_segment = mesh.getConnectivity(el); /// update facet connectivity auto it = std::find(conn_segment.begin(), conn_segment.end(), old_node); *it = new_node; doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node; doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node; } } } + mesh_facets.sendEvent(event); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, GhostType gt_facet, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (type_facet != _segment_3) return; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); UInt old_nb_facet = mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); /// this ones matter only for segments in 3D Array> * el_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; if (third_dim_segments) { el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_facet, gt_facet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_facet, gt_facet); } std::vector middle_nodes; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt node = conn_facet(old_facet, 2); if (!mesh.isPureGhostNode(node)) middle_nodes.push_back(node); } UInt n = doubled_nodes.size(); doubleNodes(mesh, middle_nodes, doubled_nodes); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt old_node = conn_facet(old_facet, 2); if (mesh.isPureGhostNode(old_node)) continue; UInt new_node = doubled_nodes(n, 1); UInt new_facet = old_nb_facet + facet; conn_facet(new_facet, 2) = new_node; if (third_dim_segments) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_facet(new_facet)); updateElementalConnectivity(mesh, old_node, new_node, (*el_to_subfacet_double)(facet), &(*f_to_subfacet_double)(facet)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_facet(new_facet)); } ++n; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); /// update subfacet_to_facet vector ElementType type_facet = _not_defined; GhostType gt_facet = _casper; Array * subfacet_to_facet = nullptr; UInt nb_subfacet_per_facet = 0; UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) - nb_subfacet_to_double; Array> * facet_list = nullptr; if (facet_mode) facet_list = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); else facet_list = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); Element old_subfacet_el{type_subfacet, 0, gt_subfacet}; Element new_subfacet_el{type_subfacet, 0, gt_subfacet}; for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { old_subfacet_el.element = sf_to_double(sf); new_subfacet_el.element = old_nb_subfacet + sf; for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) { Element & facet = (*facet_list)(sf)[f]; if (facet.type != type_facet || facet.ghost_type != gt_facet) { type_facet = facet.type; gt_facet = facet.ghost_type; subfacet_to_facet = &mesh_facets.getSubelementToElement(type_facet, gt_facet); nb_subfacet_per_facet = subfacet_to_facet->getNbComponent(); } Element * sf_update = std::find( subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet, subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, old_subfacet_el); AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet != subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, "Subfacet not found"); *sf_update = new_subfacet_el; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); Array> & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); Array> * facet_to_subfacet_double = nullptr; if (facet_mode) { facet_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); } else { facet_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); } UInt old_nb_subfacet = facet_to_subfacet.size(); facet_to_subfacet.resize(old_nb_subfacet + nb_subfacet_to_double); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) facet_to_subfacet(old_nb_subfacet + sf) = (*facet_to_subfacet_double)(sf); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (spatial_dimension == 1) return; + NewElementsEvent event; + for (auto gt_subfacet : ghost_types) { for (auto type_subfacet : mesh_facets.elementTypes(0, gt_subfacet)) { auto & sf_to_double = mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); if (nb_subfacet_to_double == 0) continue; AKANTU_DEBUG_ASSERT( type_subfacet == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); auto & conn_subfacet = mesh_facets.getConnectivity(type_subfacet, gt_subfacet); UInt old_nb_subfacet = conn_subfacet.size(); UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double; conn_subfacet.resize(new_nb_subfacet); + Element subfacet{type_subfacet, 0, gt_subfacet}; + for ( auto el : arange(old_nb_subfacet, new_nb_subfacet)) { + subfacet.element = el; + event.getList().push_back(subfacet); + } + std::vector nodes_to_double; UInt old_nb_doubled_nodes = doubled_nodes.size(); /// double nodes for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_subfacet = sf_to_double(sf); nodes_to_double.push_back(conn_subfacet(old_subfacet)); } doubleNodes(mesh, nodes_to_double, doubled_nodes); /// add new nodes in connectivity for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt new_subfacet = old_nb_subfacet + sf; UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); conn_subfacet(new_subfacet) = new_node; } /// update facet and element connectivity Array> & f_to_subfacet_double = mesh_facets.getData>("facets_to_subfacet_double", type_subfacet, gt_subfacet); Array> & el_to_subfacet_double = mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); Array> * sf_to_subfacet_double = nullptr; if (spatial_dimension == 3) sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0); UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); updateElementalConnectivity(mesh, old_node, new_node, el_to_subfacet_double(sf), &f_to_subfacet_double(sf)); updateElementalConnectivity(mesh_facets, old_node, new_node, f_to_subfacet_double(sf)); if (spatial_dimension == 3) updateElementalConnectivity(mesh_facets, old_node, new_node, (*sf_to_subfacet_double)(sf)); } if (spatial_dimension == 2) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true); } else if (spatial_dimension == 3) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false); } } } + mesh_facets.sendEvent(event); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::flipFacets( Mesh & mesh_facets, const ElementTypeMapArray & remote_global_connectivities, GhostType gt_facet) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); /// get global connectivity for local mesh ElementTypeMapArray local_global_connectivities( "local_global_connectivity", mesh_facets.getID(), mesh_facets.getMemoryID()); local_global_connectivities.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _ghost_type = gt_facet, _with_nb_nodes_per_element = true, _with_nb_element = true); mesh_facets.getGlobalConnectivity(local_global_connectivities); /// loop on every facet for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & connectivity = mesh_facets.getConnectivity(type_facet, gt_facet); auto & local_global_connectivity = local_global_connectivities(type_facet, gt_facet); const auto & remote_global_connectivity = remote_global_connectivities(type_facet, gt_facet); auto & element_per_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); auto & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); auto nb_nodes_per_facet = connectivity.getNbComponent(); auto nb_nodes_per_P1_facet = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet)); for (auto && data : zip(make_view(connectivity, nb_nodes_per_facet), make_view(local_global_connectivity, nb_nodes_per_facet), make_view(remote_global_connectivity, nb_nodes_per_facet), make_view(subfacet_to_facet, subfacet_to_facet.getNbComponent()), make_view(element_per_facet))) { auto & conn = std::get<0>(data); auto & local_gconn = std::get<1>(data); const auto & remote_gconn = std::get<2>(data); /// skip facet if connectivities are the same if (local_gconn == remote_gconn) continue; /// re-arrange connectivity auto conn_tmp = conn; auto begin = local_gconn.begin(); auto end = local_gconn.end(); std::transform(remote_gconn.begin(), remote_gconn.end(), conn.begin(), [&](auto && gnode) { auto it = std::find(begin, end, gnode); AKANTU_DEBUG_ASSERT(it != end, "Node not found"); return conn_tmp(it - begin); }); /// if 3D, check if facets are just rotated if (spatial_dimension == 3) { auto begin = remote_gconn.storage(); /// find first node auto it = std::find(begin, begin + remote_gconn.size(), local_gconn(0)); UInt n, start = it - begin; /// count how many nodes in the received connectivity follow /// the same order of those in the local connectivity for (n = 1; n < nb_nodes_per_P1_facet && local_gconn(n) == remote_gconn((start + n) % nb_nodes_per_P1_facet); ++n) ; /// skip the facet inversion if facet is just rotated if (n == nb_nodes_per_P1_facet) continue; } /// update data to invert facet auto & element_per_facet = std::get<4>(data); std::swap(element_per_facet[0], element_per_facet[1]); auto & subfacets_of_facet = std::get<3>(data); std::swap(subfacets_of_facet(0), subfacets_of_facet(1)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if (mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call " "the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID(), mesh.getMemoryID()); barycenters.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = _all_dimensions); // mesh.initElementTypeMapArray(barycenters, spatial_dimension, // _all_dimensions); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(_all_dimensions, ghost_type)) { element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & barycenters_arr = barycenters(type, ghost_type); barycenters_arr.resize(nb_element); auto bary = barycenters_arr.begin(spatial_dimension); auto bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { element.element = el; mesh.getBarycenter(element, *bary); } } } MeshAccessor mesh_accessor(mesh); for (Int sp(spatial_dimension); sp >= 1; --sp) { if (mesh.getNbElement(sp) == 0) continue; for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(sp, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type) .resize(mesh.getNbElement(type, ghost_type)); mesh_accessor.getSubelementToElement(type, ghost_type).set(ElementNull); } for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { mesh_accessor.getElementToSubelement(type, ghost_type) .resize(mesh.getNbElement(type, ghost_type)); mesh.getElementToSubelement(type, ghost_type).clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (auto ghost_type : ghost_types) { facet_element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { facet_element.type = type; auto & element_to_subelement = mesh.getElementToSubelement(type, ghost_type); const auto & connectivity = mesh.getConnectivity(type, ghost_type); for (auto && data : enumerate( make_view(connectivity, mesh.getNbNodesPerElement(type)))) { const auto & facet = std::get<1>(data); facet_element.element = std::get<0>(data); std::map element_seen_counter; auto nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(type)); // count the number of node in common between the facet and the other // element connected to the nodes of the facet for (auto node : arange(nb_nodes_per_facet)) { for (auto & elem : nodes_to_elements.getRow(facet(node))) { auto cit = element_seen_counter.find(elem); if (cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } // check which are the connected elements std::vector connected_elements; for (auto && cit : element_seen_counter) { if (cit.second == nb_nodes_per_facet) connected_elements.push_back(cit.first); } // add the connected elements as sub-elements for (auto & connected_element : connected_elements) { element_to_subelement(facet_element.element) .push_back(connected_element); } // add the element as sub-element to the connected elements for (auto & connected_element : connected_elements) { Vector subelements_to_element = mesh.getSubelementToElement(connected_element); // find the position where to insert the element auto it = std::find(subelements_to_element.begin(), subelements_to_element.end(), ElementNull); AKANTU_DEBUG_ASSERT( it != subelements_to_element.end(), "The element " << connected_element << " seems to have too many facets!! (" << (it - subelements_to_element.begin()) << " < " << mesh.getNbFacetsPerElement(connected_element.type) << ")"); *it = facet_element; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MeshUtils::findElementsAroundSubfacet( const Mesh & mesh_facets, const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & element_list, std::vector & facet_list, std::vector * subfacet_list) { AKANTU_DEBUG_IN(); bool facet_matched = false; element_list.clear(); facet_list.clear(); if (third_dim_points) { subfacet_list->clear(); } element_list.push_back(starting_element); std::queue elements_to_check; elements_to_check.push(starting_element); /// keep going as long as there are elements to check while (not elements_to_check.empty()) { /// check current element Element & current_element = elements_to_check.front(); const Vector facets_to_element = mesh_facets.getSubelementToElement(current_element); // for every facet of the element for (auto & current_facet : facets_to_element) { if (current_facet == ElementNull) continue; if (current_facet == end_facet) facet_matched = true; // facet already listed if (std::find(facet_list.begin(), facet_list.end(), current_facet) != facet_list.end()) continue; // subfacet_connectivity is not in the connectivity of current_facet; if ((std::find(facet_list.begin(), facet_list.end(), current_facet) != facet_list.end()) or not hasElement(mesh_facets.getConnectivity(current_facet), subfacet_connectivity)) continue; facet_list.push_back(current_facet); if (third_dim_points) { const Vector subfacets_of_facet = mesh_facets.getSubelementToElement(current_facet); /// check subfacets for (const auto & current_subfacet : subfacets_of_facet) { if (current_subfacet == ElementNull) continue; if ((std::find(subfacet_list->begin(), subfacet_list->end(), current_subfacet) == subfacet_list->end()) and hasElement(mesh_facets.getConnectivity(current_subfacet), subfacet_connectivity)) subfacet_list->push_back(current_subfacet); } } /// consider opposing element const auto & elements_to_facet = mesh_facets.getElementToSubelement(current_facet); UInt opposing = 0; if (elements_to_facet[0] == current_element) opposing = 1; auto & opposing_element = elements_to_facet[opposing]; /// skip null elements since they are on a boundary if (opposing_element == ElementNull) continue; /// skip this element if already added if (std::find(element_list.begin(), element_list.end(), opposing_element) != element_list.end()) continue; /// only regular elements have to be checked if (opposing_element.kind() == _ek_regular) elements_to_check.push(opposing_element); element_list.push_back(opposing_element); AKANTU_DEBUG_ASSERT( hasElement( mesh_facets.getMeshParent().getConnectivity(opposing_element), subfacet_connectivity), "Subfacet doesn't belong to this element"); } /// erased checked element from the list elements_to_check.pop(); } AKANTU_DEBUG_OUT(); return facet_matched; } /* -------------------------------------------------------------------------- */ void MeshUtils::updateElementalConnectivity( Mesh & mesh, UInt old_node, UInt new_node, const std::vector & element_list, const std::vector * #if defined(AKANTU_COHESIVE_ELEMENT) facet_list #endif ) { AKANTU_DEBUG_IN(); for (auto & element : element_list) { if (element.type == _not_defined) continue; Vector connectivity = mesh.getConnectivity(element); #if defined(AKANTU_COHESIVE_ELEMENT) if (element.kind() == _ek_cohesive) { AKANTU_DEBUG_ASSERT( facet_list != nullptr, "Provide a facet list in order to update cohesive elements"); const Vector facets = mesh.getMeshFacets().getSubelementToElement(element); auto facet_nb_nodes = connectivity.size() / 2; /// loop over cohesive element's facets for (const auto & facet : enumerate(facets)) { /// skip facets if not present in the list if (std::find(facet_list->begin(), facet_list->end(), std::get<1>(facet)) == facet_list->end()) { continue; } auto n = std::get<0>(facet); auto begin = connectivity.begin() + n * facet_nb_nodes; auto end = begin + facet_nb_nodes; auto it = std::find(begin, end, old_node); AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element"); *it = new_node; } } else #endif { auto it = std::find(connectivity.begin(), connectivity.end(), old_node); AKANTU_DEBUG_ASSERT(it != connectivity.end(), "Node not found in current element"); /// update connectivity *it = new_node; } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/common/boundary_condition/boundary_condition_tmpl.hh b/src/model/common/boundary_condition/boundary_condition_tmpl.hh index ebb81e561..c64b4b44e 100644 --- a/src/model/common/boundary_condition/boundary_condition_tmpl.hh +++ b/src/model/common/boundary_condition/boundary_condition_tmpl.hh @@ -1,230 +1,230 @@ /** * @file boundary_condition_tmpl.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Tue Feb 20 2018 * * @brief implementation of the applyBC * * @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.hh" #include "element_group.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_BOUNDARY_CONDITION_TMPL_HH__ #define __AKANTU_BOUNDARY_CONDITION_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template void BoundaryCondition::initBC(ModelType & model, Array & primal, Array & dual) { this->model = &model; this->primal = &primal; this->dual = &dual; } /* -------------------------------------------------------------------------- */ template void BoundaryCondition::initBC(ModelType & model, Array & primal, Array & primal_increment, Array & dual) { this->initBC(model, primal, dual); this->primal_increment = &primal_increment; } /* -------------------------------------------------------------------------- */ /* Partial specialization for DIRICHLET functors */ template template struct BoundaryCondition::TemplateFunctionWrapper< FunctorType, BC::Functor::_dirichlet> { static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance) { auto & model = bc_instance.getModel(); auto & primal = bc_instance.getPrimal(); const auto & coords = model.getMesh().getNodes(); auto & boundary_flags = model.getBlockedDOFs(); UInt dim = model.getMesh().getSpatialDimension(); auto primal_iter = primal.begin(primal.getNbComponent()); auto coords_iter = coords.begin(dim); auto flags_iter = boundary_flags.begin(boundary_flags.getNbComponent()); for (auto n : group.getNodeGroup()) { Vector flag(flags_iter[n]); Vector primal(primal_iter[n]); Vector coords(coords_iter[n]); func(n, flag, primal, coords); } } }; /* -------------------------------------------------------------------------- */ /* Partial specialization for NEUMANN functors */ template template struct BoundaryCondition::TemplateFunctionWrapper< FunctorType, BC::Functor::_neumann> { static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance) { UInt dim = bc_instance.getModel().getSpatialDimension(); switch (dim) { case 1: { AKANTU_TO_IMPLEMENT(); break; } case 2: case 3: { applyBC(func, group, bc_instance, _not_ghost); applyBC(func, group, bc_instance, _ghost); break; } } } static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance, GhostType ghost_type) { auto & model = bc_instance.getModel(); auto & dual = bc_instance.getDual(); const auto & mesh = model.getMesh(); const auto & nodes_coords = mesh.getNodes(); const auto & fem_boundary = model.getFEEngineBoundary(); - + UInt dim = model.getSpatialDimension(); UInt nb_degree_of_freedom = dual.getNbComponent(); IntegrationPoint quad_point; quad_point.ghost_type = ghost_type; // Loop over the boundary element types for (auto && type : group.elementTypes(dim - 1, ghost_type)) { const auto & element_ids = group.getElements(type, ghost_type); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type, ghost_type); UInt nb_elements = element_ids.size(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array dual_before_integ(nb_elements * nb_quad_points, nb_degree_of_freedom, 0.); Array quad_coords(nb_elements * nb_quad_points, dim); const auto & normals_on_quad = fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type); fem_boundary.interpolateOnIntegrationPoints( nodes_coords, quad_coords, dim, type, ghost_type, element_ids); auto normals_begin = normals_on_quad.begin(dim); decltype(normals_begin) normals_iter; auto quad_coords_iter = quad_coords.begin(dim); auto dual_iter = dual_before_integ.begin(nb_degree_of_freedom); quad_point.type = type; for (auto el : element_ids) { quad_point.element = el; normals_iter = normals_begin + el * nb_quad_points; for (auto && q : arange(nb_quad_points)) { quad_point.num_point = q; func(quad_point, *dual_iter, *quad_coords_iter, *normals_iter); ++dual_iter; ++quad_coords_iter; ++normals_iter; } } Array dual_by_shapes(nb_elements * nb_quad_points, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type, ghost_type, element_ids); Array dual_by_shapes_integ( nb_elements, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.integrate(dual_by_shapes, dual_by_shapes_integ, nb_degree_of_freedom * nb_nodes_per_element, type, ghost_type, element_ids); // assemble the result into force vector model.getDOFManager().assembleElementalArrayLocalArray( dual_by_shapes_integ, dual, type, ghost_type, 1., element_ids); } } }; /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func) { auto bit = model->getMesh().getGroupManager().element_group_begin(); auto bend = model->getMesh().getGroupManager().element_group_end(); for (; bit != bend; ++bit) applyBC(func, *bit); } /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func, const std::string & group_name) { try { const ElementGroup & element_group = model->getMesh().getElementGroup(group_name); applyBC(func, element_group); } catch (akantu::debug::Exception & e) { AKANTU_EXCEPTION("Error applying a boundary condition onto \"" << group_name << "\"! [" << e.what() << "]"); } } /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func, const ElementGroup & element_group) { #if !defined(AKANTU_NDEBUG) if (element_group.getDimension() != model->getSpatialDimension() - 1) AKANTU_DEBUG_WARNING("The group " << element_group.getName() << " does not contain only boundaries elements"); #endif TemplateFunctionWrapper::applyBC(func, element_group, *this); } #endif /* __AKANTU_BOUNDARY_CONDITION_TMPL_HH__ */ } // namespace akantu diff --git a/src/model/common/integration_scheme/integration_scheme.cc b/src/model/common/integration_scheme/integration_scheme.cc index 2b5075a56..a42a9f2f8 100644 --- a/src/model/common/integration_scheme/integration_scheme.cc +++ b/src/model/common/integration_scheme/integration_scheme.cc @@ -1,76 +1,77 @@ /** * @file integration_scheme.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Jan 31 2018 * * @brief Common interface to all interface schemes * * @section LICENSE * * Copyright (©) 2015-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 "integration_scheme.hh" #include "dof_manager.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ IntegrationScheme::IntegrationScheme(DOFManager & dof_manager, const ID & dof_id, UInt order) : Parsable(ParserType::_integration_scheme, dof_id), dof_manager(dof_manager), dof_id(dof_id), order(order) {} /* -------------------------------------------------------------------------- */ /// standard input stream operator for SolutionType std::istream & operator>>(std::istream & stream, IntegrationScheme::SolutionType & type) { std::string str; stream >> str; if (str == "displacement") type = IntegrationScheme::_displacement; else if (str == "temperature") type = IntegrationScheme::_temperature; else if (str == "velocity") type = IntegrationScheme::_velocity; else if (str == "temperature_rate") type = IntegrationScheme::_temperature_rate; else if (str == "acceleration") type = IntegrationScheme::_acceleration; else { stream.setstate(std::ios::failbit); } return stream; } -// void IntegrationScheme::assembleJacobian(const SolutionType & /*type*/, Real /*delta_t*/) { +// void IntegrationScheme::assembleJacobian(const SolutionType & /*type*/, Real +// /*delta_t*/) { // auto & J = dof_manager.getLumpedMatrix("J"); // auto & M = dof_manager.getLumpedMatrix("M"); // if (J.release() == M.release()) return; // J = M; // } -} // akantu +} // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 5f9ba3174..0dd47c922 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,480 +1,481 @@ /** * @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, const MemoryID & memory_id) - : ContactDetector(mesh, mesh.getNodes(), id, memory_id) { - - } - +ContactDetector::ContactDetector(Mesh & mesh, const ID & id, + const MemoryID & memory_id) + : ContactDetector(mesh, mesh.getNodes(), id, memory_id) {} + /* -------------------------------------------------------------------------- */ -ContactDetector::ContactDetector(Mesh & mesh, Array positions, const ID & id, - const MemoryID & memory_id) - : Memory(id, memory_id), - Parsable(ParserType::_contact_detector, id), - mesh(mesh) { +ContactDetector::ContactDetector(Mesh & mesh, Array positions, + const ID & id, const MemoryID & memory_id) + : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), + mesh(mesh), positions(0, mesh.getSpatialDimension()) { - AKANTU_DEBUG_IN(); + AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); - this->positions = positions; - - this->mesh.fillNodesToElements(this->spatial_dimension - 1); + this->positions.copy(positions); this->parseSection(); - this->computeMaximalDetectionDistance(); - + // this->mesh.fillNodesToElements(this->spatial_dimension - 1); + + // this->computeMaximalDetectionDistance(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactDetector::parseSection() { const Parser & parser = getStaticParser(); const ParserSection & section = - *(parser.getSubSections(ParserType::_contact_detector).first); + *(parser.getSubSections(ParserType::_contact_detector).first); auto type = section.getParameterValue("type"); - + if (type == "implicit") { this->detection_type = _implicit; - } - else if (type == "explicit"){ + } else if (type == "explicit") { this->detection_type = _explicit; - } - else { + } else { AKANTU_ERROR("Unknown detection type : " << type); } - - surfaces[Surface::master] = section.getParameterValue("master"); - surfaces[Surface::slave ] = section.getParameterValue("slave"); + + // surfaces[Surface::master] = + // section.getParameterValue("master"); surfaces[Surface::slave] = + // section.getParameterValue("slave"); two_pass_algorithm = section.getParameterValue("two_pass_algorithm"); } - + /* -------------------------------------------------------------------------- */ void ContactDetector::search(std::map & contact_map) { + this->mesh.fillNodesToElements(this->spatial_dimension - 1); + this->computeMaximalDetectionDistance(); + contact_pairs.clear(); - + SpatialGrid master_grid(spatial_dimension); SpatialGrid slave_grid(spatial_dimension); - + this->globalSearch(slave_grid, master_grid); - + this->localSearch(slave_grid, master_grid); - //if (two_pass_algorithm) { + // if (two_pass_algorithm) { // this->localSearch(master_grid, slave_grid); //} - + this->constructContactMap(contact_map); } - + /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(SpatialGrid & slave_grid, - SpatialGrid & master_grid) { - - auto & master_list = - mesh.getElementGroup(surfaces[Surface::master]).getNodeGroup().getNodes(); - - auto & slave_list = - mesh.getElementGroup(surfaces[Surface::slave]).getNodeGroup().getNodes(); - + SpatialGrid & master_grid) { + + auto & master_list = node_selector->getMasterList(); + auto & slave_list = node_selector->getSlaveList(); + 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 ); - + + 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); - + master_grid.setCenter(center); master_grid.setSpacing(spacing); this->constructGrid(master_grid, bbox_intersection, master_list); slave_grid.setCenter(center); slave_grid.setSpacing(spacing); this->constructGrid(slave_grid, bbox_intersection, slave_list); - + // 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 elements 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. } /* -------------------------------------------------------------------------- */ void ContactDetector::localSearch(SpatialGrid & slave_grid, - SpatialGrid & master_grid) { + SpatialGrid & master_grid) { // 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 - + // elements + /// find the closet master node for each slave node for (auto && cell_id : slave_grid) { - /// loop over all the slave nodes of the current cell - for (auto && slave_node: slave_grid.getCell(cell_id)) { + /// loop over all the slave nodes of the current cell + for (auto && slave_node : slave_grid.getCell(cell_id)) { bool pair_exists = false; - + Vector pos(spatial_dimension); - for (UInt s: arange(spatial_dimension)) - pos(s) = this->positions(slave_node, s); - + for (UInt s : arange(spatial_dimension)) + pos(s) = this->positions(slave_node, s); + Real closet_distance = std::numeric_limits::max(); UInt closet_master_node; /// loop over all the neighboring cells of the current cell for (auto && neighbor_cell : cell_id.neighbors()) { - /// loop over the data of neighboring cells from master grid - for (auto && master_node : master_grid.getCell(neighbor_cell)) { - - if (slave_node == master_node) - continue; - - Vector pos2(spatial_dimension); - for (UInt s: arange(spatial_dimension)) - pos2(s) = this->positions(master_node, s); - - Real distance = pos.distance(pos2); - - if (distance <= closet_distance) { - closet_master_node = master_node; - closet_distance = distance; - pair_exists = true; - } - } + /// loop over the data of neighboring cells from master grid + for (auto && master_node : master_grid.getCell(neighbor_cell)) { + + /// check for self contact + if (slave_node == master_node) + continue; + + Vector pos2(spatial_dimension); + for (UInt s : arange(spatial_dimension)) + pos2(s) = this->positions(master_node, s); + + Real distance = pos.distance(pos2); + + if (distance <= closet_distance) { + closet_master_node = master_node; + closet_distance = distance; + pair_exists = true; + } + } } - if (pair_exists) - contact_pairs.push_back( std::make_pair(slave_node, - closet_master_node)); + if (pair_exists) + contact_pairs.push_back(std::make_pair(slave_node, closet_master_node)); } - } + } } /* -------------------------------------------------------------------------- */ -void ContactDetector::constructContactMap(std::map & contact_map) { +void ContactDetector::constructContactMap( + std::map & contact_map) { auto surface_dimension = spatial_dimension - 1; std::map previous_contact_map = contact_map; contact_map.clear(); - - auto get_index = [&](auto & gaps, auto & projections) { + auto get_index = [&](auto & gaps, auto & projections) { UInt index; - Real gap_min = std::numeric_limits::max(); + Real gap_min = std::numeric_limits::max(); UInt counter = 0; - for (auto && values : zip(gaps, - make_view(projections, surface_dimension))) { - auto & gap = std::get<0>(values); + for (auto && values : + zip(gaps, make_view(projections, surface_dimension))) { + auto & gap = std::get<0>(values); auto & projection = std::get<1>(values); /// todo adhoc fix to assign a master element in case the /// projection does not lie in the extended element. As it is /// tolerance based bool is_valid = this->checkValidityOfProjection(projection); - + if (is_valid and gap <= gap_min) { - gap_min = gap; - index = counter; + gap_min = gap; + index = counter; } counter++; } if (index >= gaps.size()) { auto gap_min_it = std::min_element(gaps.begin(), gaps.end()); auto index_it = std::find(gaps.begin(), gaps.end(), *gap_min_it); - index = *index_it; + index = *index_it; } - + return index; }; - auto get_connectivity = [&](auto & slave, auto & master) { - Vector master_conn = const_cast(this->mesh).getConnectivity(master); + Vector master_conn = + const_cast(this->mesh).getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { - elem_conn[i] = master_conn[i-1]; - } + elem_conn[i] = master_conn[i - 1]; + } return elem_conn; }; for (auto & pairs : contact_pairs) { - const auto & slave_node = pairs.first; + const auto & slave_node = pairs.first; const auto & master_node = pairs.second; - Array elements; - this->mesh.getAssociatedElements(master_node, elements); + Array all_elements; + this->mesh.getAssociatedElements(master_node, all_elements); - Array gaps(elements.size(), 1, "gaps"); - Array normals(elements.size(), spatial_dimension, "normals"); - Array projections(elements.size(), surface_dimension, "projections"); + Array elements; + this->filterBoundaryElements(all_elements, elements); - this->computeOrthogonalProjection(slave_node, elements, - normals, gaps, projections); + std::cerr << "---" << slave_node << std::endl; + for (auto sl_el : elements) { + std::cerr << sl_el << std::endl; + } + Array gaps(elements.size(), 1, "gaps"); + Array normals(elements.size(), spatial_dimension, "normals"); + Array projections(elements.size(), surface_dimension, "projections"); + + this->computeOrthogonalProjection(slave_node, elements, normals, gaps, + projections); + auto index = get_index(gaps, projections); auto connectivity = get_connectivity(slave_node, elements[index]); - + contact_map[slave_node].setMaster(elements[index]); contact_map[slave_node].setGap(gaps[index]); - contact_map[slave_node].setNormal(Vector(normals.begin(spatial_dimension)[index], true)); - contact_map[slave_node].setProjection(Vector(projections.begin(surface_dimension)[index], true)); + contact_map[slave_node].setNormal( + Vector(normals.begin(spatial_dimension)[index], true)); + contact_map[slave_node].setProjection( + Vector(projections.begin(surface_dimension)[index], true)); contact_map[slave_node].setConnectivity(connectivity); - + Matrix tangents(surface_dimension, spatial_dimension); this->computeTangentsOnElement(contact_map[slave_node].master, - contact_map[slave_node].projection, - tangents); + contact_map[slave_node].projection, + tangents); + + std::cerr << contact_map[slave_node].normal < e_z(3); e_z[0] = 0.; e_z[1] = 0.; e_z[2] = 1.; Vector tangent(3); - tangent[0] = tangents(0, 0); - tangent[1] = tangents(0, 1); - tangent[2] = 0.; - + tangent[0] = tangents(0, 0); + tangent[1] = tangents(0, 1); + tangent[2] = 0.; + auto exp_normal = e_z.crossProduct(tangent); - + auto & cal_normal = contact_map[slave_node].normal; auto ddot = cal_normal.dot(exp_normal); if (ddot < 0) { - tangents *= -1.0; + tangents *= -1.0; } break; } case 3: { auto tang_trans = tangents.transpose(); auto tang1 = Vector(tang_trans(0)); auto tang2 = Vector(tang_trans(1)); auto tang1_cross_tang2 = tang1.crossProduct(tang2); auto exp_normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); - + auto & cal_normal = contact_map[slave_node].normal; auto ddot = cal_normal.dot(exp_normal); if (ddot < 0) { - tang_trans(1) *= -1.0; + tang_trans(1) *= -1.0; } tangents = tang_trans.transpose(); break; } default: break; } - + contact_map[slave_node].setTangent(tangents); auto search = previous_contact_map.find(slave_node); if (search != previous_contact_map.end()) { - auto previous_projection = previous_contact_map[slave_node].getPreviousProjection(); + auto previous_projection = + previous_contact_map[slave_node].getPreviousProjection(); contact_map[slave_node].setPreviousProjection(previous_projection); auto previous_traction = previous_contact_map[slave_node].getTraction(); contact_map[slave_node].setTraction(previous_traction); - } - else { + } else { Vector previous_projection(surface_dimension, 0.); contact_map[slave_node].setPreviousProjection(previous_projection); Vector previous_traction(surface_dimension, 0.); contact_map[slave_node].setTraction(previous_traction); } } contact_pairs.clear(); } - + /* -------------------------------------------------------------------------- */ -void ContactDetector::computeOrthogonalProjection(const UInt & node, - const Array & elements, - Array & normals, Array & gaps, - Array & projections) { +void ContactDetector::computeOrthogonalProjection( + const UInt & node, const Array & elements, Array & normals, + Array & gaps, Array & projections) { Vector query(spatial_dimension); - for (UInt s: arange(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(projections, spatial_dimension - 1))) { - + zip(elements, gaps, make_view(normals, spatial_dimension), + make_view(projections, spatial_dimension - 1))) { + const auto & element = std::get<0>(values); - auto & gap = std::get<1>(values); - auto & normal = std::get<2>(values); - auto & projection = std::get<3>(values); - + auto & gap = std::get<1>(values); + auto & normal = std::get<2>(values); + auto & projection = std::get<3>(values); + this->computeNormalOnElement(element, normal); - + Vector real_projection(spatial_dimension); - this->computeProjectionOnElement(element, normal, query, - projection, real_projection); + this->computeProjectionOnElement(element, normal, query, projection, + real_projection); gap = this->computeGap(query, real_projection, normal); } - } - /* -------------------------------------------------------------------------- */ -void ContactDetector::computeProjectionOnElement(const Element & element, - const Vector & normal, - const Vector & query, - Vector & natural_projection, - Vector & real_projection) { - +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)); + Vector point(coords(0)); Real alpha = (query - point).dot(normal); real_projection = query - alpha * normal; - - this->computeNaturalProjection(element, real_projection, natural_projection); + + this->computeNaturalProjection(element, real_projection, natural_projection); } /* -------------------------------------------------------------------------- */ -void ContactDetector::computeNaturalProjection(const Element & element, Vector & real_projection, - Vector & natural_projection) { +void ContactDetector::computeNaturalProjection( + 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(); + UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(this->positions, 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) + 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 } - - /* -------------------------------------------------------------------------- */ -void ContactDetector::computeTangentsOnElement(const Element & el, Vector & projection, Matrix & tangents) { +void ContactDetector::computeTangentsOnElement(const Element & el, + Vector & projection, + Matrix & tangents) { + + const ElementType & type = el.type; - const ElementType & type = el.type; - UInt nb_nodes_master = Mesh::getNbNodesPerElement(type); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); - -#define GET_SHAPES_NATURAL(type) \ + +#define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); -#undef GET_SHAPES_NATURAL +#undef GET_SHAPES_NATURAL -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ +#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 coords(spatial_dimension, nb_nodes_master); coordinatesOfElement(el, coords); - + tangents.mul(shapes_derivatives, coords); auto temp_tangents = tangents.transpose(); - for (UInt i = 0; i < spatial_dimension -1; ++i) { + for (UInt i = 0; i < spatial_dimension - 1; ++i) { auto temp = Vector(temp_tangents(i)); temp_tangents(i) = temp.normalize(); } tangents = temp_tangents.transpose(); } - - -} // akantu + +} // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.hh b/src/model/contact_mechanics/contact_detector.hh index 833f6b183..2ca8dcce2 100644 --- a/src/model/contact_mechanics/contact_detector.hh +++ b/src/model/contact_mechanics/contact_detector.hh @@ -1,211 +1,210 @@ /** * @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 "aka_grid_dynamic.hh" +#include "aka_memory.hh" +#include "contact_element.hh" +#include "element_class.hh" +#include "element_group.hh" +#include "fe_engine.hh" #include "mesh.hh" #include "mesh_io.hh" -#include "fe_engine.hh" +#include "node_selector.hh" #include "parsable.hh" -#include "element_group.hh" -#include "contact_element.hh" -#include "element_class.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_HH__ #define __AKANTU_CONTACT_DETECTOR_HH__ - namespace akantu { -enum class Surface { - master, - slave -}; - -/* -------------------------------------------------------------------------- */ +enum class Surface { master, slave }; -class ContactDetector : - private Memory, public Parsable { +/* -------------------------------------------------------------------------- */ + +class ContactDetector : private Memory, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: - ContactDetector(Mesh &, const ID & id = "contact_detector", - const MemoryID & memory_id = 0); + const MemoryID & memory_id = 0); + + ContactDetector(Mesh &, Array positions, + const ID & id = "contact_detector", + const MemoryID & memory_id = 0); - ContactDetector(Mesh &, Array positions, const ID & id = "contact_detector", - const MemoryID & memory_id = 0); - ~ContactDetector() = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: - /// performs all search steps + /// performs all search steps void search(std::map &); /// performs global spatial search to construct spatial grids void globalSearch(SpatialGrid &, SpatialGrid &); /// performs local search to find closet master node to a slave node void localSearch(SpatialGrid &, SpatialGrid &); /// constructs contact map for a given pair of slave and master node void constructContactMap(std::map &); - + private: /// reads the input file to get contact detection options void parseSection(); /// computes the orthogonal projection on master elements - void computeOrthogonalProjection(const UInt & , const Array &, - Array &, Array &, Array &); - + void computeOrthogonalProjection(const UInt &, const Array &, + Array &, Array &, Array &); + /// computes tangents on a given natural coordinate - void computeTangentsOnElement(const Element &, Vector &, Matrix &); + void computeTangentsOnElement(const Element &, Vector &, + Matrix &); /// computes projection of a query point on an element void computeProjectionOnElement(const Element &, const Vector &, - const Vector &, Vector &, Vector &); + const Vector &, Vector &, + Vector &); /// computes natural projection of a real projection - void computeNaturalProjection(const Element &, Vector &, Vector &); + void computeNaturalProjection(const Element &, Vector &, + Vector &); /* ------------------------------------------------------------------------ */ /* Inline Methods */ /* ------------------------------------------------------------------------ */ public: /// checks whether the natural projection is valid or not - inline bool checkValidityOfProjection(Vector & ); - + inline bool checkValidityOfProjection(Vector &); + /// extracts the coordinates of an element - inline void coordinatesOfElement(const Element & , Matrix & ); + inline void coordinatesOfElement(const Element &, Matrix &); - /// computes the optimal cell size for grid - inline void computeCellSpacing(Vector & ); + /// computes the optimal cell size for grid + inline void computeCellSpacing(Vector &); /// constructs a grid containing nodes lying within bounding box inline void constructGrid(SpatialGrid &, BBox &, const Array &); /// constructs the bounding box based on nodes list inline void constructBoundingBox(BBox &, const Array &); /// get the surface id - template - inline std::string getSurfaceId(); - + template inline std::string getSurfaceId(); + /// set the surface id - template - inline void setSurfaceId(const std::string); - - /// computes the maximum in radius for a given mesh + template inline void setSurfaceId(const std::string); + + /// computes the maximum in radius for a given mesh inline void computeMaximalDetectionDistance(); /// constructs the connectivity for a contact element inline Vector constructConnectivity(UInt &, const Element &); /// computes normal on an element - inline void computeNormalOnElement(const Element &, Vector & ); - + inline void computeNormalOnElement(const Element &, Vector &); + /// extracts vectors which forms the plane of element - inline void vectorsAlongElement(const Element &, Matrix & ); + inline void vectorsAlongElement(const Element &, Matrix &); /// computes the gap between slave and its projection on master /// surface inline Real computeGap(Vector &, Vector &, Vector &); - + + /// filter boundary elements + inline void filterBoundaryElements(Array & elements, + Array & boundary_elements); + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// returns the maximum detection distance AKANTU_GET_MACRO(MaximumDetectionDistance, max_dd, Real); - - /// sets the maximum detection distance AKANTU_SET_MACRO(MaximumDetectionDistance, max_dd, Real); /// returns the bounding box extension AKANTU_GET_MACRO(MaximumBoundingBox, max_bb, Real); - - /// sets the bounding box extension AKANTU_SET_MACRO(MaximumBoundingBox, max_bb, Real); - /// returns the positions - AKANTU_GET_MACRO(Positions, positions, Array); - - /// sets the positions + AKANTU_GET_MACRO_NOT_CONST(Positions, positions, Array &); AKANTU_SET_MACRO(Positions, positions, Array); - + AKANTU_GET_MACRO_NOT_CONST(NodeSelector, *node_selector, NodeSelector &); + AKANTU_SET_MACRO(NodeSelector, node_selector, std::shared_ptr); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// map to contain ids for surfaces - std::map surfaces; - + std::map surfaces; + 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}; - + + /// node selector for selecting master and slave nodes + std::shared_ptr node_selector; + /// contact pair slave node to closet master node - std::vector > contact_pairs; - + std::vector> contact_pairs; + /// contains the updated positions of the nodes Array positions; - /// type of detection explicit/implicit + /// type of detection explicit/implicit DetectionType detection_type; /// bool two_pass_algorithm; }; - + } // namespace akantu #include "contact_detector_inline_impl.cc" #endif /* __AKANTU_CONTACT_DETECTOR_HH__ */ diff --git a/src/model/contact_mechanics/contact_detector_inline_impl.cc b/src/model/contact_mechanics/contact_detector_inline_impl.cc index 74d6c69b5..0bd81b5f6 100644 --- a/src/model/contact_mechanics/contact_detector_inline_impl.cc +++ b/src/model/contact_mechanics/contact_detector_inline_impl.cc @@ -1,268 +1,325 @@ /** * @file contact_detection.hh * * @author Mohit Pundir * * @date creation: Mon Apr 29 2019 * @date last modification: Mon Apr 29 2019 * * @brief inine implementation of the contact detector 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 "contact_detector.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__ #define __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ -inline bool ContactDetector::checkValidityOfProjection(Vector & projection) { +inline bool +ContactDetector::checkValidityOfProjection(Vector & projection) { UInt nb_xi_inside = 0; Real epsilon = 1e-3; for (auto xi : projection) { - if (xi >= -1.0 - epsilon and xi <= 1.0 + epsilon) + if (xi >= -1.0 - epsilon and xi <= 1.0 + epsilon) nb_xi_inside++; - } + } - if (nb_xi_inside == projection.size()) + if (nb_xi_inside == projection.size()) return true; - + return false; } /* -------------------------------------------------------------------------- */ -inline void ContactDetector::coordinatesOfElement(const Element & el, Matrix & coords) { +inline 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]; + .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)) { + for (UInt s : arange(spatial_dimension)) { coords(s, n) = this->positions(node, s); } } } - + /* -------------------------------------------------------------------------- */ inline void ContactDetector::computeCellSpacing(Vector & spacing) { - for (UInt s: arange(spatial_dimension)) + for (UInt s : arange(spatial_dimension)) spacing(s) = std::sqrt(2.0) * max_dd; } -/* -------------------------------------------------------------------------- */ -inline void ContactDetector::constructBoundingBox(BBox & bbox, const Array & nodes_list) { - +/* -------------------------------------------------------------------------- */ +inline 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); + 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); auto & lower_bound = bbox.getLowerBounds(); auto & upper_bound = bbox.getUpperBounds(); - for (UInt s: arange(spatial_dimension)) { + for (UInt s : arange(spatial_dimension)) { lower_bound(s) -= this->max_bb; upper_bound(s) += this->max_bb; } - + AKANTU_DEBUG_INFO("BBox" << bbox); } - /* -------------------------------------------------------------------------- */ -inline void ContactDetector::constructGrid(SpatialGrid & grid, BBox & bbox, - const Array & nodes_list) { +inline 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)) { + 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); } - + /* -------------------------------------------------------------------------- */ -template -inline std::string ContactDetector::getSurfaceId() { +template inline std::string ContactDetector::getSurfaceId() { return surfaces[id]; } /* -------------------------------------------------------------------------- */ -template +template inline void ContactDetector::setSurfaceId(const std::string name) { surfaces[id] = name; } - + /* -------------------------------------------------------------------------- */ inline void ContactDetector::computeMaximalDetectionDistance() { AKANTU_DEBUG_IN(); Real elem_size; Real max_elem_size = std::numeric_limits::min(); - auto & master_group = - mesh.getElementGroup(surfaces[Surface::master]); + /*auto & master_group = + mesh.getElementGroup(surfaces[Surface::master]); + + for (auto type : master_group.elementTypes(spatial_dimension - 1, _not_ghost, + _ek_regular)) { - for (auto type: - master_group.elementTypes(spatial_dimension - 1, _not_ghost, _ek_regular)) { - - const auto & element_ids = master_group.getElements(type); + const auto & element_ids = master_group.getElements(type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - + Element elem; elem.type = type; for (auto el : element_ids) { elem.element = el; Matrix elem_coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(elem, elem_coords); elem_size = FEEngine::getElementInradius(elem_coords, type); max_elem_size = std::max(max_elem_size, elem_size); } - AKANTU_DEBUG_INFO("The maximum element size : " - << max_elem_size ); + AKANTU_DEBUG_INFO("The maximum element size : " << max_elem_size); + }*/ + + auto & master_nodes = this->node_selector->getMasterList(); + + for (auto & master : master_nodes) { + Array elements; + this->mesh.getAssociatedElements(master, elements); + + for (auto element : elements) { + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type); + Matrix elem_coords(spatial_dimension, nb_nodes_per_element); + this->coordinatesOfElement(element, elem_coords); + + elem_size = FEEngine::getElementInradius(elem_coords, element.type); + max_elem_size = std::max(max_elem_size, elem_size); + } } + AKANTU_DEBUG_INFO("The maximum element size : " << max_elem_size); + this->max_dd = max_elem_size; this->max_bb = max_elem_size; - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -inline Vector ContactDetector::constructConnectivity(UInt & slave, const Element & master) { - - Vector master_conn = const_cast(this->mesh).getConnectivity(master); +inline Vector +ContactDetector::constructConnectivity(UInt & slave, const Element & master) { + + Vector master_conn = + const_cast(this->mesh).getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { - elem_conn[i] = master_conn[i-1]; - } + elem_conn[i] = master_conn[i - 1]; + } return elem_conn; } /* -------------------------------------------------------------------------- */ -inline void ContactDetector::computeNormalOnElement(const Element & element, Vector & normal) { - +inline 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); } } // to ensure that normal is always outwards from master surface const auto & element_to_subelement = - mesh.getElementToSubelement(element.type)(element.element); + mesh.getElementToSubelement(element.type)(element.element); Vector outside(spatial_dimension); mesh.getBarycenter(element, outside); - + + // check if mesh facets exists for cohesive elements contact Vector inside(spatial_dimension); - mesh.getBarycenter(element_to_subelement[0], inside); + if (mesh.isMeshFacets()) { + mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); + } else { + mesh.getBarycenter(element_to_subelement[0], inside); + } Vector inside_to_outside = outside - inside; auto projection = inside_to_outside.dot(normal); if (projection < 0) { - normal *= -1.0; + normal *= -1.0; } } /* -------------------------------------------------------------------------- */ -inline void ContactDetector::vectorsAlongElement(const Element & el, Matrix & vectors) { +inline 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); } } - } /* -------------------------------------------------------------------------- */ -inline Real ContactDetector::computeGap(Vector & slave, Vector & master, Vector & normal) { +inline Real ContactDetector::computeGap(Vector & slave, + Vector & master, + Vector & normal) { Vector slave_to_master(spatial_dimension); slave_to_master = master - slave; Real gap = slave_to_master.norm(); - + auto projection = slave_to_master.dot(normal); - // if slave node is beneath th master surface + // if slave node is beneath the master surface if (projection > 0) { gap *= -1.0; } return gap; } - - - -/* -------------------------------------------------------------------------- */ -} //akantu +/* -------------------------------------------------------------------------- */ +inline void +ContactDetector::filterBoundaryElements(Array & elements, + Array & boundary_elements) { + + for (auto elem : elements) { + + // to ensure that normal is always outwards from master surface + const auto & element_to_subelement = + mesh.getElementToSubelement(elem.type)(elem.element); + + UInt nb_subelements_regular = 0; + for (auto subelem : element_to_subelement) { + if (subelem == ElementNull) + continue; + + if (subelem.kind() == _ek_regular) { + ++nb_subelements_regular; + } + } + + auto nb_subelements = element_to_subelement.size(); + + if (nb_subelements_regular < nb_subelements) { + boundary_elements.push_back(elem); + } + } +} + +/* -------------------------------------------------------------------------- */ + +} // namespace akantu #endif /* __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__ */ diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 3d8b9a611..ed5ea2e8f 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,652 +1,641 @@ /** * @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 "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( - Mesh & mesh, Array & positions, UInt dim, const ID & id, - const MemoryID & memory_id, std::shared_ptr dof_manager, - const ModelType model_type) - : Model(mesh, model_type, dof_manager, dim, id, memory_id), - current_positions(positions) { + Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, + std::shared_ptr dof_manager, const ModelType model_type) + : Model(mesh, model_type, dof_manager, dim, id, memory_id) { 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, + Model::spatial_dimension - 1, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); - this->detector = std::make_unique( - this->mesh, current_positions, id + ":contact_detector"); + this->detector = + std::make_unique(this->mesh, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -ContactMechanicsModel::ContactMechanicsModel( - Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, - std::shared_ptr dof_manager, const ModelType model_type) - : ContactMechanicsModel(mesh, mesh.getNodes(), dim, id, memory_id, - dof_manager, model_type) {} - /* -------------------------------------------------------------------------- */ 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) { // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->previous_gaps, 1, "previous_gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->normals, spatial_dimension, "normals"); this->allocNodalField(this->tangents, spatial_dimension, "tangents"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_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(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); UInt nb_nodes = mesh.getNbNodes(); this->internal_force->clear(); internal_force->resize(nb_nodes, 0.); // assemble the forces due to contact 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); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { this->detector->search(this->contact_map); for (auto & entry : contact_map) { auto & element = entry.second; if (element.gap < 0) element.gap = std::abs(element.gap); else element.gap = -element.gap; } this->assembleFieldsFromContactMap(); this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search(Array & increment) { this->contact_map.clear(); this->detector->search(this->contact_map); for (auto & entry : contact_map) { auto & element = entry.second; const auto & connectivity = element.connectivity; const auto & normal = element.normal; auto slave_node = connectivity[0]; Vector u_slave(spatial_dimension); auto master_node = connectivity[1]; Vector u_master(spatial_dimension); for (UInt s : arange(spatial_dimension)) { u_slave(s) = increment(slave_node, s); u_master(s) = increment(master_node, s); } // todo check this initial guess auto u = (u_master.norm() > u_slave.norm()) ? u_master : u_slave; Real uv = Math::vectorDot(u.storage(), normal.storage(), spatial_dimension); if (element.gap - uv <= 0) { element.gap = std::abs(element.gap - uv); } else { element.gap = 0.0; } } this->assembleFieldsFromContactMap(); this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleFieldsFromContactMap() { UInt nb_nodes = mesh.getNbNodes(); this->previous_gaps->copy(*(this->gaps)); this->gaps->clear(); gaps->resize(nb_nodes, 0.); normals->resize(nb_nodes, 0.); tangents->resize(nb_nodes, 0.); if (this->contact_map.empty()) return; for (auto & entry : contact_map) { const auto & element = entry.second; auto connectivity = element.connectivity; auto node = connectivity(0); (*gaps)[node] = element.gap; for (UInt i = 0; i < spatial_dimension; ++i) { (*normals)(node, i) = element.normal[i]; (*tangents)(node, i) = element.tangents(0, i); } } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas() { UInt nb_nodes = mesh.getNbNodes(); this->nodal_area->clear(); this->external_force->clear(); nodal_area->resize(nb_nodes, 0.); external_force->resize(nb_nodes, 0.); - ///todo change it for self contact as master is also slave therefore - ///double the area will be considered - this->applyBC( - BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), - this->detector->getSurfaceId()); - + auto & fem_boundary = this->getFEEngineBoundary(); + fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); + fem_boundary.computeNormalsOnIntegrationPoints(_ghost); + this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), - this->detector->getSurfaceId()); + mesh.getElementGroup("contact_surface")); for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); auto & force = std::get<1>(tuple); for (auto & f : force) area += pow(f, 2); area = sqrt(area); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ 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) { 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")); } for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr 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->internal_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["normals"] = this->normals; real_nodal_fields["tangents"] = this->tangents; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["previous_gaps"] = this->previous_gaps; real_nodal_fields["areas"] = this->nodal_area; if (padding_flag) { return this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); } else { return this->mesh.createNodalField(real_nodal_fields[field_name], group_name); } std::shared_ptr field; return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldReal(const std::string & /*field_name*/, const std::string & /*group_name*/, 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; 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/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh index 9eea9f2d5..4d63fc857 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,333 +1,325 @@ /** * @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 "boundary_condition.hh" #include "contact_detector.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.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, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* 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, std::shared_ptr dof_manager = nullptr, const ModelType model_type = ModelType::_contact_mechanics_model); - ContactMechanicsModel( - Mesh & mesh, Array & positions, - UInt spatial_dimension = _all_dimensions, - const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, - std::shared_ptr dof_manager = nullptr, - const ModelType model_type = ModelType::_contact_mechanics_model); - ~ContactMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize all internal arrays for resolutions void initResolutions(); /// initialize the modelType void initModel() override; /// call back for the solver, computes the force residual void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// 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; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Contact Detection */ /* ------------------------------------------------------------------------ */ public: void search(); void search(Array &); void computeNodalAreas(); void assembleFieldsFromContactMap(); /* ------------------------------------------------------------------------ */ /* 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 */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: - std::shared_ptr createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - std::shared_ptr createNodalFieldBool(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; + std::shared_ptr + createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; + + std::shared_ptr + createNodalFieldBool(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); 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); /* ------------------------------------------------------------------------ */ /* 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: /// contact detection class friend class ContactDetector; /// contact resolution class friend class Resolution; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the ContactMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the ContactMechanicsModel::increment vector \warn only consistent /// if ContactMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the ContactMechanics::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the ContactMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the ContactMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the ContactMechanics::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the ContactMechanics::gaps (contact gaps) AKANTU_GET_MACRO(Gaps, *gaps, Array &); /// get the ContactMechanics::normals (normals on slave nodes) AKANTU_GET_MACRO(Normals, *normals, Array &); /// get the ContactMechanics::areas (nodal areas) AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); - /// get the ContactMechanics::internal_force vector (internal forces) - AKANTU_GET_MACRO(CurrentPositions, current_positions, Array &); + /// get contact detector + AKANTU_GET_MACRO_NOT_CONST(ContactDetector, *detector, ContactDetector &); /// get the contat map inline std::map & getContactMap() { return contact_map; } /// inline void setPositions(Array positions) { detector->setPositions(positions); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// tells if the resolutions are instantiated bool are_resolutions_instantiated; /// displacements array Array * displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// contact forces array Array * internal_force{nullptr}; /// external forces array Array * external_force{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; /// array to store gap between slave and master Array * gaps{nullptr}; /// array to store gap from previous iteration Array * previous_gaps{nullptr}; /// array to store normals from master to slave Array * normals{nullptr}; /// array to store tangents on the master element Array * tangents{nullptr}; /// array to store nodal areas Array * nodal_area{nullptr}; - /// array of current position used during update residual - Array & current_positions; - /// 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 "parser.hh" #include "resolution.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/contact_mechanics/node_selector.cc b/src/model/contact_mechanics/node_selector.cc new file mode 100644 index 000000000..49f68aee2 --- /dev/null +++ b/src/model/contact_mechanics/node_selector.cc @@ -0,0 +1,234 @@ +/** + * @file node_selector.cc + * + * @author Mohit Pundir + * + * @date creation: Fri Jun 21 2019 + * @date last modification: Fri Jun 21 2019 + * + * @brief Node selector for contact detector + * + * @section LICENSE + * + * Copyright (©) 2015-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 "node_selector.hh" +#include "model.hh" +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +NodeSelector::NodeSelector(const Model & model) + : Parsable(ParserType::_contact_detector), mesh(model.getMesh()) {} + +/* -------------------------------------------------------------------------- */ +PhysicalSurfaceNodeSelector::PhysicalSurfaceNodeSelector(const Model & model) + : NodeSelector(model) { + + const Parser & parser = getStaticParser(); + + const ParserSection & section = + *(parser.getSubSections(ParserType::_contact_detector).first); + + master = section.getParameterValue("master"); + slave = section.getParameterValue("slave"); + + auto & group = mesh.createElementGroup("contact_surface"); + group.append(mesh.getElementGroup(master)); + group.append(mesh.getElementGroup(slave)); + + group.optimize(); +} + +/* -------------------------------------------------------------------------- */ +Array & PhysicalSurfaceNodeSelector::getMasterList() { + return mesh.getElementGroup(master).getNodeGroup().getNodes(); +} + +/* -------------------------------------------------------------------------- */ +Array & PhysicalSurfaceNodeSelector::getSlaveList() { + return mesh.getElementGroup(slave).getNodeGroup().getNodes(); +} + +#if defined(AKANTU_COHESIVE_ELEMENT) +/* -------------------------------------------------------------------------- */ +CohesiveSurfaceNodeSelector::CohesiveSurfaceNodeSelector(const Model & model) + : NodeSelector(model), mesh_facets(model.getMesh().getMeshFacets()) { + this->mesh.registerEventHandler(*this, _ehp_lowest); +} + +/* -------------------------------------------------------------------------- */ +void CohesiveSurfaceNodeSelector::onNodesAdded(const Array & new_nodes, + const NewNodesEvent & event) { + + if (not aka::is_of_type(event)) + return; + + const auto & cohesive_event = aka::as_type(event); + const auto & old_nodes = cohesive_event.getOldNodesList(); + + UInt nb_new_nodes = new_nodes.size(); + UInt nb_old_nodes = old_nodes.size(); + + //new_nodes_list.reserve(nb_new_nodes + nb_old_nodes); + + for (auto n : arange(nb_new_nodes)) { + new_nodes_list.push_back(new_nodes(n)); + } + + for (auto n : arange(nb_old_nodes)) { + new_nodes_list.push_back(old_nodes(n)); + } + + mesh_facets.fillNodesToElements(mesh.getSpatialDimension() - 1); + + auto & group = + mesh_facets.createElementGroup("contact_surface", _all_dimensions, true); + + for (auto node : new_nodes_list) { + Array all_elements; + mesh_facets.getAssociatedElements(node, all_elements); + + Array mesh_facet_elements; + this->filterBoundaryElements(all_elements, mesh_facet_elements); + + for (auto nb_elem : arange(mesh_facet_elements.size())) + group.add(mesh_facet_elements[nb_elem], true); + } + group.optimize(); +} + +/* -------------------------------------------------------------------------- */ +void CohesiveSurfaceNodeSelector::filterBoundaryElements( + Array & elements, Array & boundary_elements) { + + for (auto elem : elements) { + + // to ensure that normal is always outwards from master surface + const auto & element_to_subelement = + mesh_facets.getElementToSubelement(elem.type)(elem.element); + + UInt nb_subelements_regular = 0; + for (auto subelem : element_to_subelement) { + if (subelem == ElementNull) + continue; + + if (subelem.kind() == _ek_regular) { + ++nb_subelements_regular; + } + } + + auto nb_subelements = element_to_subelement.size(); + + if (nb_subelements_regular < nb_subelements) { + boundary_elements.push_back(elem); + } + } +} + +/* -------------------------------------------------------------------------- */ +Array & CohesiveSurfaceNodeSelector::getMasterList() { + return this->getNewNodesList(); +} + +/* -------------------------------------------------------------------------- */ +Array & CohesiveSurfaceNodeSelector::getSlaveList() { + return this->getNewNodesList(); +} + +/* -------------------------------------------------------------------------- */ +AllSurfaceNodeSelector::AllSurfaceNodeSelector(const Model & model) + : NodeSelector(model), mesh_facets(model.getMesh().getMeshFacets()) { + this->mesh.registerEventHandler(*this, _ehp_lowest); + + const Parser & parser = getStaticParser(); + + const ParserSection & section = + *(parser.getSubSections(ParserType::_contact_detector).first); + + master = section.getParameterValue("master"); + slave = section.getParameterValue("slave"); + + auto & group = mesh_facets.createElementGroup("contact_surface"); + group.append(mesh.getElementGroup(master)); + group.append(mesh.getElementGroup(slave)); + + group.optimize(); +} + +/* -------------------------------------------------------------------------- */ +void AllSurfaceNodeSelector::onNodesAdded(const Array & new_nodes, + const NewNodesEvent & event) { + + if (not aka::is_of_type(event)) + return; + + const auto & cohesive_event = aka::as_type(event); + const auto & old_nodes = cohesive_event.getOldNodesList(); + + UInt nb_new_nodes = new_nodes.size(); + UInt nb_old_nodes = old_nodes.size(); + + new_nodes_list.reserve(nb_new_nodes + nb_old_nodes); + + auto & slave_group = mesh.getElementGroup(slave).getNodeGroup(); + auto & master_group = mesh.getElementGroup(master).getNodeGroup(); + + for (auto n : arange(nb_new_nodes)) { + new_nodes_list.push_back(new_nodes(n)); + slave_group.add(new_nodes(n)); + master_group.add(new_nodes(n)); + } + + for (auto n : arange(nb_old_nodes)) { + new_nodes_list.push_back(old_nodes(n)); + slave_group.add(old_nodes(n)); + master_group.add(old_nodes(n)); + } + + mesh_facets.fillNodesToElements(mesh.getSpatialDimension() - 1); + + auto & group = mesh_facets.getElementGroup("contact_surface"); + + for (auto node : new_nodes_list) { + Array mesh_facet_elements; + mesh_facets.getAssociatedElements(node, mesh_facet_elements); + + for (auto nb_elem : arange(mesh_facet_elements.size())) + group.add(mesh_facet_elements[nb_elem], true); + } + + group.optimize(); +} + +/* -------------------------------------------------------------------------- */ +Array & AllSurfaceNodeSelector::getMasterList() { + return mesh.getElementGroup(master).getNodeGroup().getNodes(); +} + +/* -------------------------------------------------------------------------- */ +Array & AllSurfaceNodeSelector::getSlaveList() { + return mesh.getElementGroup(slave).getNodeGroup().getNodes(); +} + +#endif + +} // namespace akantu diff --git a/src/model/contact_mechanics/node_selector.hh b/src/model/contact_mechanics/node_selector.hh new file mode 100644 index 000000000..b4fb17d06 --- /dev/null +++ b/src/model/contact_mechanics/node_selector.hh @@ -0,0 +1,197 @@ +/** + * @file node_selector.hh + * + * @author Mohit Pundir + * + * @date creation: Fri Jun 21 2019 + * @date last modification: Fri Jun 21 2019 + * + * @brief Node selectors for contact detection + * + * @section LICENSE + * + * Copyright (©) 2015-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 "mesh.hh" +#include "mesh_utils.hh" +#include "parsable.hh" +#if defined(AKANTU_COHESIVE_ELEMENT) +#include "cohesive_element_inserter.hh" +#endif + +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_NODE_SELECTOR_HH__ +#define __AKANTU_NODE_SELECTOR_HH__ + +namespace akantu { +class Model; +class GlobalIdsUpdater; +} // namespace akantu +namespace akantu { + +/** + * main class to assign nodes for contact detection + */ +class NodeSelector : public MeshEventHandler, public Parsable { + /* ------------------------------------------------------------------------ */ + /* Constructor/Destructor */ + /* ------------------------------------------------------------------------ */ +public: + NodeSelector(const Model & model); + + virtual ~NodeSelector() = default; + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + virtual Array & getMasterList() = 0; + + virtual Array & getSlaveList() = 0; + + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ +protected: + Mesh & mesh; +}; + +/* -------------------------------------------------------------------------- */ + +class PhysicalSurfaceNodeSelector : public NodeSelector { + /* ------------------------------------------------------------------------ */ + /* Constructor/Destructor */ + /* ------------------------------------------------------------------------ */ +public: + PhysicalSurfaceNodeSelector(const Model & model); + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + Array & getMasterList() override; + + Array & getSlaveList() override; + + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ +protected: + std::string master; + std::string slave; +}; + +/* -------------------------------------------------------------------------- */ +#if defined(AKANTU_COHESIVE_ELEMENT) +class CohesiveSurfaceNodeSelector : public NodeSelector { + /* ------------------------------------------------------------------------ */ + /* Constructor/Destructor */ + /* ------------------------------------------------------------------------ */ +public: + CohesiveSurfaceNodeSelector(const Model & model); + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +protected: + void onNodesAdded(const Array & nodes_list, + const NewNodesEvent & event) override; + + void filterBoundaryElements(Array & elements, + Array & boundary_elements); + +public: + Array & getMasterList() override; + + Array & getSlaveList() override; + + /* ------------------------------------------------------------------------ + */ + /* */ + /* ------------------------------------------------------------------------ + */ + AKANTU_GET_MACRO_NOT_CONST(NewNodesList, new_nodes_list, Array &); + AKANTU_GET_MACRO(NewNodesList, new_nodes_list, const Array &); + + /* ------------------------------------------------------------------------ + */ + /* Members */ + /* ------------------------------------------------------------------------ + */ +protected: + Mesh & mesh_facets; + + Array new_nodes_list; +}; + +class AllSurfaceNodeSelector : public NodeSelector { + /* ------------------------------------------------------------------------ + */ + /* Constructor/Destructor */ + /* ------------------------------------------------------------------------ + */ +public: + AllSurfaceNodeSelector(const Model & model); + + /* ------------------------------------------------------------------------ + */ + /* Methods */ + /* ------------------------------------------------------------------------ + */ +protected: + void onNodesAdded(const Array & nodes_list, + const NewNodesEvent & event) override; + +public: + Array & getMasterList() override; + + Array & getSlaveList() override; + + /* ------------------------------------------------------------------------ + */ + /* */ + /* ------------------------------------------------------------------------ + */ + AKANTU_GET_MACRO_NOT_CONST(NewNodesList, new_nodes_list, Array &); + AKANTU_GET_MACRO(NewNodesList, new_nodes_list, const Array &); + + /* ------------------------------------------------------------------------ + */ + /* Members */ + /* ------------------------------------------------------------------------ + */ +protected: + std::string master; + + std::string slave; + + Mesh & mesh_facets; + + Array new_nodes_list; +}; + +#endif + +} // namespace akantu + +#endif /* __AKANTU_NODE_SELECTOR_HH__ */ diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index 97603c60d..00196ffd4 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,262 +1,270 @@ /** * @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" #include "sparse_matrix.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("slave", slave, std::string(), _pat_parsable | _pat_readable); registerParam("master", master, std::string(), _pat_parsable | _pat_readable); registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable, "Friciton Coefficient"); - registerParam("two_pass_algorithm", two_pass_algorithm, bool(false), _pat_parsable | _pat_modifiable, - "Two pass algorithm"); - + registerParam("two_pass_algorithm", two_pass_algorithm, bool(false), + _pat_parsable | _pat_modifiable, "Two pass algorithm"); } /* -------------------------------------------------------------------------- */ 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 auto slave_nodes = model.getMesh().getElementGroup(slave).getNodeGroup().getNodes(); + const auto slave_nodes = + model.getContactDetector().getNodeSelector().getSlaveList(); + this->assembleInternalForces(slave_nodes); AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -void Resolution::assembleInternalForces(const Array & local_nodes) { +/* -------------------------------------------------------------------------- */ +void Resolution::assembleInternalForces(const Array & slave_nodes) { AKANTU_DEBUG_IN(); auto & internal_force = const_cast &>(model.getInternalForce()); auto & nodal_area = const_cast &>(model.getNodalArea()); auto & contact_map = model.getContactMap(); - for (auto & slave : local_nodes) { + for (auto & slave : slave_nodes) { if (contact_map.find(slave) == contact_map.end()) continue; - + auto & element = contact_map[slave]; - + const auto & conn = element.connectivity; - + Vector contact_force(conn.size() * spatial_dimension); - + Vector n(conn.size() * spatial_dimension); ResolutionUtils::computeN(n, element); - - computeNormalForce(contact_force, n, element); - if(mu != 0) { - Array t_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); - Array n_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); - Array d_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); + computeNormalForce(contact_force, n, element); + if (mu != 0) { + + Array t_alpha(conn.size() * spatial_dimension, + spatial_dimension - 1); + Array n_alpha(conn.size() * spatial_dimension, + spatial_dimension - 1); + Array d_alpha(conn.size() * spatial_dimension, + spatial_dimension - 1); ResolutionUtils::computeTalpha(t_alpha, element); ResolutionUtils::computeNalpha(n_alpha, element); ResolutionUtils::computeDalpha(d_alpha, n_alpha, t_alpha, element); - + computeFrictionalForce(contact_force, d_alpha, element); } - + ResolutionUtils::assembleToInternalForce(contact_force, internal_force, - nodal_area, element); + nodal_area, element); } AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ void Resolution::assembleStiffnessMatrix(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); auto & stiffness = const_cast(model.getDOFManager().getMatrix("K")); - - const auto local_nodes = - model.getMesh().getElementGroup(name).getNodeGroup().getNodes(); - auto & nodal_area = - const_cast &>(model.getNodalArea()); + // const auto local_nodes = + // model.getMesh().getElementGroup(name).getNodeGroup().getNodes(); + + const auto slave_nodes = + model.getContactDetector().getNodeSelector().getSlaveList(); + + auto & nodal_area = const_cast &>(model.getNodalArea()); auto & contact_map = model.getContactMap(); - for (auto & slave : local_nodes) { + for (auto & slave : slave_nodes) { if (contact_map.find(slave) == contact_map.end()) continue; - + auto & element = contact_map[slave]; - + const auto & conn = element.connectivity; Matrix kc(conn.size() * spatial_dimension, - conn.size() * spatial_dimension); - + conn.size() * spatial_dimension); + Matrix m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1); ResolutionUtils::computeMetricTensor(m_alpha_beta, element.tangents); // normal tangent moduli Vector n(conn.size() * spatial_dimension); Array t_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); Array n_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); Array d_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); - ResolutionUtils::computeN( n, element); - ResolutionUtils::computeTalpha( t_alpha, element); - ResolutionUtils::computeNalpha( n_alpha, element); - ResolutionUtils::computeDalpha( d_alpha, n_alpha, t_alpha, element); + ResolutionUtils::computeN(n, element); + ResolutionUtils::computeTalpha(t_alpha, element); + ResolutionUtils::computeNalpha(n_alpha, element); + ResolutionUtils::computeDalpha(d_alpha, n_alpha, t_alpha, element); computeNormalModuli(kc, n_alpha, d_alpha, n, element); // frictional tangent moduli - if(mu != 0) { + if (mu != 0) { Array t_alpha_beta(conn.size() * spatial_dimension, - (spatial_dimension - 1) * (spatial_dimension -1)); + (spatial_dimension - 1) * + (spatial_dimension - 1)); Array p_alpha(conn.size() * spatial_dimension, - spatial_dimension - 1); + spatial_dimension - 1); Array n_alpha_beta(conn.size() * spatial_dimension, - (spatial_dimension - 1) * (spatial_dimension -1)); + (spatial_dimension - 1) * + (spatial_dimension - 1)); computeFrictionalTraction(m_alpha_beta, element); ResolutionUtils::computeTalphabeta(t_alpha_beta, element); ResolutionUtils::computeNalphabeta(n_alpha_beta, element); ResolutionUtils::computePalpha(p_alpha, element); auto phi = computeNablaOfDisplacement(element); - - computeFrictionalModuli(kc, t_alpha_beta, n_alpha_beta, - n_alpha, d_alpha, phi, n, element); + + computeFrictionalModuli(kc, t_alpha_beta, n_alpha_beta, n_alpha, d_alpha, + phi, n, element); } - + std::vector equations; UInt nb_degree_of_freedom = model.getSpatialDimension(); std::vector areas; for (UInt i : arange(conn.size())) { UInt n = conn[i]; for (UInt j : arange(nb_degree_of_freedom)) { equations.push_back(n * nb_degree_of_freedom + j); areas.push_back(nodal_area[n]); } } for (UInt i : arange(kc.rows())) { UInt row = equations[i]; for (UInt j : arange(kc.cols())) { UInt col = equations[j]; kc(i, j) *= areas[i]; stiffness.add(row, col, kc(i, j)); } } } AKANTU_DEBUG_OUT(); } - + /* -------------------------------------------------------------------------- */ Matrix Resolution::computeNablaOfDisplacement(ContactElement & element) { const auto & type = element.master.type; const auto & conn = element.connectivity; - + auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); Matrix values(spatial_dimension, nb_nodes_per_element); - + auto & displacement = model.getDisplacement(); for (UInt n : arange(nb_nodes_per_element)) { UInt node = conn[n]; for (UInt s : arange(spatial_dimension)) { values(s, n) = displacement(node, s); } } - //Matrix shape_second_derivatives(surface_dimension * surface_dimension, - // nb_nodes_per_element); - + // Matrix shape_second_derivatives(surface_dimension * + // surface_dimension, nb_nodes_per_element); + //#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ //ElementClass::computeDN2DS2(element.projection, shape_second_derivatives) // AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); //#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL - Matrix nabla_u(surface_dimension * surface_dimension, spatial_dimension); - //nabla_u.mul(shape_second_derivatives, values); + Matrix nabla_u(surface_dimension * surface_dimension, + spatial_dimension); + // nabla_u.mul(shape_second_derivatives, values); return nabla_u; } - + } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_cohesive_contact.cc b/src/model/model_couplers/coupler_solid_cohesive_contact.cc index 4c3f37be1..535af65dc 100644 --- a/src/model/model_couplers/coupler_solid_cohesive_contact.cc +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -1,567 +1,564 @@ /** * @file coupler_solid_cohesive_contact.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics and conatct 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 "coupler_solid_cohesive_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidCohesiveContact::CouplerSolidCohesiveContact( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject( "CohesiveFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_cohesive_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_cohesive_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModelCohesive(mesh, Model::spatial_dimension, "solid_mechanics_model_cohesive", 0, this->dof_manager); contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, "contact_mechanics_model", 0, this->dof_manager); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidCohesiveContact::~CouplerSolidCohesiveContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initModel() { getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); // getFEEngine().initShapeFunctions(_not_ghost); // getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidCohesiveContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initSolver(TimeStepSolverType, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); } /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidCohesiveContact::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidCohesiveContact::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidCohesiveContact::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_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 CouplerSolidCohesiveContact::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & pair : contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) external_force(node, s) = contact_force(node, s); } } break; } default: break; } /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & pair : contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) external_force(node, s) = contact_force(node, s); } } break; } default: break; } if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::predictor() { switch (method) { case _explicit_dynamic_contact: { Array displacement(0, Model::spatial_dimension); - Array current_positions(0, Model::spatial_dimension); - auto positions = mesh.getNodes(); - current_positions.copy(positions); + auto & current_positions = contact->getContactDetector().getPositions(); + current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); - // const auto deltas = this->getDOFManager().getSolution("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } - contact->setPositions(current_positions); contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::corrector() { switch (method) { case _implicit_contact: case _explicit_contact: { Array displacement(0, Model::spatial_dimension); Array current_positions(0, Model::spatial_dimension); auto positions = mesh.getNodes(); current_positions.copy(positions); auto us = this->getDOFManager().getDOFs("displacement"); const auto deltas = this->getDOFManager().getSolution("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), deltas, make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & delta = std::get<1>(tuple); const auto & bld = std::get<2>(tuple); auto & cp = std::get<3>(tuple); if (not bld) cp += u + delta; } contact->setPositions(current_positions); contact->search(); break; } default: break; } /*auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); std::stringstream filename; filename << "out" << "-00" << step << ".csv"; std::ofstream outfile(filename.str()); outfile << "x,gap,residual" << std::endl; auto & contact_map = contact->getContactMap(); for (auto & pair: contact_map) { auto & connectivity = pair.second.connectivity; auto node = connectivity(0); if (pair.second.gap > 0) { outfile << positions(node, 0) << "," << pair.second.gap << "," << external_force(node, 1) + internal_force(node, 1) << std::endl; } } outfile.close(); step++;*/ } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidCohesiveContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); switch (method) { case _implicit_contact: { contact->assembleStiffnessMatrix(); break; } default: break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { return solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal( const std::string & field_name, const std::string & group_name, bool padding_flag) { if (field_name == "contact_force" or field_name == "normals" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { return contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { return solid->createNodalFieldReal(field_name, group_name, padding_flag); } std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool( const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldBool(field_name, group_name, padding_flag); std::shared_ptr field; return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc index 5daf92cad..760d95c6e 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,550 +1,527 @@ /** - * @file coupler_solid_contact_explicit.cc + * @file coupler_solid_contact.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics and conatct 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 "coupler_solid_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidContact::CouplerSolidContact( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("CouplerSolidContact", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModel(mesh, Model::spatial_dimension, "solid_mechanics_model", 0, this->dof_manager); contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, "contact_mechanics_model", 0, this->dof_manager); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidContact::~CouplerSolidContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initSolver(TimeStepSolverType, NonLinearSolverType) {} /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidContact::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidContact::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidContact::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_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 CouplerSolidContact::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & pair : contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) external_force(node, s) = contact_force(node, s); } } break; } default: break; } /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { for (auto & pair : contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) external_force(node, s) = contact_force(node, s); } } break; } default: break; } if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::predictor() { switch (method) { case _explicit_dynamic_contact: { Array displacement(0, Model::spatial_dimension); - Array current_positions(0, Model::spatial_dimension); + /*Array current_positions(0, Model::spatial_dimension); auto positions = mesh.getNodes(); - current_positions.copy(positions); + current_positions.copy(positions);*/ + auto & current_positions = contact->getContactDetector().getPositions(); + current_positions.copy(mesh.getNodes()); + auto us = this->getDOFManager().getDOFs("displacement"); - // const auto deltas = this->getDOFManager().getSolution("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } - contact->setPositions(current_positions); + //contact->setPositions(current_positions); contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::corrector() { switch (method) { case _implicit_contact: case _explicit_contact: { Array displacement(0, Model::spatial_dimension); Array current_positions(0, Model::spatial_dimension); auto positions = mesh.getNodes(); current_positions.copy(positions); auto us = this->getDOFManager().getDOFs("displacement"); const auto deltas = this->getDOFManager().getSolution("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), deltas, make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & delta = std::get<1>(tuple); const auto & bld = std::get<2>(tuple); auto & cp = std::get<3>(tuple); if (not bld) cp += u + delta; } contact->setPositions(current_positions); contact->search(); break; } default: break; } - - /*auto & internal_force = solid->getInternalForce(); - auto & external_force = solid->getExternalForce(); - - - - std::stringstream filename; - filename << "out" << "-00" << step << ".csv"; - - std::ofstream outfile(filename.str()); - - outfile << "x,gap,residual" << std::endl; - - auto & contact_map = contact->getContactMap(); - for (auto & pair: contact_map) { - auto & connectivity = pair.second.connectivity; - auto node = connectivity(0); - if (pair.second.gap > 0) { - outfile << positions(node, 0) << "," << pair.second.gap << "," - << external_force(node, 1) + internal_force(node, 1) << std::endl; - } - } - - outfile.close(); - step++;*/ } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); switch (method) { case _implicit_contact: { contact->assembleStiffnessMatrix(); break; } default: break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { return solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { if (field_name == "contact_force" or field_name == "normals" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { return contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { return solid->createNodalFieldReal(field_name, group_name, padding_flag); } std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldBool(field_name, group_name, padding_flag); std::shared_ptr field; return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/test/test_model/test_contact_mechanics_model/material_hertz.dat b/test/test_model/test_contact_mechanics_model/material_hertz.dat new file mode 100644 index 000000000..ca9ddb641 --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/material_hertz.dat @@ -0,0 +1,30 @@ +material elastic [ + name = bot_body + rho = 1.0 # density + E = 1.0 # young's modulu + nu = 0. # poisson's ratio +] + +material elastic [ + name = top_body + rho = 1.0 # density + E = 1.0 # young's modulu + nu = 0. # poisson's ratio +] + +contact_detector [ + type = implicit + master = contact_top + slave = contact_bottom + two_pass_algorithm = false +] + +contact_resolution penalty [ + name = contact_bottom + slave = contact_bottom + master = contact_bottom + mu = 0.0 + epsilon = 100.0 + epsilon_t = 100 + two_pass_algorithm = false +] \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt b/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt index dcf2508b6..f7d40de97 100644 --- a/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt +++ b/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt @@ -1,59 +1,33 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Wed Dec 18 2018 # @date last modification: Wed Apr 29 2019 # # @brief configuration for contact detection 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_mesh (implicit_2d implicit_2d.geo 2 1) -add_mesh (spheres_2d spheres_2d.geo 2 1) -add_mesh (implicit_3d implicit_3d.geo 3 1) -add_mesh (paraview_2d paraview_2d.geo 2 1) +add_mesh (node node.geo 2 1) -register_test(test_implicit_3d - SOURCES test_detection_implicit_3d.cc - DEPENDS implicit_3d - FILES_TO_COPY options.dat - PACKAGE contact_mechanics - ) - -register_test(test_implicit_2d - SOURCES test_detection_implicit_2d.cc - DEPENDS implicit_2d - FILES_TO_COPY options.dat - PACKAGE contact_mechanics - ) - - -register_test(test_paraview - SOURCES test_paraview.cc - DEPENDS paraview_2d - FILES_TO_COPY material.dat - PACKAGE contact_mechanics - ) - - -register_test(test_spheres_2d - SOURCES test_detection_implicit_2d.cc - DEPENDS spheres_2d +register_test(test_node_selector + SOURCES test_node_selector.cc + DEPENDS node FILES_TO_COPY options.dat PACKAGE contact_mechanics ) diff --git a/test/test_model/test_contact_mechanics_model/test_detector/extrinsic_2d.geo b/test/test_model/test_contact_mechanics_model/test_detector/extrinsic_2d.geo deleted file mode 100644 index 4bd8b36d5..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/extrinsic_2d.geo +++ /dev/null @@ -1,29 +0,0 @@ -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_detector/extrinsic_3d.geo b/test/test_model/test_contact_mechanics_model/test_detector/extrinsic_3d.geo deleted file mode 100644 index 40fb74194..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/extrinsic_3d.geo +++ /dev/null @@ -1,55 +0,0 @@ -cl1 = 0.05; -cl2 = 0.1; -cl3 = 0.2; -cl4 = 0.1; -Dy = 0.099; -Dz = 1; -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, 0.6-Dy, 0, cl2}; -Point(5) = {0, 0.6-Dy, -0.5, cl2}; -Point(6) = {0, 0.6-Dy, 0.5, cl2}; - -Point(7) = {0, 0.1, 0, cl4}; -Point(8) = {0.5, 0.1, 0, cl4}; -Point(9) = {-0.5, 0.1, 0, cl4}; -Point(10) = {0, 0.1, -0.5, cl4}; -Point(11) = {0, 0.1, 0.5, cl4}; - -Circle(1) = {3, 4, 1}; -Circle(2) = {1, 4, 2}; -Circle(3) = {6, 4, 1}; -Circle(4) = {1, 4, 5}; -Circle(5) = {3, 4, 6}; -Circle(6) = {6, 4, 2}; -Circle(7) = {2, 4, 5}; -Circle(8) = {5, 4, 3}; -Line Loop(1) = {3, 2, -6}; -Ruled Surface(1) = {1}; -Line Loop(2) = {-2, -7, 4}; -Ruled Surface(2) = {2}; -Line Loop(3) = {-4, -8, -1}; -Ruled Surface(3) = {3}; -Line Loop(4) = {1, -3, -5}; -Ruled Surface(4) = {4}; -Line Loop(5) = {6, 7, 8, 5}; -Plane Surface(5) = {5}; -Surface Loop(1) = {3, 2, 1, 4, 5}; -Volume(1) = {1}; - -Circle(9) = {9, 7, 11}; -Circle(10) = {11, 7, 8}; -Circle(11) = {8, 7, 10}; -Circle(12) = {10, 7, 9}; -Line Loop(6) = {10, 11, 12, 9}; -Plane Surface(6) = {6}; -Extrude {0, -cl3, 0} { - Surface{6}; -} - -Physical Surface("top_surface") = {5}; -Physical Surface("contact_surface") = {6}; -Physical Volume("elastic") = {1}; -Physical Volume("contact_surface") = {2}; -Physical Surface("rigid") = {1, 2, 4, 3}; diff --git a/test/test_model/test_contact_mechanics_model/test_detector/implicit_3d.geo b/test/test_model/test_contact_mechanics_model/test_detector/implicit_3d.geo deleted file mode 100644 index e5ed5b200..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/implicit_3d.geo +++ /dev/null @@ -1,85 +0,0 @@ -cl1 = 0.02; -cl2 = 0.02; -cl3 = 0.05; -cl4 = 0.02; -radius = 0.1; -depth = -0.1; -y = 0.0; -Dz = 1; - -Point(1) = {0, 0, 0, cl1}; -Point(2) = {radius/2, radius/2, 0, cl2}; -Point(3) = {-radius/2, radius/2, 0, cl2}; -Point(4) = {0, radius/2, 0, cl2}; -Point(5) = {0, radius/2, -radius/2, cl2}; -Point(6) = {0, radius/2, radius/2, cl2}; - -Point(12) = {0, y, 0, cl4}; -Point(13) = {radius, y, radius, cl4}; -Point(14) = {-radius, y, radius, cl4}; -Point(15) = {-radius, y, -radius, cl4}; -Point(16) = {radius, y, -radius, cl4}; - -Point(17) = {radius, y + depth, radius, cl2}; -Point(18) = {-radius, y + depth, radius, cl2}; -Point(19) = {-radius, y + depth, -radius, cl2}; -Point(20) = {radius, y + depth, -radius, cl2}; - - -Circle(1) = {3, 4, 1}; -Circle(2) = {1, 4, 2}; -Circle(3) = {6, 4, 1}; -Circle(4) = {1, 4, 5}; -Circle(5) = {3, 4, 6}; -Circle(6) = {6, 4, 2}; -Circle(7) = {2, 4, 5}; -Circle(8) = {5, 4, 3}; - -Line Loop(1) = {3, 2, -6}; -Ruled Surface(1) = {1}; -Line Loop(2) = {-2, -7, 4}; -Ruled Surface(2) = {2}; -Line Loop(3) = {-4, -8, -1}; -Ruled Surface(3) = {3}; -Line Loop(4) = {1, -3, -5}; -Ruled Surface(4) = {4}; -Line Loop(5) = {6, 7, 8, 5}; -Plane Surface(5) = {5}; -Surface Loop(1) = {3, 2, 1, 4, 5}; -Volume(1) = {1}; - -Line(95) = {17, 18}; -Line(96) = {18, 19}; -Line(97) = {19, 20}; -Line(98) = {20, 17}; -Line(100) = {13, 14}; -Line(101) = {13, 17}; -Line(102) = {18, 14}; -Line(103) = {14, 15}; -Line(104) = {19, 15}; -Line(105) = {15, 16}; -Line(106) = {20, 16}; -Line(107) = {16, 13}; - -Line Loop(7) = {95, 96, 97, 98}; -Plane Surface(8) = {7}; -Line Loop(8) = {100, 103, 105, 107}; -Plane Surface(9) = {8}; -Line Loop(9) = {-106, 98, -101, -107}; -Plane Surface(10) = {9}; -Line Loop(10) = {101, 95, 102, -100}; -Plane Surface(11) = {10}; -Line Loop(11) = {-102, 96, 104, -103}; -Plane Surface(12) = {11}; -Line Loop(12) = {-104, 97, 106, -105}; -Plane Surface(13) = {12}; - -Surface Loop(101) = {8, 9, 10, 11, 12, 13}; -Volume(3) = {101}; - -Physical Surface("top_surface") = {5}; -Physical Surface("flat") = {9}; -Physical Surface("bottom_surface") = {8}; -Physical Volume("top_body") = {1}; -Physical Volume("bot_body") = {3}; -Physical Surface("curved") = {1, 2, 4, 3}; diff --git a/test/test_model/test_contact_mechanics_model/test_detector/material.dat b/test/test_model/test_contact_mechanics_model/test_detector/material.dat deleted file mode 100644 index 4114e603f..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/material.dat +++ /dev/null @@ -1,26 +0,0 @@ -material elastic [ - name = bot_body - rho = 7800 # density - E = 1e5 # young's modulu - nu = 0. # poisson's ratio -] - -material elastic [ - name = top_body - rho = 7800 # density - E = 1 # young's modulu - nu = 0. # poisson's ratio -] - -contact_detector [ - type = implicit - master = contact_bottom - slave = contact_top -] - -contact_resolution penalty [ - name = contact_top - mu = 0.0 - epsilon = 1e6 - epsilon_t = 100 -] \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_detector/implicit_2d.geo b/test/test_model/test_contact_mechanics_model/test_detector/node.geo similarity index 73% rename from test/test_model/test_contact_mechanics_model/test_detector/implicit_2d.geo rename to test/test_model/test_contact_mechanics_model/test_detector/node.geo index b908ced78..fa05b0b70 100644 --- a/test/test_model/test_contact_mechanics_model/test_detector/implicit_2d.geo +++ b/test/test_model/test_contact_mechanics_model/test_detector/node.geo @@ -1,33 +1,33 @@ cl1 = 0.0025; -cl2 = 0.005; -cl3 = 0.0025; +cl2 = 0.0025; +cl3 = 0.005; Dy = 0.0; radius = 0.1; -y = 0.0; +y = 0.1; epsilon = 0.0; Point(1) = {0, y, 0, cl1}; Point(2) = {radius, radius + y, 0, cl2}; Point(3) = {-radius, radius + y, 0, cl2}; -Point(4) = {0.25, y -epsilon, 0, cl3}; -Point(5) = {-0.25, y-epsilon, 0, cl3}; +Point(4) = {0.25, 0.1 -epsilon, 0, cl3}; +Point(5) = {-0.25, 0.1-epsilon, 0, cl3}; Point(6) = {-0.25, -0.2, 0, 10*cl2}; Point(7) = {0.25, -0.2, 0, 10*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("flat") = {6}; -Physical Line("curved") = {1, 2}; +Physical Line("contact_bottom") = {6}; +Physical Line("contact_top") = {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/test/test_model/test_contact_mechanics_model/test_detector/options.dat b/test/test_model/test_contact_mechanics_model/test_detector/options.dat index 2e0f21809..ca9ddb641 100644 --- a/test/test_model/test_contact_mechanics_model/test_detector/options.dat +++ b/test/test_model/test_contact_mechanics_model/test_detector/options.dat @@ -1,5 +1,30 @@ +material elastic [ + name = bot_body + rho = 1.0 # density + E = 1.0 # young's modulu + nu = 0. # poisson's ratio +] + +material elastic [ + name = top_body + rho = 1.0 # density + E = 1.0 # young's modulu + nu = 0. # poisson's ratio +] + contact_detector [ type = implicit - master = flat - slave = curved + master = contact_top + slave = contact_bottom + two_pass_algorithm = false +] + +contact_resolution penalty [ + name = contact_bottom + slave = contact_bottom + master = contact_bottom + mu = 0.0 + epsilon = 100.0 + epsilon_t = 100 + two_pass_algorithm = false ] \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_detector/paraview_2d.geo b/test/test_model/test_contact_mechanics_model/test_detector/paraview_2d.geo deleted file mode 100644 index ec7b08464..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/paraview_2d.geo +++ /dev/null @@ -1,43 +0,0 @@ -cl1 = 0.001; -cl2 = 0.005; -cl3 = 0.005; -Dy = 0.0; -radius = 0.1; -y = 0.1; -epsilon = 1e-10; -Point(1) = {0, y + epsilon, 0, cl1}; -Point(2) = {radius, radius + y + epsilon, 0, cl2}; -Point(3) = {-radius, radius + y + epsilon, 0, cl2}; - -Point(11) = {0, y, 0, cl2}; -Point(12) = {radius, -radius + y, 0, cl2}; -Point(13) = {-radius,- radius + y, 0, cl2}; - -Point(8) = {0, radius + y, 0, cl2}; -Point(18) = {0, -radius + y, 0, cl2}; - -Circle(1) = {3, 8, 1}; -Circle(2) = {1, 8, 2}; - -Circle(11) = {13, 18, 11}; -Circle(12) = {11, 18, 12}; - -Line(3) = {2, 8}; -Line(13) = {8, 3}; - -Line(30) = {12, 18}; -Line(31) = {18, 13}; - - -Line Loop(9) = {2, 3, 13, 1}; -Plane Surface(9) = {9}; - -Line Loop(19) = {-12, -30, -31, -11}; -Plane Surface(19) = {19}; - -Physical Line("contact_bottom") = {11, 12}; -Physical Line("contact_top") = {1, 2}; -Physical Line("bottom") = {30, 31}; -Physical Line("top") = {3, 13}; -Physical Surface("bot_body") = {19}; -Physical Surface("top_body") = {9}; \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_detector/test_detection_extrinsic_3d.cc b/test/test_model/test_contact_mechanics_model/test_detector/test_detection_extrinsic_3d.cc deleted file mode 100644 index 615a0c662..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/test_detection_extrinsic_3d.cc +++ /dev/null @@ -1,95 +0,0 @@ -/** - * @file test_contact_detection.cc - * - * @author Mohit Pundir - * - * @date creation: Wed Dec 18 2018 - * @date last modification: Wed Dec 18 2018 - * - * @brief Test for extrinsic detection 3D - * - * @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 "contact_detector.hh" -#include "contact_element.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -Real getAnalyticalGap(Mesh &, UInt &); - -const Real radius = 0.5; -const UInt spatial_dimension = 3; - -using Elements = std::vector; - -int main(int argc, char *argv[]) { - - /*initialize("options.dat", argc, argv); - - Mesh mesh(spatial_dimension); - mesh.read("extrinsic_3d.msh"); - - Elements contact_elements; - - ContactDetector detector(mesh); - detector.search(contact_elements); - - Real epsilon = 1e-5; - for (auto & element: contact_elements) { - auto node = element.getSlave(); - auto gap = element.getGap(); - auto normal = element.getNormal(); - - auto analytical_gap = getAnalyticalGap(mesh, node); - std::cerr << "analytical = " << analytical_gap << " computed = " << gap << std::endl; - std::cerr << " normal = " << normal << std::endl; - //if (abs(analytical_gap - gap) <= epsilon) { - // return EXIT_FAILURE; - //} - }*/ - - - return EXIT_SUCCESS; -} - - -Real getAnalyticalGap(Mesh & mesh, UInt & node) { - Vector pos(spatial_dimension); - Vector center(spatial_dimension); - - center(0) = 0.0; - center(1) = 0.501; - center(2) = 0.0; - - - auto & positions = mesh.getNodes(); - for (UInt s: arange(spatial_dimension)) { - pos(s) = positions(node, s); - } - - Real distance = Math::distance_3d(pos.storage(), center.storage()); - Real analytical_gap = distance - radius; - - return analytical_gap; -} diff --git a/test/test_model/test_contact_mechanics_model/test_detector/test_detection_implicit_2d.cc b/test/test_model/test_contact_mechanics_model/test_detector/test_detection_implicit_2d.cc deleted file mode 100644 index d5534891c..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/test_detection_implicit_2d.cc +++ /dev/null @@ -1,266 +0,0 @@ -/** - * @file test_contact_detection.cc - * - * @author Mohit Pundir - * - * @date creation: Wed Dec 18 2018 - * @date last modification: Wed Dec 18 2018 - * - * @brief Test for extrinsic detection 2D - * - * @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 "contact_detector.hh" -#include "contact_element.hh" -#include "aka_grid_dynamic.hh" -#include -#include - -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - - -const Real radius = 0.1; -const UInt spatial_dimension = 2; - - -auto analyticalCurvedSlave(Mesh & mesh, const UInt & node) { - - auto & positions = mesh.getNodes(); - - Real analytical_gap = positions(node, 1); - - Vector normal(spatial_dimension); - normal[0] = 0.0; - normal[1] = 1.0; - - Vector tangent(spatial_dimension); - tangent[0] = -1.0; - tangent[1] = 0.0; - - return std::make_tuple(analytical_gap, normal, tangent); -} - - -auto analyticalCurvedMaster(Mesh & mesh, const UInt & node) { - - auto & positions = mesh.getNodes(); - - Vector slave_point(spatial_dimension); - slave_point[0] = positions(node, 0); - slave_point[1] = positions(node, 1); - - Real slope = -radius/slave_point[0]; - - Real sign = slave_point[0] < 0 ? -1.0 : 1.0; - - Vector master_point(spatial_dimension); - master_point[0] = sign* radius/std::sqrt(1 + slope * slope); - master_point[1] = slope * master_point[0] + radius; - - auto distance = slave_point - master_point; - Real analytical_gap = Math::norm(spatial_dimension, distance.storage()); - - auto normal = distance.normalize(); - - Vector normal_3d(spatial_dimension + 1); - normal_3d[0] = normal[0]; - normal_3d[1] = normal[1]; - normal_3d[2] = 0.0; - - Vector outward_3d(spatial_dimension + 1); - outward_3d[0] = 0.0; - outward_3d[1] = 0.0; - outward_3d[2] = 1.0; - - auto tangent_3d = outward_3d.crossProduct(normal_3d); - - Vector tangent(spatial_dimension); - tangent[0] = tangent_3d[0]; - tangent[1] = tangent_3d[1]; - - return std::make_tuple(analytical_gap, normal, tangent); -} - - - - -auto checkCurvedSlave(int argc, char *argv[]) { - - initialize("options.dat", argc, argv); - - Mesh mesh(spatial_dimension); - mesh.read("implicit_2d.msh"); - - std::map contact_map; - - ContactDetector detector(mesh); - - detector.setSurfaceId("curved"); - detector.setSurfaceId("flat"); - - SpatialGrid master_grid(spatial_dimension); - SpatialGrid slave_grid(spatial_dimension); - - detector.globalSearch(slave_grid, master_grid); - detector.localSearch(slave_grid, master_grid); - detector.constructContactMap(contact_map); - - for (auto & entry : contact_map) { - const auto & slave = entry.first; - const auto & element = entry.second; - const auto & gap = element.gap; - const auto & normal = element.normal; - const auto & tangent = element.tangents; - - Real analytical_gap; - Vector analytical_normal, analytical_tangent; - - std::tie(analytical_gap, analytical_normal, analytical_tangent) - = analyticalCurvedSlave(mesh, slave); - - Real tolerance = 1e-8; - - auto gap_error = std::abs(gap - analytical_gap); - if (gap_error > tolerance) { - std::cerr << "gap error: " << gap_error << " > " << tolerance - << std::endl; - std::cerr << "gap: " << gap << std::endl - << "analytical gap: " << analytical_gap << std::endl; - return EXIT_FAILURE; - } - - auto normal_error = normal - analytical_normal; - if (std::abs(normal_error[0]) > tolerance or std::abs(normal_error[1]) > tolerance) { - std::cerr << "normal error: " << normal_error << " > " << tolerance - << std::endl; - std::cerr << "normal: " << normal << std::endl - << "analytical normal: " << analytical_normal << std::endl; - return EXIT_FAILURE; - } - - auto tangent_trans = tangent.transpose(); - auto tang = Vector(tangent_trans(0)); - - - auto tangent_error = tang - analytical_tangent; - if (std::abs(tangent_error[0]) > tolerance or std::abs(tangent_error[1]) > tolerance) { - std::cerr << "tangent error: " << tangent_error << " > " << tolerance - << std::endl; - std::cerr << "tangent: " << tang << std::endl - << "analytical tangent: " << analytical_tangent << std::endl; - return EXIT_FAILURE; - } - - } - - return EXIT_SUCCESS; -} - - -auto checkCurvedMaster(int argc, char *argv[]) { - - initialize("options.dat", argc, argv); - - Mesh mesh(spatial_dimension); - mesh.read("implicit_2d.msh"); - - std::map contact_map; - - ContactDetector detector(mesh); - - detector.setSurfaceId("flat"); - detector.setSurfaceId("curved"); - - SpatialGrid master_grid(spatial_dimension); - SpatialGrid slave_grid(spatial_dimension); - - detector.globalSearch(slave_grid, master_grid); - detector.localSearch(slave_grid, master_grid); - detector.constructContactMap(contact_map); - - for (auto & entry : contact_map) { - const auto & slave = entry.first; - const auto & element = entry.second; - const auto & gap = element.gap; - const auto & normal = element.normal; - const auto & tangent = element.tangents; - - Real analytical_gap; - Vector analytical_normal, analytical_tangent; - - std::tie(analytical_gap, analytical_normal, analytical_tangent) - = analyticalCurvedMaster(mesh, slave); - - Real tolerance = 1e-2; - - auto gap_error = std::abs(gap - analytical_gap); - if (gap_error > tolerance) { - std::cerr << "slave node: " << slave << std::endl; - std::cerr << "gap error: " << gap_error << " > " << tolerance - << std::endl; - std::cerr << "gap: " << gap << std::endl - << "analytical gap: " << analytical_gap << std::endl; - return EXIT_FAILURE; - } - - auto normal_error = normal - analytical_normal; - if (std::abs(normal_error[0]) > tolerance or std::abs(normal_error[1]) > tolerance) { - std::cerr << "normal error: " << normal_error << " > " << tolerance - << std::endl; - std::cerr << "normal: " << normal << std::endl - << "analytical normal: " << analytical_normal << std::endl; - return EXIT_FAILURE; - } - - auto tangent_trans = tangent.transpose(); - auto tang = Vector(tangent_trans(0)); - - auto tangent_error = tang - analytical_tangent; - if (std::abs(tangent_error[0]) > tolerance or std::abs(tangent_error[1]) > tolerance) { - std::cerr << "tangent error: " << tangent_error << " > " << tolerance - << std::endl; - std::cerr << "tangent: " << tang << std::endl - << "analytical tangent: " << analytical_tangent << std::endl; - return EXIT_FAILURE; - } - - } - - return EXIT_SUCCESS; -} - - -int main(int argc, char *argv[]) -{ - - //checkCurvedSlave(argc, argv); - checkCurvedMaster(argc, argv); - - return EXIT_SUCCESS; -} - - - - diff --git a/test/test_model/test_contact_mechanics_model/test_detector/test_detection_implicit_3d.cc b/test/test_model/test_contact_mechanics_model/test_detector/test_detection_implicit_3d.cc deleted file mode 100644 index 4cd4fb36b..000000000 --- a/test/test_model/test_contact_mechanics_model/test_detector/test_detection_implicit_3d.cc +++ /dev/null @@ -1,143 +0,0 @@ -/** - * @file test_contact_detection.cc - * - * @author Mohit Pundir - * - * @date creation: Wed Dec 18 2018 - * @date last modification: Wed Dec 18 2018 - * - * @brief Test for extrinsic detection 2D - * - * @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 "contact_detector.hh" -#include "contact_element.hh" -#include "aka_grid_dynamic.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - - -const Real radius = 0.1; -const UInt spatial_dimension = 3; - -auto analyticalCurvedSlave(Mesh & mesh, const UInt & node) { - - auto & positions = mesh.getNodes(); - - Real analytical_gap = positions(node, 1); - - Vector normal(spatial_dimension); - normal[0] = 0.0; - normal[1] = 1.0; - normal[2] = 0.0; - - Vector tangent(spatial_dimension); - tangent[0] = -1.0; - tangent[1] = 0.0; - tangent[2] = 0.0; - - return std::make_tuple(analytical_gap, normal, tangent); -} - - -auto checkCurvedSlave(int argc, char *argv[]) { - - initialize("options.dat", argc, argv); - - Mesh mesh(spatial_dimension); - mesh.read("implicit_3d.msh"); - - std::map contact_map; - - ContactDetector detector(mesh); - - detector.setSurfaceId("curved"); - detector.setSurfaceId("flat"); - - SpatialGrid master_grid(spatial_dimension); - SpatialGrid slave_grid(spatial_dimension); - - detector.globalSearch(slave_grid, master_grid); - detector.localSearch(slave_grid, master_grid); - detector.constructContactMap(contact_map); - - for (auto & entry : contact_map) { - const auto & slave = entry.first; - const auto & element = entry.second; - const auto & gap = element.gap; - const auto & normal = element.normal; - const auto & tangent = element.tangents; - - Real analytical_gap; - Vector analytical_normal, analytical_tangent; - - std::tie(analytical_gap, analytical_normal, analytical_tangent) - = analyticalCurvedSlave(mesh, slave); - - Real tolerance = 1e-8; - - auto gap_error = std::abs(gap - analytical_gap); - if (gap_error > tolerance) { - std::cerr << "gap error: " << gap_error << " > " << tolerance - << std::endl; - std::cerr << "gap: " << gap << std::endl - << "analytical gap: " << analytical_gap << std::endl; - return EXIT_FAILURE; - } - - auto normal_error = normal - analytical_normal; - if (std::abs(normal_error[0]) > tolerance or std::abs(normal_error[1]) > tolerance - or std::abs(normal_error[2]) > tolerance) { - std::cerr << "normal error: " << normal_error << " > " << tolerance - << std::endl; - std::cerr << "normal: " << normal << std::endl - << "analytical normal: " << analytical_normal << std::endl; - return EXIT_FAILURE; - } - - auto tangent_trans = tangent.transpose(); - auto tang_1 = Vector(tangent_trans(0)); - auto tang_2 = Vector(tangent_trans(1)); - - auto tangent_error_1 = tang_1 - analytical_tangent; - auto tangent_error_2 = tang_2 - analytical_tangent; - - if (std::abs(tangent_error_1[0]) > tolerance or std::abs(tangent_error_1[1]) > tolerance - or std::abs(tangent_error_1[2]) > tolerance) { - std::cerr << "tangent error: " << tangent_error_1 << " > " << tolerance - << std::endl; - std::cerr << "tangent: " << tang_1 << std::endl - << "analytical tangent: " << analytical_tangent << std::endl; - return EXIT_FAILURE; - } - - } -} - - -int main(int argc, char *argv[]) -{ - checkCurvedSlave(argc, argv); - return 0; -} diff --git a/test/test_model/test_contact_mechanics_model/test_detector/test_detection_intrinsic_3d.cc b/test/test_model/test_contact_mechanics_model/test_detector/test_detection_intrinsic_3d.cc deleted file mode 100644 index e69de29bb..000000000 diff --git a/test/test_model/test_contact_mechanics_model/test_detector/test_node_selector.cc b/test/test_model/test_contact_mechanics_model/test_detector/test_node_selector.cc new file mode 100644 index 000000000..c0c176614 --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_detector/test_node_selector.cc @@ -0,0 +1,27 @@ +#include "aka_common.hh" +#include "node_selector.hh" +#include "contact_mechanics_model.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char * argv[]) { + + const UInt spatial_dimension = 3; + + initialize("options.dat", argc, argv); + + Mesh mesh(spatial_dimension); + mesh.read("node.msh"); + + ContactMechanicsModel model(mesh); + + PhysicalSurfaceNodeSelector selector(model); + auto & slave = selector.getSlaveList(); + + for (auto & s : slave) { + std::cerr << s << std::endl; + } + + return 0; +} diff --git a/test/test_model/test_contact_mechanics_model/test_detector/test_paraview.cc b/test/test_model/test_contact_mechanics_model/test_hertz_2d.cc similarity index 57% rename from test/test_model/test_contact_mechanics_model/test_detector/test_paraview.cc rename to test/test_model/test_contact_mechanics_model/test_hertz_2d.cc index abfa4955e..de07574bc 100644 --- a/test/test_model/test_contact_mechanics_model/test_detector/test_paraview.cc +++ b/test/test_model/test_contact_mechanics_model/test_hertz_2d.cc @@ -1,111 +1,124 @@ /** - * @file test_contact_mechanics_model.cc + * @file test_hertz_2d.cc * * @author Mohit Pundir * * @date creation: Tue Apr 30 2019 * @date last modification: Tue Apr 30 2019 * - * @brief Test for contact mechanics model class + * @brief Hertz Test for contact mechanics model 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 "solid_mechanics_model.hh" #include "contact_mechanics_model.hh" #include "coupler_solid_contact.hh" #include "non_linear_solver.hh" +#include "dumper_text.hh" +#include "dumper_variable.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { + + UInt max_steps = 5000; + UInt imposing_steps = 1000; + UInt damping_interval = 10; + Real damping_ratio = 0.9; + Real pressure = 0.01; - const UInt spatial_dimension = 2; - initialize("material.dat", argc, argv); + std::string mesh_file = "hertz_2d.msh"; + std::string material_file = "material_hertz.dat"; - auto increment = 1e-3; - auto nb_steps = 10; + const UInt spatial_dimension = 2; + initialize(material_file, argc, argv); + Mesh mesh(spatial_dimension); - mesh.read("paraview_2d.msh"); - + mesh.read(mesh_file); + CouplerSolidContact coupler(mesh); auto & solid = coupler.getSolidMechanicsModel(); auto & contact = coupler.getContactMechanicsModel(); - solid.initFull( _analysis_method = _static); - contact.initFull(_analysis_method = _implicit_contact); + auto && selector = std::make_shared>( + "physical_names",solid); - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bot_body"); - solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "bot_body"); + solid.initFull( _analysis_method = _explicit_lumped_mass); + contact.initFull(_analysis_method = _explicit_dynamic_contact); - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top"); - - coupler.initFull(_analysis_method = _implicit_contact); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top_body"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top_body"); + + coupler.initFull(_analysis_method = _explicit_dynamic_contact); + + Vector weight = {0, pressure}; + solid.applyBC(BC::Neumann::FromSameDim(weight), "bottom"); - coupler.setBaseName("test-circles-2d"); + Real time_step = solid.getStableTimeStep(); + time_step *= 0.8; + + coupler.setTimeStep(time_step); + + coupler.setBaseName("test-hertz-2d"); coupler.addDumpFieldVector("displacement"); coupler.addDumpFieldVector("normals"); coupler.addDumpFieldVector("contact_force"); coupler.addDumpFieldVector("external_force"); coupler.addDumpFieldVector("internal_force"); coupler.addDumpField("gaps"); + coupler.addDumpField("areas"); coupler.addDumpField("blocked_dofs"); coupler.addDumpField("grad_u"); coupler.addDumpField("stress"); - auto & before_assembly = - const_cast(coupler.getDOFManager().getNewMatrix("K", _symmetric)); - - solid.assembleStiffnessMatrix(); + auto & velocity = solid.getVelocity(); - auto & solid_assembly = - const_cast(coupler.getDOFManager().getMatrix("K")); - solid_assembly.saveMatrix("solid_assembly.mtx"); - - auto & displacement = solid.getDisplacement(); - - contact.search(displacement); - contact.assembleStiffnessMatrix(); + auto & gaps = contact.getGaps(); + auto & areas = contact.getNodalArea(); - auto contact_map = contact.getContactMap(); - auto nb_contacts = contact_map.size(); - - auto & contact_assembly = - const_cast(coupler.getDOFManager().getMatrix("K")); - contact_assembly.saveMatrix("contact_assembly.mtx"); - - solid.assembleInternalForces(); - contact.assembleInternalForces(); - - coupler.dump(); + for (UInt s : arange(max_steps)) { + std::cerr << "Step " << s << std::endl; + + coupler.solveStep(); + + if (s % damping_interval == 0) { + for (auto & v : make_view(velocity)) + v *= damping_ratio; + } + if (s % 100 == 0) { + coupler.dump(); + } + + } finalize(); return EXIT_SUCCESS; + }