diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 2161dc89e..4bb30c67d 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,504 +1,505 @@ /** * @file aka_common.hh * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Thu Jan 21 2016 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ using MemoryID = UInt; // using Surface = std::string; // using SurfacePair= std::pair; // using SurfacePairList = std::list; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; #define AKANTU_PP_ENUM(s, data, i, elem) \ BOOST_PP_TUPLE_REM() \ elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data))) } // namespace akantu #if (defined(__GNUC__) || defined(__GNUG__)) #define AKA_GCC_VERSION \ (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) #if AKA_GCC_VERSION < 60000 #define AKANTU_ENUM_HASH(type_name) \ namespace std { \ template <> struct hash<::akantu::type_name> { \ using argument_type = ::akantu::type_name; \ size_t operator()(const argument_type & e) const noexcept { \ auto ue = underlying_type_t(e); \ return uh(ue); \ } \ \ private: \ const hash> uh{}; \ }; \ } #else #define AKANTU_ENUM_HASH(type_name) #endif // AKA_GCC_VERSION #endif // GNU #include "aka_element_classes_info.hh" namespace akantu { #define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) #define AKANTU_PP_TYPE_TO_STR(s, data, elem) \ ({BOOST_PP_CAT(data::_, elem), BOOST_PP_STRINGIZE(elem)}) #define AKANTU_PP_STR_TO_TYPE(s, data, elem) \ ({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data::_, elem)}) #define AKANTU_ENUM_DECLARE(type_name, list) \ enum class type_name { \ BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \ }; #define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ } \ AKANTU_ENUM_HASH(type_name) \ namespace aka { \ inline std::string to_string(const ::akantu::type_name & type) { \ static std::unordered_map<::akantu::type_name, std::string> convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, \ ::akantu::type_name, list))}; \ return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << aka::to_string(type); \ return stream; \ } #define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ inline std::istream & operator>>(std::istream & stream, type_name & type) { \ std::string str; \ stream >> str; \ static std::unordered_map convert{ \ BOOST_PP_SEQ_FOR_EACH_I( \ AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, type_name, list))}; \ type = convert.at(str); \ return stream; \ } /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpacialDirection { _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) + (structural_mechanics_model) \ + (embedded_model) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) /// 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 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent _nls_lumped, ///< Case of lumped mass or equivalent matrix _nls_auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; /// Define the node/dof type enum NodeType : Int { _nt_pure_gost = -3, _nt_master = -2, _nt_normal = -1 }; /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_size, //--- SolidMechanicsModel tags --- _gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _gst_smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _gst_smm_uv, ///< synchronization of the nodal velocities and displacement _gst_smm_res, ///< synchronization of the nodal residual _gst_smm_init_mat, ///< synchronization of the data to initialize materials _gst_smm_stress, ///< synchronization of the stresses to compute the internal /// forces _gst_smmc_facets, ///< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, ///< synchronization of facet global connectivity _gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet /// synch _gst_smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _gst_giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _gst_ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gst_gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _gst_htm_capacity, ///< synchronization of the nodal heat capacity _gst_htm_temperature, ///< synchronization of the nodal temperature _gst_htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _gst_htm_phi, ///< synchronization of the nodal level set value phi _gst_htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _gst_mnl_for_average, ///< synchronization of data to average in non local /// material _gst_mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _gst_nh_criterion, // --- General tags --- _gst_test, ///< Test tag _gst_user_1, ///< tag for user simulations _gst_user_2, ///< tag for user simulations _gst_material_id, ///< synchronization of the material ids _gst_for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _gst_cf_nodal, ///< synchronization of disp, velo, and current position _gst_cf_incr, ///< synchronization of increment // --- Solver tags --- _gst_solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; } 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 /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ struct TensorTrait {}; /* -------------------------------------------------------------------------- */ template using is_tensor = std::is_base_of; /* -------------------------------------------------------------------------- */ template using is_scalar = std::is_arithmetic; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ // #if defined(__INTEL_COMPILER) // /// remark #981: operands are evaluated in unspecified order // #pragma warning(disable : 981) // /// remark #383: value copied to temporary, reference to temporary used // #pragma warning(disable : 383) // #endif // defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ } // namespace akantu #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic number * is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/geometry/aabb_primitives/aabb_primitive.cc b/src/geometry/aabb_primitives/aabb_primitive.cc index 3470de999..56b56a7bd 100644 --- a/src/geometry/aabb_primitives/aabb_primitive.cc +++ b/src/geometry/aabb_primitives/aabb_primitive.cc @@ -1,49 +1,49 @@ /** * @file aabb_primitive.cc * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Jan 04 2013 * @date last modification: Thu Jan 14 2016 * * @brief Macro classe (primitive) for AABB CGAL algos * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #include "aabb_primitive.hh" namespace akantu { Triangle_primitive::Point Triangle_primitive::reference_point() const { return primitive.vertex(0); } Line_arc_primitive::Point Line_arc_primitive::reference_point() const { Real x = to_double(primitive.source().x()); Real y = to_double(primitive.source().y()); Real z = to_double(primitive.source().z()); - return Spherical::Point_3(x, y, z); + return cgal::Spherical::Point_3(x, y, z); } } // akantu diff --git a/src/geometry/aabb_primitives/aabb_primitive.hh b/src/geometry/aabb_primitives/aabb_primitive.hh index 4a9a9f506..2fea28c96 100644 --- a/src/geometry/aabb_primitives/aabb_primitive.hh +++ b/src/geometry/aabb_primitives/aabb_primitive.hh @@ -1,90 +1,90 @@ /** * @file aabb_primitive.hh * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Mar 13 2015 * @date last modification: Thu Jan 14 2016 * * @brief Macro classe (primitive) for AABB CGAL algos * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AABB_PRIMITIVE_HH__ #define __AKANTU_AABB_PRIMITIVE_HH__ #include "aka_common.hh" #include "line_arc.hh" #include "triangle.hh" #include "tetrahedron.hh" #include "mesh_geom_common.hh" namespace akantu { /** * This macro defines a class that is used in the CGAL AABB tree algorithm. * All the `typedef`s and methods are required by the AABB module. * * The member variables are * - the id of the element associated to the primitive * - the geometric primitive of the element * * @param name the name of the primitive type * @param kernel the name of the kernel used */ #define AKANTU_AABB_CLASS(name, kernel) \ class name##_primitive { \ typedef std::list< name >::iterator Iterator; \ \ public: \ typedef UInt Id; \ typedef kernel::Point_3 Point; \ typedef kernel::name##_3 Datum; \ \ public: \ name##_primitive() : meshId(0), primitive() {} \ name##_primitive(Iterator it) : meshId(it->id()), primitive(*it) {} \ \ public: \ const Datum & datum() const { return primitive; } \ Point reference_point() const; \ const Id & id() const { return meshId; } \ \ protected: \ Id meshId; \ name primitive; \ \ } // If the primitive is supported by CGAL::intersection() then the // implementation process is really easy with this macro -AKANTU_AABB_CLASS(Triangle, Cartesian); -AKANTU_AABB_CLASS(Line_arc, Spherical); +AKANTU_AABB_CLASS(Triangle, cgal::Cartesian); +AKANTU_AABB_CLASS(Line_arc, cgal::Spherical); #undef AKANTU_AABB_CLASS } // akantu #endif // __AKANTU_AABB_PRIMITIVE_HH__ diff --git a/src/geometry/geom_helper_functions.hh b/src/geometry/geom_helper_functions.hh index 14e933bf3..0d68d20e2 100644 --- a/src/geometry/geom_helper_functions.hh +++ b/src/geometry/geom_helper_functions.hh @@ -1,104 +1,112 @@ /** * @file geom_helper_functions.hh * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Jan 04 2013 * @date last modification: Thu Jan 14 2016 * * @brief Helper functions for the computational geometry algorithms * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef _AKANTU_GEOM_HELPER_FUNCTIONS_HH__ #define _AKANTU_GEOM_HELPER_FUNCTIONS_HH__ #include "aka_common.hh" #include "aka_math.hh" #include "tree_type_helper.hh" #include "mesh_geom_common.hh" namespace akantu { /// Fuzzy compare of two points template inline bool comparePoints(const Point & a, const Point & b) { - return Math::are_float_equal(a.x(), b.x()) && - Math::are_float_equal(a.y(), b.y()) && - Math::are_float_equal(a.z(), b.z()); + return Math::are_float_equal(a.x(), b.x()) && + Math::are_float_equal(a.y(), b.y()) && + Math::are_float_equal(a.z(), b.z()); } template <> -inline bool comparePoints(const Spherical::Circular_arc_point_3 & a, - const Spherical::Circular_arc_point_3 & b) { - return Math::are_float_equal(CGAL::to_double(a.x()), CGAL::to_double(b.x())) && - Math::are_float_equal(CGAL::to_double(a.y()), CGAL::to_double(b.y())) && - Math::are_float_equal(CGAL::to_double(a.z()), CGAL::to_double(b.z())); +inline bool comparePoints(const cgal::Spherical::Circular_arc_point_3 & a, + const cgal::Spherical::Circular_arc_point_3 & b) { + return Math::are_float_equal(CGAL::to_double(a.x()), + CGAL::to_double(b.x())) && + Math::are_float_equal(CGAL::to_double(a.y()), + CGAL::to_double(b.y())) && + Math::are_float_equal(CGAL::to_double(a.z()), CGAL::to_double(b.z())); } /// Fuzzy compare of two segments template -inline bool compareSegments(const CGAL::Segment_3 & a, const CGAL::Segment_3 & b) { - return (comparePoints(a.source(), b.source()) && comparePoints(a.target(), b.target())) || - (comparePoints(a.source(), b.target()) && comparePoints(a.target(), b.source())); +inline bool compareSegments(const CGAL::Segment_3 & a, + const CGAL::Segment_3 & b) { + return (comparePoints(a.source(), b.source()) && + comparePoints(a.target(), b.target())) || + (comparePoints(a.source(), b.target()) && + comparePoints(a.target(), b.source())); } -typedef Cartesian K; - /// Compare segment pairs -inline bool compareSegmentPairs(const std::pair & a, const std::pair & b) { +inline bool +compareSegmentPairs(const std::pair & a, + const std::pair & b) { return compareSegments(a.first, b.first); } /// Pair ordering operator based on first member struct segmentPairsLess { - inline bool operator()(const std::pair & a, const std::pair & b) { - return CGAL::compare_lexicographically(a.first.min(), b.first.min()) || CGAL::compare_lexicographically(a.first.max(), b.first.max()); + inline bool + operator()(const std::pair & a, + const std::pair & b) { + return CGAL::compare_lexicographically(a.first.min(), b.first.min()) || + CGAL::compare_lexicographically(a.first.max(), b.first.max()); } }; /* -------------------------------------------------------------------------- */ /* Predicates */ /* -------------------------------------------------------------------------- */ /// Predicate used to determine if two segments are equal class IsSameSegment { public: - IsSameSegment(const K::Segment_3 & segment): - segment(segment) - {} + IsSameSegment(const cgal::Cartesian::Segment_3 & segment) + : segment(segment) {} - bool operator()(const std::pair & test_pair) { + bool + operator()(const std::pair & test_pair) { return compareSegments(segment, test_pair.first); } + protected: - const K::Segment_3 segment; + const cgal::Cartesian::Segment_3 segment; }; } // akantu #endif // _AKANTU_GEOM_HELPER_FUNCTIONS_HH__ - diff --git a/src/geometry/mesh_geom_common.hh b/src/geometry/mesh_geom_common.hh index 72289b9c9..c514a67e4 100644 --- a/src/geometry/mesh_geom_common.hh +++ b/src/geometry/mesh_geom_common.hh @@ -1,58 +1,57 @@ /** * @file mesh_geom_common.hh * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Jan 04 2013 * @date last modification: Thu Jan 14 2016 * * @brief Common file for MeshGeom module * * @section LICENSE * * Copyright (©) 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 . * */ #ifndef __AKANTU_MESH_GEOM_COMMON_HH__ #define __AKANTU_MESH_GEOM_COMMON_HH__ #include "aka_common.hh" #include #include #include #include #include #include -#include namespace akantu { -typedef CGAL::Simple_cartesian Cartesian; +namespace cgal { + using Cartesian = CGAL::Simple_cartesian; -typedef CGAL::Quotient NT; -typedef CGAL::Spherical_kernel_3, CGAL::Algebraic_kernel_for_spheres_2_3 > Spherical; - -//typedef CGAL::Lazy_exact_nt > NT; -//typedef CGAL::Spherical_kernel_3, CGAL::Algebraic_kernel_for_spheres_2_3 > Spherical; + using Spherical = CGAL::Spherical_kernel_3< + CGAL::Simple_cartesian>, + CGAL::Algebraic_kernel_for_spheres_2_3>>; +} // cgal } // akantu #endif // __AKANTU_MESH_GEOM_COMMON_HH__ diff --git a/src/geometry/mesh_geom_factory_tmpl.hh b/src/geometry/mesh_geom_factory_tmpl.hh index dfd7f01e3..f7ece3cfa 100644 --- a/src/geometry/mesh_geom_factory_tmpl.hh +++ b/src/geometry/mesh_geom_factory_tmpl.hh @@ -1,259 +1,259 @@ /** * @file mesh_geom_factory_tmpl.hh * * @author Lucas Frerot * @author Clement Roux * @author Marco Vocialta * * @date creation: Fri Feb 27 2015 * @date last modification: Thu Jan 14 2016 * * @brief Class for constructing the CGAL primitives of a mesh * * @section LICENSE * * Copyright (©) 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 . * */ #ifndef __AKANTU_MESH_GEOM_FACTORY_TMPL_HH__ #define __AKANTU_MESH_GEOM_FACTORY_TMPL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh_geom_common.hh" #include "mesh_geom_factory.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template MeshGeomFactory::MeshGeomFactory(Mesh & mesh) : MeshGeomAbstract(mesh), data_tree(NULL), primitive_list() {} template MeshGeomFactory::~MeshGeomFactory() { delete data_tree; } /** * This function loops over the elements of `type` in the mesh and creates the * AABB tree of geometrical primitves (`data_tree`). */ template void MeshGeomFactory::constructData(GhostType ghost_type) { AKANTU_DEBUG_IN(); primitive_list.clear(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, ghost_type); const Array & nodes = mesh.getNodes(); UInt el_index = 0; Array::const_vector_iterator it = connectivity.begin(nb_nodes_per_element); Array::const_vector_iterator end = connectivity.end(nb_nodes_per_element); Matrix node_coordinates(dim, nb_nodes_per_element); // This loop builds the list of primitives for (; it != end ; ++it, ++el_index) { const Vector & el_connectivity = *it; for (UInt i = 0 ; i < nb_nodes_per_element ; i++) for (UInt j = 0 ; j < dim ; j++) node_coordinates(j, i) = nodes(el_connectivity(i), j); // the unique elemental id assigned to the primitive is the // linearized element index over ghost type addPrimitive(node_coordinates, el_index); } delete data_tree; // This condition allows the use of the mesh geom module // even if types are not compatible with AABB tree algorithm if (TreeTypeHelper::is_valid) data_tree = new typename TreeTypeHelper::tree(primitive_list.begin(), primitive_list.end()); AKANTU_DEBUG_OUT(); } template void MeshGeomFactory::addPrimitive(const Matrix & node_coordinates, UInt id) { this->addPrimitive(node_coordinates, id, this->primitive_list); } -// (2D, _triangle_3) decomposed into Triangle +// (2D, _triangle_3) decomposed into Triangle template<> -inline void MeshGeomFactory<2, _triangle_3, Triangle, Cartesian>::addPrimitive( +inline void MeshGeomFactory<2, _triangle_3, Triangle, cgal::Cartesian>::addPrimitive( const Matrix & node_coordinates, UInt id, - TreeTypeHelper, Cartesian>::container_type & list) { + TreeTypeHelper, cgal::Cartesian>::container_type & list) { - TreeTypeHelper, Cartesian>::point_type + TreeTypeHelper, cgal::Cartesian>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), 0.), b(node_coordinates(0, 1), node_coordinates(1, 1), 0.), c(node_coordinates(0, 2), node_coordinates(1, 2), 0.); - Triangle t(a, b, c); + Triangle t(a, b, c); t.setId(id); list.push_back(t); } -// (2D, _triangle_6) decomposed into Triangle +// (2D, _triangle_6) decomposed into Triangle template<> -inline void MeshGeomFactory<2, _triangle_6, Triangle, Cartesian>::addPrimitive( +inline void MeshGeomFactory<2, _triangle_6, Triangle, cgal::Cartesian>::addPrimitive( const Matrix & node_coordinates, UInt id, - TreeTypeHelper, Cartesian>::container_type & list) { + TreeTypeHelper, cgal::Cartesian>::container_type & list) { - TreeTypeHelper, Cartesian>::point_type + TreeTypeHelper, cgal::Cartesian>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), 0.), b(node_coordinates(0, 1), node_coordinates(1, 1), 0.), c(node_coordinates(0, 2), node_coordinates(1, 2), 0.); - Triangle t(a, b, c); + Triangle t(a, b, c); t.setId(id); list.push_back(t); } -// (2D, _triangle_3) decomposed into Line_arc +// (2D, _triangle_3) decomposed into Line_arc template<> -inline void MeshGeomFactory<2, _triangle_3, Line_arc, Spherical>::addPrimitive( +inline void MeshGeomFactory<2, _triangle_3, Line_arc, cgal::Spherical>::addPrimitive( const Matrix & node_coordinates, UInt id, - TreeTypeHelper, Spherical>::container_type & list) { + TreeTypeHelper, cgal::Spherical>::container_type & list) { - TreeTypeHelper, Spherical>::point_type + TreeTypeHelper, cgal::Spherical>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), 0.), b(node_coordinates(0, 1), node_coordinates(1, 1), 0.), c(node_coordinates(0, 2), node_coordinates(1, 2), 0.); /*std::cout << "elem " << id << " node 1 : x_node=" << node_coordinates(0, 0) << ", x_arc_node=" << a.x() << ", y_node=" << node_coordinates(1, 0) << ", y_arc_node=" << a.y() << std::endl ; std::cout << "elem " << id << " node 2 : x_node=" << node_coordinates(0, 1) << ", x_arc_node=" << b.x() << ", y_node=" << node_coordinates(1, 1) << ", y_arc_node=" << b.y() << std::endl ; std::cout << "elem " << id << " node 2 : x_node=" << node_coordinates(0, 2) << ", x_arc_node=" << c.x() << ", y_node=" << node_coordinates(1, 2) << ", y_arc_node=" << c.y() << std::endl ;*/ - CGAL::Line_3 l1(a, b), l2(b, c), l3(c, a); - Line_arc s1(l1,a, b), s2(l2, b, c), s3(l3, c, a); + CGAL::Line_3 l1(a, b), l2(b, c), l3(c, a); + Line_arc s1(l1,a, b), s2(l2, b, c), s3(l3, c, a); s1.setId(id); s1.setSegId(0); s2.setId(id); s2.setSegId(1); s3.setId(id); s3.setSegId(2); list.push_back(s1); list.push_back(s2); list.push_back(s3); } #if defined(AKANTU_IGFEM) -// (2D, _igfem_triangle_4) decomposed into Line_arc +// (2D, _igfem_triangle_4) decomposed into Line_arc template<> -inline void MeshGeomFactory<2, _igfem_triangle_4, Line_arc, Spherical>::addPrimitive( +inline void MeshGeomFactory<2, _igfem_triangle_4, Line_arc, cgal::Spherical>::addPrimitive( const Matrix & node_coordinates, UInt id, - TreeTypeHelper, Spherical>::container_type & list) { + TreeTypeHelper, cgal::Spherical>::container_type & list) { - TreeTypeHelper, Spherical>::point_type + TreeTypeHelper, cgal::Spherical>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), 0.), b(node_coordinates(0, 1), node_coordinates(1, 1), 0.), c(node_coordinates(0, 2), node_coordinates(1, 2), 0.); - CGAL::Line_3 l1(a, b), l2(b, c), l3(c, a); - Line_arc s1(l1,a, b), s2(l2, b, c), s3(l3, c, a); + CGAL::Line_3 l1(a, b), l2(b, c), l3(c, a); + Line_arc s1(l1,a, b), s2(l2, b, c), s3(l3, c, a); s1.setId(id); s1.setSegId(0); s2.setId(id); s2.setSegId(1); s3.setId(id); s3.setSegId(2); list.push_back(s1); list.push_back(s2); list.push_back(s3); } -// (2D, _igfem_triangle_4) decomposed into Line_arc +// (2D, _igfem_triangle_4) decomposed into Line_arc template<> -inline void MeshGeomFactory<2, _igfem_triangle_5, Line_arc, Spherical>::addPrimitive( +inline void MeshGeomFactory<2, _igfem_triangle_5, Line_arc, cgal::Spherical>::addPrimitive( const Matrix & node_coordinates, UInt id, - TreeTypeHelper, Spherical>::container_type & list) { + TreeTypeHelper, cgal::Spherical>::container_type & list) { - TreeTypeHelper, Spherical>::point_type + TreeTypeHelper, cgal::Spherical>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), 0.), b(node_coordinates(0, 1), node_coordinates(1, 1), 0.), c(node_coordinates(0, 2), node_coordinates(1, 2), 0.); - CGAL::Line_3 l1(a, b), l2(b, c), l3(c, a); - Line_arc s1(l1,a, b), s2(l2, b, c), s3(l3, c, a); + CGAL::Line_3 l1(a, b), l2(b, c), l3(c, a); + Line_arc s1(l1,a, b), s2(l2, b, c), s3(l3, c, a); s1.setId(id); s1.setSegId(0); s2.setId(id); s2.setSegId(1); s3.setId(id); s3.setSegId(2); list.push_back(s1); list.push_back(s2); list.push_back(s3); } #endif -// (3D, _tetrahedron_4) decomposed into Triangle +// (3D, _tetrahedron_4) decomposed into Triangle template<> -inline void MeshGeomFactory<3, _tetrahedron_4, Triangle, Cartesian>::addPrimitive( +inline void MeshGeomFactory<3, _tetrahedron_4, Triangle, cgal::Cartesian>::addPrimitive( const Matrix & node_coordinates, UInt id, - TreeTypeHelper, Cartesian>::container_type & list) { + TreeTypeHelper, cgal::Cartesian>::container_type & list) { - TreeTypeHelper, Cartesian>::point_type + TreeTypeHelper, cgal::Cartesian>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), node_coordinates(2, 0)), b(node_coordinates(0, 1), node_coordinates(1, 1), node_coordinates(2, 1)), c(node_coordinates(0, 2), node_coordinates(1, 2), node_coordinates(2, 2)), d(node_coordinates(0, 3), node_coordinates(1, 3), node_coordinates(2, 3)); - Triangle + Triangle t1(a, b, c), t2(b, c, d), t3(c, d, a), t4(d, a, b); t1.setId(id); t2.setId(id); t3.setId(id); t4.setId(id); list.push_back(t1); list.push_back(t2); list.push_back(t3); list.push_back(t4); } } // akantu #endif // __AKANTU_MESH_GEOM_FACTORY_TMPL_HH__ diff --git a/src/geometry/mesh_segment_intersector.hh b/src/geometry/mesh_segment_intersector.hh index 8ef17445d..6d913f561 100644 --- a/src/geometry/mesh_segment_intersector.hh +++ b/src/geometry/mesh_segment_intersector.hh @@ -1,103 +1,103 @@ /** * @file mesh_segment_intersector.hh * * @author Lucas Frerot * @author Clement Roux * @author Marco Vocialta * * @date creation: Wed Apr 29 2015 * @date last modification: Thu Jan 14 2016 * * @brief Computation of mesh intersection with segments * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_SEGMENT_INTERSECTOR_HH__ #define __AKANTU_MESH_SEGMENT_INTERSECTOR_HH__ #include "aka_common.hh" #include "mesh_geom_intersector.hh" #include "mesh_geom_common.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -/// Here, we know what kernel we have to use -typedef Cartesian K; - -template -class MeshSegmentIntersector : public MeshGeomIntersector, K::Segment_3, K> { +template +class MeshSegmentIntersector + : public MeshGeomIntersector, + cgal::Cartesian::Segment_3, cgal::Cartesian> { + using K = cgal::Cartesian; /// Parent class type typedef MeshGeomIntersector, K::Segment_3, K> parent_type; /// Result of intersection function type typedef typename IntersectionTypeHelper, K>, K::Segment_3>::intersection_type result_type; /// Pair of segments and element id typedef std::pair pair_type; public: /// Construct from mesh explicit MeshSegmentIntersector(Mesh & mesh, Mesh & result_mesh); /// Destructor virtual ~MeshSegmentIntersector(); public: /** * @brief Computes the intersection of the mesh with a segment * * @param query the segment to compute the intersections with the mesh */ virtual void computeIntersectionQuery(const K::Segment_3 & query); /// Compute intersection points between the mesh and a query virtual void computeMeshQueryIntersectionPoint(const K::Segment_3 & query, UInt nb_old_nodes); /// Compute the embedded mesh virtual void buildResultFromQueryList(const std::list & query_list); void setPhysicalName(const std::string & other) { current_physical_name = other; } protected: /// Compute segments from intersection list void computeSegments(const std::list & intersections, std::set & segments, const K::Segment_3 & query); protected: /// Result mesh Mesh & result_mesh; /// Physical name of the current batch of queries std::string current_physical_name; }; } // akantu #include "mesh_segment_intersector_tmpl.hh" #endif // __AKANTU_MESH_SEGMENT_INTERSECTOR_HH__ diff --git a/src/geometry/mesh_sphere_intersector.hh b/src/geometry/mesh_sphere_intersector.hh index 63bcea07b..776bd29c2 100644 --- a/src/geometry/mesh_sphere_intersector.hh +++ b/src/geometry/mesh_sphere_intersector.hh @@ -1,105 +1,119 @@ /** * @file mesh_sphere_intersector.hh * * @author Aurelia Isabel Cuba Ramos * @author Lucas Frerot * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jun 23 2015 * @date last modification: Thu Jan 14 2016 * * @brief Computation of mesh intersection with sphere(s) * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ #define __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ #include "aka_common.hh" #include "mesh_geom_intersector.hh" #include "mesh_geom_common.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -/// Here, we know what kernel we have to use -typedef Spherical SK; +template +class MeshSphereIntersector + : public MeshGeomIntersector, + cgal::Spherical::Sphere_3, cgal::Spherical> { + using SK = cgal::Spherical; + using K = cgal::Cartesian; -template -class MeshSphereIntersector : public MeshGeomIntersector, SK::Sphere_3, SK> { /// Parent class type - typedef MeshGeomIntersector, SK::Sphere_3, SK> parent_type; + typedef MeshGeomIntersector, SK::Sphere_3, SK> + parent_type; /// Result of intersection function type - typedef typename IntersectionTypeHelper, K>, K::Segment_3>::intersection_type result_type; + typedef typename IntersectionTypeHelper, K>, + K::Segment_3>::intersection_type + result_type; /// Pair of intersection points and element id typedef std::pair pair_type; public: /// Construct from mesh explicit MeshSphereIntersector(Mesh & mesh); /// Destructor virtual ~MeshSphereIntersector(); public: /// Construct the primitive tree object virtual void constructData(GhostType ghost_type = _not_ghost); /** * @brief Computes the intersection of the mesh with a sphere * * @param query (sphere) to compute the intersections with the mesh */ - virtual void computeIntersectionQuery(const SK::Sphere_3 & query){ - AKANTU_DEBUG_ERROR("This function is not implemented for spheres (It was to generic and has been replaced by computeMeshQueryIntersectionPoint"); + virtual void computeIntersectionQuery(const SK::Sphere_3 & /*query*/) { + AKANTU_DEBUG_ERROR("This function is not implemented for spheres (It was " + "to generic and has been replaced by " + "computeMeshQueryIntersectionPoint"); } - /// Compute intersection points between the mesh primitives (segments) and a query (surface in 3D or a curve in 2D), double intersection points for the same primitives are not considered. A maximum is set to the number of intersection nodes per element: 2 in 2D and 4 in 3D - virtual void computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query, UInt nb_old_nodes); + /// Compute intersection points between the mesh primitives (segments) and a + /// query (surface in 3D or a curve in 2D), double intersection points for the + /// same primitives are not considered. A maximum is set to the number of + /// intersection nodes per element: 2 in 2D and 4 in 3D + virtual void computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query, + UInt nb_old_nodes); /// Build the IGFEM mesh - virtual void buildResultFromQueryList(const std::list & query){ - AKANTU_DEBUG_ERROR("This function is no longer implemented to split geometrical operations and dedicated result construction"); + virtual void + buildResultFromQueryList(const std::list & /*query*/) { + AKANTU_DEBUG_ERROR("This function is no longer implemented to split " + "geometrical operations and dedicated result " + "construction"); } /// Set the tolerance void setToleranceIntersectionOnNode(UInt tol) { this->tol_intersection_on_node = tol; } protected: - /// tolerance for which the intersection is considered on the mesh node (relative to the segment lenght) + /// tolerance for which the intersection is considered on the mesh node + /// (relative to the segment lenght) Real tol_intersection_on_node; - }; - + } // akantu #include "mesh_sphere_intersector_tmpl.hh" #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_HH__ diff --git a/src/geometry/mesh_sphere_intersector_tmpl.hh b/src/geometry/mesh_sphere_intersector_tmpl.hh index ff815fd40..ab6dc705a 100644 --- a/src/geometry/mesh_sphere_intersector_tmpl.hh +++ b/src/geometry/mesh_sphere_intersector_tmpl.hh @@ -1,196 +1,195 @@ /** * @file mesh_sphere_intersector_tmpl.hh * * @author Aurelia Isabel Cuba Ramos * @author Lucas Frerot * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jun 23 2015 * @date last modification: Thu Jan 14 2016 * * @brief Computation of mesh intersection with spheres * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ #define __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ #include "aka_common.hh" #include "mesh_geom_common.hh" #include "tree_type_helper.hh" #include "mesh_sphere_intersector.hh" -#include "static_communicator.hh" namespace akantu { template MeshSphereIntersector::MeshSphereIntersector(Mesh & mesh): parent_type(mesh), tol_intersection_on_node(1e-10) { #if defined(AKANTU_IGFEM) if( (type == _triangle_3) || (type == _igfem_triangle_4) || (type == _igfem_triangle_5) ){ const_cast(this->nb_seg_by_el) = 3; } else { AKANTU_DEBUG_ERROR("Not ready for mesh type " << type); } #else if( (type != _triangle_3) ) AKANTU_DEBUG_ERROR("Not ready for mesh type " << type); #endif // initialize the intersection pointsss array with the spatial dimension this->intersection_points = new Array(0,dim); // A maximum is set to the number of intersection nodes per element to limit the size of new_node_per_elem: 2 in 2D and 4 in 3D this->new_node_per_elem = new Array(0, 1 + 4 * (dim-1)); } template MeshSphereIntersector::~MeshSphereIntersector() { delete this->new_node_per_elem; delete this->intersection_points; } template void MeshSphereIntersector::constructData(GhostType ghost_type) { this->new_node_per_elem->resize(this->mesh.getNbElement(type, ghost_type)); this->new_node_per_elem->clear(); MeshGeomIntersector, SK::Sphere_3, SK>::constructData(ghost_type); } template void MeshSphereIntersector:: computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query, UInt nb_old_nodes) { /// function to replace computeIntersectionQuery in a more generic geometry module version // The newNodeEvent is not send from this method who only compute the intersection points AKANTU_DEBUG_IN(); Array & nodes = this->mesh.getNodes(); UInt nb_node = nodes.size() + this->intersection_points->size(); // Tolerance for proximity checks should be defined by user Real global_tolerance = Math::getTolerance(); Math::setTolerance(tol_intersection_on_node); typedef boost::variant sk_inter_res; - TreeTypeHelper, Spherical>::const_iterator + TreeTypeHelper, cgal::Spherical>::const_iterator it = this->factory.getPrimitiveList().begin(), end= this->factory.getPrimitiveList().end(); for (; it != end ; ++it) { // loop on the primitives (segments) std::list s_results; CGAL::intersection(*it, query, std::back_inserter(s_results)); if (s_results.size() == 1) { // just one point if (pair_type * pair = boost::get(&s_results.front())) { if (pair->second == 1) { // not a point tangent to the sphere // the intersection point written as a vector Vector new_node(dim, 0.0); - Cartesian::Point_3 point(CGAL::to_double(pair->first.x()), + cgal::Cartesian::Point_3 point(CGAL::to_double(pair->first.x()), CGAL::to_double(pair->first.y()), CGAL::to_double(pair->first.z())); for (UInt i = 0 ; i < dim ; i++) { new_node(i) = point[i]; } /// boolean to decide wheter intersection point is on a standard node of the mesh or not bool is_on_mesh = false; /// boolean to decide if this intersection point has been already computed for a neighbor element bool is_new = true; /// check if intersection point has already been computed UInt n = nb_old_nodes; // check if we already compute this intersection and add it as a node for a neighboor element of another type Array::vector_iterator existing_node = nodes.begin(dim); for (; n < nodes.size() ; ++n) {// loop on the nodes from nb_old_nodes if (Math::are_vector_equal(dim, new_node.storage(), existing_node[n].storage())) { is_new = false; break; } } if(is_new){ Array::vector_iterator intersection_points_it = this->intersection_points->begin(dim); Array::vector_iterator intersection_points_end = this->intersection_points->end(dim); for (; intersection_points_it != intersection_points_end ; ++intersection_points_it, ++n) { if (Math::are_vector_equal(dim, new_node.storage(), intersection_points_it->storage())) { is_new = false; break; } } } // get the initial and final points of the primitive (segment) and write them as vectors - Cartesian::Point_3 source_cgal(CGAL::to_double(it->source().x()), + cgal::Cartesian::Point_3 source_cgal(CGAL::to_double(it->source().x()), CGAL::to_double(it->source().y()), CGAL::to_double(it->source().z())); - Cartesian::Point_3 target_cgal(CGAL::to_double(it->target().x()), + cgal::Cartesian::Point_3 target_cgal(CGAL::to_double(it->target().x()), CGAL::to_double(it->target().y()), CGAL::to_double(it->target().z())); Vector source(dim), target(dim); for (UInt i = 0 ; i < dim ; i++) { source(i) = source_cgal[i]; target(i) = target_cgal[i]; } // Check if we are close from a node of the primitive (segment) if (Math::are_vector_equal(dim, source.storage(), new_node.storage()) || Math::are_vector_equal(dim, target.storage(), new_node.storage())) { is_on_mesh = true; is_new = false; } if (is_new) {// if the intersection point is a new one add it to the list this->intersection_points->push_back(new_node); nb_node++; } // deduce the element id UInt element_id = it->id(); // fill the new_node_per_elem array if (!is_on_mesh) { // if the node is not on a mesh node UInt & nb_new_nodes_per_el = (*this->new_node_per_elem)(element_id, 0); nb_new_nodes_per_el += 1; AKANTU_DEBUG_ASSERT(2 * nb_new_nodes_per_el < this->new_node_per_elem->getNbComponent(), "You might have to interface crossing the same material"); (*this->new_node_per_elem)(element_id, (2 * nb_new_nodes_per_el) - 1) = n; (*this->new_node_per_elem)(element_id, 2 * nb_new_nodes_per_el) = it->segId(); } } } } } Math::setTolerance(global_tolerance); AKANTU_DEBUG_OUT(); } } // akantu #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__ diff --git a/src/geometry/tree_type_helper.hh b/src/geometry/tree_type_helper.hh index 853bc63ff..4f20fa558 100644 --- a/src/geometry/tree_type_helper.hh +++ b/src/geometry/tree_type_helper.hh @@ -1,111 +1,112 @@ /** * @file tree_type_helper.hh * * @author Lucas Frerot * * @date creation: Fri Jan 04 2013 * @date last modification: Thu Jan 14 2016 * * @brief Converts element types of a mesh to CGAL primitive types * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TREE_TYPE_HELPER_HH__ #define __AKANTU_TREE_TYPE_HELPER_HH__ #include "aka_common.hh" #include "line_arc.hh" #include "triangle.hh" #include "tetrahedron.hh" #include "aabb_primitive.hh" #include "mesh_geom_common.hh" #include #include namespace akantu { /* -------------------------------------------------------------------------- */ /// Replacement class for algorithm that can't use the AABB tree types template struct VoidTree { VoidTree(const iterator & begin, const iterator & end) {} }; /// Helper class used to ease the use of CGAL AABB tree algorithm template struct TreeTypeHelper { static const bool is_valid = false; typedef Primitive primitive_type; typedef typename std::list container_type; typedef typename container_type::iterator iterator; typedef typename container_type::const_iterator const_iterator; typedef typename CGAL::Point_3 point_type; typedef VoidTree tree; }; /// Helper class used to ease the use of intersections template struct IntersectionTypeHelper; - /** * Macro used to specialize TreeTypeHelper * @param my_primitive associated primitive type * @param my_query query_type * @param my_kernel kernel type */ -#define TREE_TYPE_HELPER_MACRO(my_primitive, my_query, my_kernel) \ - template<> \ - struct TreeTypeHelper, my_kernel> { \ - static const bool is_valid = true; \ - typedef my_primitive primitive_type; \ - typedef my_primitive##_primitive aabb_primitive_type; \ - typedef CGAL::Point_3 point_type; \ - \ - typedef std::list container_type; \ +#define TREE_TYPE_HELPER_MACRO(my_primitive, my_query, my_kernel) \ + template <> struct TreeTypeHelper, my_kernel> { \ + static const bool is_valid = true; \ + typedef my_primitive primitive_type; \ + typedef my_primitive##_primitive aabb_primitive_type; \ + typedef CGAL::Point_3 point_type; \ + \ + typedef std::list container_type; \ typedef container_type::iterator iterator; \ - typedef CGAL::AABB_traits aabb_traits_type; \ - typedef CGAL::AABB_tree tree; \ - typedef tree::Primitive_id id_type; \ - }; \ - \ - template<> \ - struct IntersectionTypeHelper, my_kernel>, my_query> { \ - typedef boost::optional< \ - TreeTypeHelper, my_kernel>::tree::Intersection_and_primitive_id::Type \ - > intersection_type; \ + typedef CGAL::AABB_traits \ + aabb_traits_type; \ + typedef CGAL::AABB_tree tree; \ + typedef tree::Primitive_id id_type; \ + }; \ + \ + template <> \ + struct IntersectionTypeHelper< \ + TreeTypeHelper, my_kernel>, my_query> { \ + typedef boost::optional, \ + my_kernel>::tree::Intersection_and_primitive_id::Type> \ + intersection_type; \ } -TREE_TYPE_HELPER_MACRO(Triangle, Cartesian::Segment_3, Cartesian); -//TREE_TYPE_HELPER_MACRO(Line_arc, Spherical::Sphere_3, Spherical); +TREE_TYPE_HELPER_MACRO(Triangle, cgal::Cartesian::Segment_3, cgal::Cartesian); +//TREE_TYPE_HELPER_MACRO(Line_arc, cgal::Spherical::Sphere_3, cgal::Spherical); #undef TREE_TYPE_HELPER_MACRO } // akantu #endif // __AKANTU_TREE_TYPE_HELPER_HH__ diff --git a/src/model/model.hh b/src/model/model.hh index b265a1dc6..fc8c01c91 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,367 +1,373 @@ /** * @file model.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Oct 16 2015 * * @brief Interface of a model * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_named_argument.hh" #include "fe_engine.hh" #include "mesh.hh" #include "model_options.hh" #include "model_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_HH__ #define __AKANTU_MODEL_HH__ namespace akantu { class SynchronizerRegistry; class Parser; class DumperIOHelper; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { class Model : public Memory, public ModelSolver, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Model(Mesh & mesh, const ModelType & type, UInt spatial_dimension = _all_dimensions, const ID & id = "model", const MemoryID & memory_id = 0); ~Model() override; using FEEngineMap = std::map>; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: virtual void initFullImpl(const ModelOptions & options); public: template std::enable_if_t::value> initFull(pack &&... _pack) { switch (this->model_type) { #ifdef AKANTU_SOLID_MECHANICS case ModelType::_solid_mechanics_model: this->initFullImpl(SolidMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_COHESIVE_ELEMENT case ModelType::_solid_mechanics_model_cohesive: this->initFullImpl(SolidMechanicsModelCohesiveOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_HEAT_TRANSFER case ModelType::_heat_transfer_model: this->initFullImpl(HeatTransferModelOptions{ use_named_args, std::forward(_pack)...}); break; +#endif +#ifdef AKANTU_EMBEDDED + case ModelType::_embedded_model: + this->initFullImpl(EmbeddedInterfaceModelOptions{ + use_named_args, std::forward(_pack)...}); + break; #endif default: this->initFullImpl(ModelOptions{use_named_args, std::forward(_pack)...}); } } template std::enable_if_t::value> initFull(pack &&... _pack) { this->initFullImpl(std::forward(_pack)...); } /// initialize a new solver if needed void initNewSolver(const AnalysisMethod & method); protected: /// get some default values for derived classes virtual std::tuple getDefaultSolverID(const AnalysisMethod & method) = 0; virtual void initModel() = 0; virtual void initFEEngineBoundary(); // /// change local equation number so that PBC is assembled properly // void changeLocalEquationNumberForPBC(std::map & pbc_pair, // UInt dimension); /// function to print the containt of the class void printself(std::ostream &, int = 0) const override{}; // /// initialize the model for PBC // void setPBC(UInt x, UInt y, UInt z); // void setPBC(SurfacePairList & surface_pairs); public: virtual void initPBC(); /// set the parser to use // void setParser(Parser & parser); /* ------------------------------------------------------------------------ */ /* Access to the dumpable interface of the boundaries */ /* ------------------------------------------------------------------------ */ /// Dump the data for a given group void dumpGroup(const std::string & group_name); void dumpGroup(const std::string & group_name, const std::string & dumper_name); /// Dump the data for all boundaries void dumpGroup(); /// Set the directory for a given group void setGroupDirectory(const std::string & directory, const std::string & group_name); /// Set the directory for all boundaries void setGroupDirectory(const std::string & directory); /// Set the base name for a given group void setGroupBaseName(const std::string & basename, const std::string & group_name); /// Get the internal dumper of a given group DumperIOHelper & getGroupDumper(const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Function for non local capabilities */ /* ------------------------------------------------------------------------ */ virtual void updateDataForNonLocalCriterion(__attribute__((unused)) ElementTypeMapReal & criterion) { AKANTU_DEBUG_TO_IMPLEMENT(); } protected: template void allocNodalField(Array *& array, UInt nb_component, const ID & name); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get id of model AKANTU_GET_MACRO(ID, id, const ID) /// get the number of surfaces AKANTU_GET_MACRO(Mesh, mesh, Mesh &) /// synchronize the boundary in case of parallel run virtual void synchronizeBoundaries(){}; /// return the fem object associated with a provided name inline FEEngine & getFEEngine(const ID & name = "") const; /// return the fem boundary object associated with a provided name virtual FEEngine & getFEEngineBoundary(const ID & name = ""); /// register a fem object associated with name template inline void registerFEEngineObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name inline void unRegisterFEEngineObject(const std::string & name); /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); /// return the fem object associated with a provided name template inline FEEngineClass & getFEEngineClass(std::string name = "") const; /// return the fem boundary object associated with a provided name template inline FEEngineClass & getFEEngineClassBoundary(std::string name = ""); /// get the pbc pairs std::map & getPBCPairs() { return pbc_pair; }; /// returns if node is slave in pbc inline bool isPBCSlaveNode(const UInt) const { throw; return false; /* TODO repair PBC*/ } /// returns the array of pbc slave nodes (boolean information) AKANTU_GET_MACRO(IsPBCSlaveNode, is_pbc_slave_node, const Array &) /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /* ------------------------------------------------------------------------ */ /* Pack and unpack helper functions */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbIntegrationPoints(const Array & elements, const ID & fem_id = ID()) const; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ void setTextModeToDumper(); virtual void addDumpGroupFieldToDumper(const std::string & field_id, dumper::Field * field, DumperIOHelper & dumper); virtual void addDumpField(const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); virtual void setBaseName(const std::string & basename); virtual void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, const ElementKind & element_kind, bool padding_flag); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual dumper::Field * createNodalFieldReal(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } virtual dumper::Field * createNodalFieldUInt(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } virtual dumper::Field * createNodalFieldBool(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } virtual dumper::Field * createElementalField(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const UInt & spatial_dimension, __attribute__((unused)) const ElementKind & kind) { return nullptr; } void setDirectory(const std::string & directory); void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory); virtual void dump(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend std::ostream & operator<<(std::ostream &, const Model &); /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// Mesh Mesh & mesh; /// Spatial dimension of the problem UInt spatial_dimension; /// the main fem object present in all models FEEngineMap fems; /// the fem object present in all models for boundaries FEEngineMap fems_boundary; /// default fem object std::string default_fem; /// pbc pairs std::map pbc_pair; /// flag per node to know is pbc slave Array is_pbc_slave_node; /// parser to the pointer to use Parser & parser; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "model_inline_impl.cc" #endif /* __AKANTU_MODEL_HH__ */ diff --git a/src/model/model_options.hh b/src/model/model_options.hh index 3c1e7c96b..a4d1d8f58 100644 --- a/src/model/model_options.hh +++ b/src/model/model_options.hh @@ -1,108 +1,137 @@ /** * @file model_options.hh * * @author Nicolas Richart * * @date creation Mon Dec 04 2017 * * @brief A Documented file. * * @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_named_argument.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_OPTIONS_HH__ #define __AKANTU_MODEL_OPTIONS_HH__ namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(analysis_method); } struct ModelOptions { explicit ModelOptions(AnalysisMethod analysis_method = _static) : analysis_method(analysis_method) {} template ModelOptions(use_named_args_t, pack &&... _pack) : ModelOptions(OPTIONAL_NAMED_ARG(analysis_method, _static)) {} virtual ~ModelOptions() = default; AnalysisMethod analysis_method; }; #ifdef AKANTU_SOLID_MECHANICS /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelOptions : public ModelOptions { explicit SolidMechanicsModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass) : ModelOptions(analysis_method) {} template SolidMechanicsModelOptions(use_named_args_t, pack &&... _pack) : SolidMechanicsModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} }; #endif /* -------------------------------------------------------------------------- */ #ifdef AKANTU_COHESIVE_ELEMENT namespace { DECLARE_NAMED_ARGUMENT(is_extrinsic); } /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelCohesiveOptions : public SolidMechanicsModelOptions { SolidMechanicsModelCohesiveOptions( AnalysisMethod analysis_method = _explicit_lumped_mass, bool extrinsic = false) : SolidMechanicsModelOptions(analysis_method), extrinsic(extrinsic) {} template SolidMechanicsModelCohesiveOptions(use_named_args_t, pack &&... _pack) : SolidMechanicsModelCohesiveOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass), OPTIONAL_NAMED_ARG(is_extrinsic, false)) {} bool extrinsic{false}; }; #endif #ifdef AKANTU_HEAT_TRANSFER /* -------------------------------------------------------------------------- */ struct HeatTransferModelOptions : public ModelOptions { explicit HeatTransferModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass) : ModelOptions(analysis_method) {} template HeatTransferModelOptions(use_named_args_t, pack &&... _pack) : HeatTransferModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} }; #endif +#ifdef AKANTU_EMBEDDED + +namespace { + DECLARE_NAMED_ARGUMENT(init_intersections); +} +/* -------------------------------------------------------------------------- */ +struct EmbeddedInterfaceModelOptions : SolidMechanicsModelOptions { + /** + * @brief Constructor for EmbeddedInterfaceModelOptions + * @param analysis_method see SolidMechanicsModelOptions + * @param init_intersections compute intersections + */ + EmbeddedInterfaceModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass, + bool init_intersections = true): + SolidMechanicsModelOptions(analysis_method), + has_intersections(init_intersections) + {} + + template + EmbeddedInterfaceModelOptions(use_named_args_t, pack &&... _pack) + : EmbeddedInterfaceModelOptions( + OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass), + OPTIONAL_NAMED_ARG(init_intersections, true)) {} + + /// Should consider reinforcements + bool has_intersections; +}; +#endif + } // akantu #endif /* __AKANTU_MODEL_OPTIONS_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh b/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh index 2d139f5b7..9c8ad1ea4 100644 --- a/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh +++ b/src/model/solid_mechanics/materials/material_embedded/embedded_internal_field.hh @@ -1,83 +1,84 @@ /** * @file embedded_internal_field.hh * * @author Lucas Frerot * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jun 30 2015 * * @brief Embedded Material internal properties * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_EMBEDDED_INTERNAL_FIELD_HH__ #define __AKANTU_EMBEDDED_INTERNAL_FIELD_HH__ namespace akantu { class Material; class FEEngine; /// This class is used for MaterialReinforcement internal fields template class EmbeddedInternalField : public InternalField { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Constructor EmbeddedInternalField(const ID & id, Material & material): InternalField(id, material, material.getModel().getFEEngine("EmbeddedInterfaceFEEngine"), material.getElementFilter()) { this->spatial_dimension = 1; } /// Copy constructor EmbeddedInternalField(const ID & id, const EmbeddedInternalField & other): InternalField(id, other) { this->spatial_dimension = 1; } void operator=(const EmbeddedInternalField & other) { InternalField::operator=(other); this->spatial_dimension = 1; } }; /// Method used to initialise the embedded internal fields from material file -template<> -inline void ParsableParamTyped< EmbeddedInternalField >::parseParam(const ParserParameter & in_param) { - ParsableParam::parseParam(in_param); +template <> +inline void ParameterTyped>::setAuto( + const ParserParameter & in_param) { + Parameter::setAuto(in_param); Real r = in_param; param.setDefaultValue(r); } } // akantu #endif // __AKANTU_EMBEDDED_INTERNAL_FIELD_HH__ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc index 3fd3dce3b..def0887e8 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc @@ -1,784 +1,766 @@ /** * @file material_reinforcement.cc * * @author Lucas Frerot * * @date creation: Wed Mar 25 2015 * @date last modification: Tue Dec 08 2015 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_voigthelper.hh" #include "material_reinforcement.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialReinforcement::MaterialReinforcement(SolidMechanicsModel & model, - UInt spatial_dimension, + UInt /*spatial_dimension*/, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Material(model, dim, mesh, fe_engine, id), // /!\ dim, not spatial_dimension ! model(NULL), stress_embedded("stress_embedded", *this, 1, fe_engine, this->element_filter), gradu_embedded("gradu_embedded", *this, 1, fe_engine, this->element_filter), directing_cosines("directing_cosines", *this, 1, fe_engine, this->element_filter), pre_stress("pre_stress", *this, 1, fe_engine, this->element_filter), area(1.0), shape_derivatives() { AKANTU_DEBUG_IN(); this->initialize(model); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initialize(SolidMechanicsModel & a_model) { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&a_model); AKANTU_DEBUG_ASSERT(this->model != NULL, "MaterialReinforcement needs an EmbeddedInterfaceModel"); this->registerParam("area", area, _pat_parsable | _pat_modifiable, "Reinforcement cross-sectional area"); this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable, "Uniform pre-stress"); this->unregisterInternal(this->stress); // Fool the AvgHomogenizingFunctor // stress.initialize(dim * dim); // Reallocate the element filter - this->element_filter.free(); - this->model->getInterfaceMesh().initElementTypeMapArray(this->element_filter, - 1, 1, false, _ek_regular); + this->element_filter.initialize(this->model->getInterfaceMesh(), + _spatial_dimension = 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialReinforcement::~MaterialReinforcement() { AKANTU_DEBUG_IN(); ElementTypeMap *>::type_iterator it = shape_derivatives.firstType(), end = shape_derivatives.lastType(); for (; it != end; ++it) { delete shape_derivatives(*it, _not_ghost); delete shape_derivatives(*it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initMaterial() { Material::initMaterial(); stress_embedded.initialize(dim * dim); gradu_embedded.initialize(dim * dim); pre_stress.initialize(1); /// We initialise the stuff that is not going to change during the simulation this->allocBackgroundShapeDerivatives(); this->initBackgroundShapeDerivatives(); this->initDirectingCosines(); } /* -------------------------------------------------------------------------- */ /** * Background shape derivatives need to be stored per background element * types but also per embedded element type, which is why they are stored * in an ElementTypeMap *>. The outer ElementTypeMap * refers to the embedded types, and the inner refers to the background types. */ template void MaterialReinforcement::allocBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); Mesh & mesh = model->getMesh(); ghost_type_t::iterator int_ghost_it = ghost_type_t::begin(); // Loop over interface ghosts for (; int_ghost_it != ghost_type_t::end(); ++int_ghost_it) { Mesh::type_iterator interface_type_it = interface_mesh.firstType(1, *int_ghost_it); Mesh::type_iterator interface_type_end = interface_mesh.lastType(1, *int_ghost_it); // Loop over interface types for (; interface_type_it != interface_type_end; ++interface_type_it) { Mesh::type_iterator background_type_it = mesh.firstType(dim, *int_ghost_it); Mesh::type_iterator background_type_end = mesh.lastType(dim, *int_ghost_it); // Loop over background types for (; background_type_it != background_type_end ; ++background_type_it) { const ElementType & int_type = *interface_type_it; const ElementType & back_type = *background_type_it; const GhostType & int_ghost = *int_ghost_it; std::string shaped_id = "embedded_shape_derivatives"; if (int_ghost == _ghost) shaped_id += ":ghost"; ElementTypeMapArray * shaped_etma = new ElementTypeMapArray(shaped_id, this->name); UInt nb_points = Mesh::getNbNodesPerElement(back_type); UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbIntegrationPoints(int_type); UInt nb_elements = element_filter(int_type, int_ghost).size(); // Alloc the background ElementTypeMapArray shaped_etma->alloc(nb_elements * nb_quad_points, dim * nb_points, back_type); // Insert the background ElementTypeMapArray in the interface // ElementTypeMap shape_derivatives(shaped_etma, int_type, int_ghost); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); Mesh & mesh = model->getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, _not_ghost); Mesh::type_iterator type_end = mesh.lastType(dim, _not_ghost); for (; type_it != type_end; ++type_it) { computeBackgroundShapeDerivatives(*type_it, _not_ghost); // computeBackgroundShapeDerivatives(*type_it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initDirectingCosines() { AKANTU_DEBUG_IN(); Mesh & mesh = model->getInterfaceMesh(); Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost); const UInt voigt_size = getTangentStiffnessVoigtSize(dim); directing_cosines.initialize(voigt_size * voigt_size); for (; type_it != type_end; ++type_it) { computeDirectingCosines(*type_it, _not_ghost); // computeDirectingCosines(*type_it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end; ++type_it) { assembleStiffnessMatrix(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end; ++type_it) { this->assembleInternalForces(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeGradU(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine") .getNbIntegrationPoints(type); Array & gradu_vec = gradu_embedded(type, ghost_type); Mesh::type_iterator back_it = model->getMesh().firstType(dim, ghost_type); Mesh::type_iterator back_end = model->getMesh().lastType(dim, ghost_type); for (; back_it != back_end; ++back_it) { UInt nodes_per_background_e = Mesh::getNbNodesPerElement(*back_it); Array & shapesd = shape_derivatives(type, ghost_type)->operator()(*back_it, ghost_type); Array * background_filter = new Array(nb_element, 1, "background_filter"); filterInterfaceBackgroundElements(*background_filter, *back_it, type, ghost_type); Array * disp_per_element = new Array(0, dim * nodes_per_background_e, "disp_elem"); FEEngine::extractNodalToElementField(model->getMesh(), model->getDisplacement(), *disp_per_element, *back_it, ghost_type, *background_filter); Array::matrix_iterator disp_it = disp_per_element->begin(dim, nodes_per_background_e); Array::matrix_iterator disp_end = disp_per_element->end(dim, nodes_per_background_e); Array::matrix_iterator shapes_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator grad_u_it = gradu_vec.begin(dim, dim); for (; disp_it != disp_end; ++disp_it) { for (UInt i = 0; i < nb_quad_points; i++, ++shapes_it, ++grad_u_it) { Matrix & B = *shapes_it; Matrix & du = *grad_u_it; Matrix & u = *disp_it; du.mul(u, B); } } delete background_filter; delete disp_per_element; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = model->getInterfaceMesh().firstType(); Mesh::type_iterator last_type = model->getInterfaceMesh().lastType(); for (; it != last_type; ++it) { computeGradU(*it, ghost_type); computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleInternalForces( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = model->getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end; ++type_it) { assembleInternalForcesInterface(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** - * Computes and assemble the residual. Residual in reinforcement is computed as - * : + * Computes and assemble the residual. Residual in reinforcement is computed as: + * * \f[ * \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s} * \f] */ template void MaterialReinforcement::assembleInternalForcesInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = getTangentStiffnessVoigtSize(dim); - Array & residual = const_cast &>(model->getResidual()); - FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); - FEEngine & background_engine = model->getFEEngine(); Array & elem_filter = element_filter(interface_type, ghost_type); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt back_dof = dim * nodes_per_background_e; - Array & shapesd = shape_derivatives(interface_type, ghost_type) - -> - operator()(background_type, ghost_type); + Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( + background_type, ghost_type); - Array * integrant = new Array(nb_quadrature_points * nb_element, - back_dof, - "integrant"); + Array integrant(nb_quadrature_points * nb_element, back_dof, + "integrant"); - Array::vector_iterator integrant_it = integrant->begin(back_dof); - Array::vector_iterator integrant_end = integrant->end(back_dof); + Array::vector_iterator integrant_it = integrant.begin(back_dof); + Array::vector_iterator integrant_end = integrant.end(back_dof); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator C_it = directing_cosines(interface_type, ghost_type) .begin(voigt_size, voigt_size); Array::matrix_iterator sigma_it = stress_embedded(interface_type, ghost_type).begin(dim, dim); Vector sigma(voigt_size); Matrix Bvoigt(voigt_size, back_dof); Vector Ct_sigma(voigt_size); for (; integrant_it != integrant_end ; ++integrant_it, ++B_it, ++C_it, ++sigma_it) { VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, nodes_per_background_e); Matrix & C = *C_it; Vector & BtCt_sigma = *integrant_it; stressTensorToVoigtVector(*sigma_it, sigma); Ct_sigma.mul(C, sigma); BtCt_sigma.mul(Bvoigt, Ct_sigma); BtCt_sigma *= area; } - Array * residual_interface = - new Array(nb_element, back_dof, "residual_interface"); - interface_engine.integrate(*integrant, *residual_interface, back_dof, + Array residual_interface(nb_element, back_dof, "residual_interface"); + interface_engine.integrate(integrant, residual_interface, back_dof, interface_type, ghost_type, elem_filter); + integrant.resize(0); - delete integrant; - - Array * background_filter = - new Array(nb_element, 1, "background_filter"); + Array background_filter(nb_element, 1, "background_filter"); - filterInterfaceBackgroundElements(*background_filter, background_type, + filterInterfaceBackgroundElements(background_filter, background_type, interface_type, ghost_type); - background_engine.assembleArray( - *residual_interface, residual, - model->getDOFSynchronizer().getLocalDOFEquationNumbers(), dim, - background_type, ghost_type, *background_filter, -1.0); - delete residual_interface; - delete background_filter; + + model->getDOFManager().assembleElementalArrayLocalArray( + residual_interface, model->getInternalForce(), background_type, ghost_type, -1., + background_filter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::filterInterfaceBackgroundElements( Array & filter, const ElementType & type, const ElementType & interface_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); filter.resize(0); filter.clear(); Array & elements = model->getInterfaceAssociatedElements(interface_type, ghost_type); Array & elem_filter = element_filter(interface_type, ghost_type); Array::scalar_iterator filter_it = elem_filter.begin(), filter_end = elem_filter.end(); for (; filter_it != filter_end; ++filter_it) { Element & elem = elements(*filter_it); if (elem.type == type) filter.push_back(elem.element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeDirectingCosines( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = this->model->getInterfaceMesh(); const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const UInt steel_dof = dim * nb_nodes_per_element; const UInt voigt_size = getTangentStiffnessVoigtSize(dim); const UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine") .getNbIntegrationPoints(type, ghost_type); Array node_coordinates(this->element_filter(type, ghost_type).size(), steel_dof); this->model->getFEEngine().template extractNodalToElementField(interface_mesh, interface_mesh.getNodes(), node_coordinates, type, ghost_type, this->element_filter(type, ghost_type)); Array::matrix_iterator directing_cosines_it = directing_cosines(type, ghost_type).begin(voigt_size, voigt_size); Array::matrix_iterator node_coordinates_it = node_coordinates.begin(dim, nb_nodes_per_element); Array::matrix_iterator node_coordinates_end = node_coordinates.end(dim, nb_nodes_per_element); for (; node_coordinates_it != node_coordinates_end; ++node_coordinates_it) { for (UInt i = 0; i < nb_quad_points; i++, ++directing_cosines_it) { Matrix & nodes = *node_coordinates_it; Matrix & cosines = *directing_cosines_it; computeDirectingCosinesOnQuad(nodes, cosines); } } // Mauro: the directing_cosines internal is defined on the quadrature points of each element AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = model->getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end; ++type_it) { assembleStiffnessMatrixInterface(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001) * \f[ * \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T * \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}} * \f] */ template void MaterialReinforcement::assembleStiffnessMatrixInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = getTangentStiffnessVoigtSize(dim); - SparseMatrix & K = const_cast(model->getStiffnessMatrix()); - - FEEngine & background_engine = model->getFEEngine(); FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); Array & elem_filter = element_filter(interface_type, ghost_type); Array & grad_u = gradu_embedded(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt back_dof = dim * nodes_per_background_e; UInt integrant_size = back_dof; grad_u.resize(nb_quadrature_points * nb_element); - Array * tangent_moduli = new Array( - nb_element * nb_quadrature_points, 1, "interface_tangent_moduli"); - tangent_moduli->clear(); - computeTangentModuli(interface_type, *tangent_moduli, ghost_type); + Array tangent_moduli(nb_element * nb_quadrature_points, 1, + "interface_tangent_moduli"); + computeTangentModuli(interface_type, tangent_moduli, ghost_type); - Array & shapesd = shape_derivatives(interface_type, ghost_type) - -> - operator()(background_type, ghost_type); + Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( + background_type, ghost_type); - Array * integrant = new Array(nb_element * nb_quadrature_points, - integrant_size * integrant_size, - "B^t*C^t*D*C*B"); - integrant->clear(); + Array integrant(nb_element * nb_quadrature_points, + integrant_size * integrant_size, "B^t*C^t*D*C*B"); /// Temporary matrices for integrant product Matrix Bvoigt(voigt_size, back_dof); Matrix DC(voigt_size, voigt_size); Matrix DCB(voigt_size, back_dof); Matrix CtDCB(voigt_size, back_dof); - Array::scalar_iterator D_it = tangent_moduli->begin(); - Array::scalar_iterator D_end = tangent_moduli->end(); + Array::scalar_iterator D_it = tangent_moduli.begin(); + Array::scalar_iterator D_end = tangent_moduli.end(); Array::matrix_iterator C_it = directing_cosines(interface_type, ghost_type) .begin(voigt_size, voigt_size); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator integrant_it = - integrant->begin(integrant_size, integrant_size); + integrant.begin(integrant_size, integrant_size); for (; D_it != D_end; ++D_it, ++C_it, ++B_it, ++integrant_it) { Real & D = *D_it; Matrix & C = *C_it; Matrix & B = *B_it; Matrix & BtCtDCB = *integrant_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(B, Bvoigt, nodes_per_background_e); DC.clear(); DC(0, 0) = D * area; DC *= C; DCB.mul(DC, Bvoigt); CtDCB.mul(C, DCB); BtCtDCB.mul(Bvoigt, CtDCB); } - delete tangent_moduli; + tangent_moduli.resize(0); - Array * K_interface = new Array( + Array K_interface( nb_element, integrant_size * integrant_size, "K_interface"); - interface_engine.integrate(*integrant, *K_interface, + interface_engine.integrate(integrant, K_interface, integrant_size * integrant_size, interface_type, ghost_type, elem_filter); - delete integrant; + integrant.resize(0); // Mauro: Here K_interface contains the local stiffness matrices, // directing_cosines contains the information about the orientation // of the reinforcements, any rotation of the local stiffness matrix // can be done here - Array * background_filter = new Array(nb_element, 1, "background_filter"); + Array background_filter(nb_element, 1, "background_filter"); - filterInterfaceBackgroundElements(*background_filter, background_type, + filterInterfaceBackgroundElements(background_filter, background_type, interface_type, ghost_type); - background_engine.assembleMatrix(*K_interface, K, dim, background_type, - ghost_type, *background_filter); + model->getDOFManager().assembleElementalMatricesToMatrix( + "K", "displacement", K_interface, background_type, ghost_type, + _symmetric, background_filter); - delete K_interface; - delete background_filter; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// In this function, type and ghost_type refer to background elements template void MaterialReinforcement::computeBackgroundShapeDerivatives( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); FEEngine & engine = model->getFEEngine(); FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); Mesh::type_iterator interface_type = interface_mesh.firstType(); Mesh::type_iterator interface_last = interface_mesh.lastType(); for (; interface_type != interface_last; ++interface_type) { Array & filter = element_filter(*interface_type, ghost_type); const UInt nb_elements = filter.size(); const UInt nb_nodes = Mesh::getNbNodesPerElement(type); const UInt nb_quad_per_element = interface_engine.getNbIntegrationPoints(*interface_type); Array quad_pos(nb_quad_per_element * nb_elements, dim, "interface_quad_points"); quad_pos.resize(nb_quad_per_element * nb_elements); interface_engine.interpolateOnIntegrationPoints( interface_mesh.getNodes(), quad_pos, dim, *interface_type, ghost_type, filter); Array & background_shapesd = shape_derivatives(*interface_type, ghost_type) -> operator()(type, ghost_type); background_shapesd.clear(); Array * background_elements = new Array( nb_elements, 1, "computeBackgroundShapeDerivatives:background_filter"); filterInterfaceBackgroundElements(*background_elements, type, *interface_type, ghost_type); Array::scalar_iterator back_it = background_elements->begin(), back_end = background_elements->end(); Array::matrix_iterator shapesd_it = background_shapesd.begin(dim, nb_nodes); Array::vector_iterator quad_pos_it = quad_pos.begin(dim); for (; back_it != back_end; ++back_it) { for (UInt i = 0; i < nb_quad_per_element; i++, ++shapesd_it, ++quad_pos_it) engine.computeShapeDerivatives(*quad_pos_it, *back_it, type, *shapesd_it, ghost_type); } delete background_elements; } AKANTU_DEBUG_OUT(); } -template Real MaterialReinforcement::getEnergy(std::string id) { +template +Real MaterialReinforcement::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); if (id == "potential") { Real epot = 0.; computePotentialEnergyByElements(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension), end = element_filter.lastType(spatial_dimension); for (; it != end; ++it) { FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); epot += interface_engine.integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); epot *= area; } return epot; } AKANTU_DEBUG_OUT(); return 0; } // /* -------------------------------------------------------------------------- // */ // template // ElementTypeMap MaterialReinforcement::getInternalDataPerElem(const // ID & field_name, // const // ElementKind // & // kind, // const // ID & // fe_engine_id) // const // { // return Material::getInternalDataPerElem(field_name, kind, // "EmbeddedInterfaceFEEngine"); // } // /* -------------------------------------------------------------------------- // */ // // Author is Guillaume Anciaux, see material.cc // template // void MaterialReinforcement::flattenInternal(const std::string & // field_id, // ElementTypeMapArray & // internal_flat, // const GhostType ghost_type, // ElementKind element_kind) // const { // AKANTU_DEBUG_IN(); // if (field_id == "stress_embedded" || field_id == "inelastic_strain") { // Material::flattenInternalIntern(field_id, // internal_flat, // 1, // ghost_type, // _ek_not_defined, // &(this->element_filter), // &(this->model->getInterfaceMesh())); // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ -INSTANTIATE_MATERIAL(MaterialReinforcement); +INSTANTIATE_MATERIAL_ONLY(MaterialReinforcement); /* -------------------------------------------------------------------------- */ } // akantu diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh index fd8c6b52e..09925013b 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh @@ -1,200 +1,192 @@ /** * @file material_reinforcement.hh * * @author Lucas Frerot * * @date creation: Fri Mar 13 2015 * @date last modification: Tue Nov 24 2015 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_REINFORCEMENT_HH__ #define __AKANTU_MATERIAL_REINFORCEMENT_HH__ #include "aka_common.hh" #include "material.hh" #include "embedded_interface_model.hh" #include "embedded_internal_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /** * @brief Material used to represent embedded reinforcements * * This class is used for computing the reinforcement stiffness matrix * along with the reinforcement residual. Room is made for constitutive law, * but actual use of contitutive laws is made in MaterialReinforcementTemplate. * * Be careful with the dimensions in this class : * - this->spatial_dimension is always 1 * - the template parameter dim is the dimension of the problem */ template class MaterialReinforcement : virtual public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Constructor MaterialReinforcement(SolidMechanicsModel & model, UInt spatial_dimension, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor - virtual ~MaterialReinforcement(); + ~MaterialReinforcement() override; protected: void initialize(SolidMechanicsModel & a_model); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Init the material - virtual void initMaterial(); + void initMaterial() override; /// Init the background shape derivatives void initBackgroundShapeDerivatives(); /// Init the cosine matrices void initDirectingCosines(); /// Assemble stiffness matrix - virtual void assembleStiffnessMatrix(GhostType ghost_type); - - /// Update the residual - virtual void updateResidual(GhostType ghost_type = _not_ghost); - - /// Assembled the residual - virtual void assembleResidual(GhostType ghost_type); + void assembleStiffnessMatrix(GhostType ghost_type) override; /// Compute all the stresses ! - virtual void computeAllStresses(GhostType ghost_type); - - /// Compute the stiffness parameter for elements of a type - virtual void computeTangentModuli(const ElementType & type, - Array & tangent, - GhostType ghost_type) = 0; + void computeAllStresses(GhostType ghost_type) override; /// Compute energy - virtual Real getEnergy(std::string id); + Real getEnergy(const std::string & id) override; + + /// Assemble the residual of one type of element (typically _segment_2) + void assembleInternalForces(GhostType ghost_type) override; // virtual ElementTypeMap getInternalDataPerElem(const ID & field_name, // const ElementKind & // kind, // const ID & // fe_engine_id) const; // /// Reimplementation of Material's function to accomodate for interface // mesh // virtual void flattenInternal(const std::string & field_id, // ElementTypeMapArray & internal_flat, // const GhostType ghost_type = _not_ghost, // ElementKind element_kind = _ek_not_defined) // const; /* ------------------------------------------------------------------------ */ /* Protected methods */ /* ------------------------------------------------------------------------ */ protected: /// Allocate the background shape derivatives void allocBackgroundShapeDerivatives(); /// Compute the directing cosines matrix for one element type void computeDirectingCosines(const ElementType & type, GhostType ghost_type); /// Compute the directing cosines matrix on quadrature points. inline void computeDirectingCosinesOnQuad(const Matrix & nodes, Matrix & cosines); /// Assemble the stiffness matrix for an element type (typically _segment_2) void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// Assemble the stiffness matrix for background & interface types void assembleStiffnessMatrixInterface(const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type); /// Compute the background shape derivatives for a type void computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type); /// Filter elements crossed by interface of a type void filterInterfaceBackgroundElements(Array & filter, const ElementType & type, const ElementType & interface_type, GhostType ghost_type); /// Assemble the residual of one type of element (typically _segment_2) void assembleInternalForces(const ElementType & type, GhostType ghost_type); /// Assemble the residual for a pair of elements void assembleInternalForcesInterface(const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type); // TODO figure out why voigt size is 4 in 2D inline void stressTensorToVoigtVector(const Matrix & tensor, Vector & vector); inline void strainTensorToVoigtVector(const Matrix & tensor, Vector & vector); /// Compute gradu on the interface quadrature points virtual void computeGradU(const ElementType & type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Embedded model EmbeddedInterfaceModel * model; /// Stress in the reinforcement InternalField stress_embedded; /// Gradu of concrete on reinforcement InternalField gradu_embedded; /// C matrix on quad InternalField directing_cosines; /// Prestress on quad InternalField pre_stress; /// Cross-sectional area Real area; /// Background mesh shape derivatives ElementTypeMap *> shape_derivatives; }; #include "material_reinforcement_inline_impl.cc" } // akantu #endif // __AKANTU_MATERIAL_REINFORCEMENT_HH__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 7cded9c8a..f015fb6d5 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,562 +1,562 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Model of Solid 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 . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials - void + virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::force vector (external forces) AKANTU_GET_MACRO(Force, *external_force, Array &); /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); AKANTU_SET_MACRO(MaterialSelector, material_selector, std::shared_ptr); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// lumped mass array Array * mass{nullptr}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force{nullptr}; /// internal forces array Array * internal_force{nullptr}; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs{nullptr}; /// array of current position used during update residual Array * current_position{nullptr}; UInt current_position_release{0}; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; #ifdef SWIGPYTHON public: #endif /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; #ifdef SWIGPYTHON protected: #endif /// class defining of to choose a material std::shared_ptr material_selector; /// flag defining if the increment must be computed or not bool increment_flag; /// tells if the material are instantiated bool are_materials_instantiated; using flatten_internal_map = std::map, ElementTypeMapArray *>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.hh b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.hh index 889a63c4e..60492c373 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.hh @@ -1,91 +1,93 @@ /** * @file embedded_interface_intersector.hh * * @author Lucas Frerot * * @date creation: Fri May 01 2015 * @date last modification: Mon Dec 14 2015 * * @brief Class that loads the interface from mesh and computes intersections * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_EMBEDDED_INTERFACE_INTERSECTOR_HH__ #define __AKANTU_EMBEDDED_INTERFACE_INTERSECTOR_HH__ #include "aka_common.hh" #include "mesh_geom_common.hh" #include "mesh_geom_abstract.hh" #include "mesh_segment_intersector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -typedef Cartesian K; +namespace { + using K = cgal::Cartesian; +} /** * @brief Computes the intersections of the reinforcements defined in the primitive mesh * * The purpose of this class is to look for reinforcements in the primitive mesh, which * should be defined by physical groups with the same names as the reinforcement materials * in the model. * * It then constructs the CGAL primitives from the elements of those reinforcements * and computes the intersections with the background mesh, to create an `interface_mesh`, * which is in turn used by the EmbeddedInterfaceModel. * * @see MeshSegmentIntersector, MeshGeomAbstract * @see EmbeddedInterfaceModel */ class EmbeddedInterfaceIntersector : public MeshGeomAbstract { public: /// Construct from mesh and a reinforcement mesh explicit EmbeddedInterfaceIntersector(Mesh & mesh, const Mesh & primitive_mesh); /// Destructor virtual ~EmbeddedInterfaceIntersector(); public: /// Generate the interface mesh virtual void constructData(GhostType ghost_type = _not_ghost); /// Create a segment with an element connectivity K::Segment_3 createSegment(const Vector & connectivity); /// Getter for interface mesh AKANTU_GET_MACRO_NOT_CONST(InterfaceMesh, interface_mesh, Mesh &); protected: /// Resulting mesh of intersection Mesh interface_mesh; /// Mesh used for primitive construction const Mesh & primitive_mesh; }; } // akantu #endif // __AKANTU_EMBEDDED_INTERFACE_INTERSECTOR_HH__ diff --git a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc index 910af71ff..cd14e54f4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc @@ -1,210 +1,149 @@ /** * @file embedded_interface_model.cc * * @author Lucas Frerot * * @date creation: Fri Mar 13 2015 * @date last modification: Mon Dec 14 2015 * * @brief Model of Solid Mechanics with embedded interfaces * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ -#include "aka_common.hh" - #include "embedded_interface_model.hh" #include "material_reinforcement.hh" +#include "mesh_iterators.hh" +#include "integrator_gauss.hh" +#include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER -# include "dumper_paraview.hh" +# include "dumper_iohelper_paraview.hh" # include "dumpable_inline_impl.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { -const EmbeddedInterfaceModelOptions - default_embedded_interface_model_options(_explicit_lumped_mass, false, false); - /* -------------------------------------------------------------------------- */ EmbeddedInterfaceModel::EmbeddedInterfaceModel(Mesh & mesh, Mesh & primitive_mesh, UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, spatial_dimension, id, memory_id), intersector(mesh, primitive_mesh), interface_mesh(NULL), primitive_mesh(primitive_mesh), interface_material_selector(NULL) { + this->model_type = ModelType::_embedded_model; + // This pointer should be deleted by ~SolidMechanicsModel() - MaterialSelector * mat_sel_pointer = - new MeshDataMaterialSelector("physical_names", *this); + auto mat_sel_pointer = + std::make_shared>("physical_names", *this); - this->setMaterialSelector(*mat_sel_pointer); + this->setMaterialSelector(mat_sel_pointer); interface_mesh = &(intersector.getInterfaceMesh()); // Create 1D FEEngine on the interface mesh registerFEEngineObject("EmbeddedInterfaceFEEngine", *interface_mesh, 1); } /* -------------------------------------------------------------------------- */ EmbeddedInterfaceModel::~EmbeddedInterfaceModel() { delete interface_material_selector; } /* -------------------------------------------------------------------------- */ -void EmbeddedInterfaceModel::initFull(const ModelOptions & options) { - const EmbeddedInterfaceModelOptions & eim_options = - dynamic_cast(options); - - // We don't want to initiate materials before shape functions are initialized - SolidMechanicsModelOptions dummy_options(eim_options.analysis_method, true); +void EmbeddedInterfaceModel::initFullImpl(const ModelOptions & options) { + const auto & eim_options = + dynamic_cast(options); // Do no initialize interface_mesh if told so - if (!eim_options.no_init_intersections) + if (eim_options.has_intersections) intersector.constructData(); - // Initialize interface FEEngine - FEEngine & engine = getFEEngine("EmbeddedInterfaceFEEngine"); - engine.initShapeFunctions(_not_ghost); - engine.initShapeFunctions(_ghost); - - SolidMechanicsModel::initFull(dummy_options); - - // This will call SolidMechanicsModel::initMaterials() last - this->initMaterials(); + SolidMechanicsModel::initFullImpl(options); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("reinforcement", id); this->mesh.addDumpMeshToDumper("reinforcement", *interface_mesh, 1, _not_ghost, _ek_regular); #endif } -/* -------------------------------------------------------------------------- */ -// This function is very similar to SolidMechanicsModel's -void EmbeddedInterfaceModel::initMaterials() { - Element element; +void EmbeddedInterfaceModel::initModel() { + // Initialize interface FEEngine + FEEngine & engine = getFEEngine("EmbeddedInterfaceFEEngine"); + engine.initShapeFunctions(_not_ghost); + engine.initShapeFunctions(_ghost); +} +/* -------------------------------------------------------------------------- */ +void EmbeddedInterfaceModel::assignMaterialToElements( + const ElementTypeMapArray * filter) { delete interface_material_selector; interface_material_selector = - new InterfaceMeshDataMaterialSelector("physical_names", *this); - - for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { - element.ghost_type = *gt; - - Mesh::type_iterator it = interface_mesh->firstType(1, *gt); - Mesh::type_iterator end = interface_mesh->lastType(1, *gt); - - for (; it != end ; ++it) { - UInt nb_element = interface_mesh->getNbElement(*it, *gt); - - element.type = *it; - - Array & mat_indexes = material_index.alloc(nb_element, 1, *it, *gt); - - for (UInt el = 0 ; el < nb_element ; el++) { - element.element = el; - UInt mat_index = (*interface_material_selector)(element); - - AKANTU_DEBUG_ASSERT(mat_index < materials.size(), - "The material selector returned an index that does not exist"); - mat_indexes(element.element) = mat_index; - materials.at(mat_index)->addElement(*it, el, *gt); - } - } - } - - SolidMechanicsModel::initMaterials(); + new InterfaceMeshDataMaterialSelector("physical_names", + *this); + + for_each_elements(mesh, + [&](auto && element) { + auto mat_index = (*interface_material_selector)(element); + material_index(element) = mat_index; + materials[mat_index]->addElement(element); + }, + _element_filter = filter); + + SolidMechanicsModel::assignMaterialToElements(filter); } -// /** -// * DO NOT REMOVE - This prevents the material reinforcement to register -// * their number of components. Problems arise with AvgHomogenizingFunctor -// * if the material reinforcement gives its number of components for a -// * field. Since AvgHomogenizingFunctor verifies that all the fields -// * have the same number of components, an exception is raised. -// */ -// ElementTypeMap EmbeddedInterfaceModel::getInternalDataPerElem(const std::string & field_name, -// const ElementKind & kind) { -// if (!(this->isInternal(field_name,kind))) AKANTU_EXCEPTION("unknown internal " << field_name); - -// for (UInt m = 0; m < materials.size() ; ++m) { -// if (materials[m]->isInternal(field_name, kind)) { -// Material * mat = NULL; - -// switch(this->spatial_dimension) { -// case 1: -// mat = dynamic_cast *>(materials[m]); -// break; - -// case 2: -// mat = dynamic_cast *>(materials[m]); -// break; - -// case 3: -// mat = dynamic_cast *>(materials[m]); -// break; -// } - -// if (mat == NULL && field_name != "stress_embedded") -// return materials[m]->getInternalDataPerElem(field_name,kind); -// else if (mat != NULL && field_name == "stress_embedded") -// return mat->getInternalDataPerElem(field_name, kind, "EmbeddedInterfaceFEEngine"); -// } -// } - -// return ElementTypeMap(); -// } - /* -------------------------------------------------------------------------- */ void EmbeddedInterfaceModel::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { #ifdef AKANTU_USE_IOHELPER dumper::Field * field = NULL; // If dumper is reinforcement, create a 1D elemental field if (dumper_name == "reinforcement") field = this->createElementalField(field_id, group_name, padding_flag, 1, element_kind); else { try { SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, element_kind, padding_flag); } catch (...) {} } if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name,group_name); Model::addDumpGroupFieldToDumper(field_id,field,dumper); } #endif } } // akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.hh b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.hh index 39edb360d..fab1f94fe 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.hh @@ -1,167 +1,150 @@ /** * @file embedded_interface_model.hh * * @author Lucas Frerot * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Dec 14 2015 * * @brief Model of Solid Mechanics with embedded interfaces * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_EMBEDDED_INTERFACE_MODEL_HH__ #define __AKANTU_EMBEDDED_INTERFACE_MODEL_HH__ #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "mesh.hh" #include "embedded_interface_intersector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -/// Options for the EmbeddedInterfaceModel -struct EmbeddedInterfaceModelOptions : SolidMechanicsModelOptions { - /** - * @brief Constructor for EmbeddedInterfaceModelOptions - * @param analysis_method see SolidMeechanicsModelOptions - * @param no_init_intersections ignores the embedded elements - * @param no_init_materials see SolidMeechanicsModelOptions - */ - EmbeddedInterfaceModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass, - bool no_init_intersections = false, - bool no_init_materials = false): - SolidMechanicsModelOptions(analysis_method, no_init_materials), - no_init_intersections(no_init_intersections) - {} - - /// Ignore the reinforcements - bool no_init_intersections; -}; - -extern const EmbeddedInterfaceModelOptions default_embedded_interface_model_options; - /** * @brief Solid mechanics model using the embedded model. * * This SolidMechanicsModel subclass implements the embedded model, * a method used to represent 1D elements in a finite elements model * (eg. reinforcements in concrete). * * In addition to the SolidMechanicsModel properties, this model has * a mesh of the 1D elements embedded in the model, and an instance of the * EmbeddedInterfaceIntersector class for the computation of the intersections of the * 1D elements with the background (bulk) mesh. * * @see MaterialReinforcement */ class EmbeddedInterfaceModel : public SolidMechanicsModel { - /// Same type as SolidMechanicsModel - typedef FEEngineTemplate MyFEEngineType; + using MyFEEngineType = SolidMechanicsModel::MyFEEngineType; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /** * @brief Constructor * * @param mesh main mesh (concrete) * @param primitive_mesh mesh of the embedded reinforcement */ EmbeddedInterfaceModel(Mesh & mesh, Mesh & primitive_mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "embedded_interface_model", const MemoryID & memory_id = 0); /// Destructor - virtual ~EmbeddedInterfaceModel(); + ~EmbeddedInterfaceModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Initialise the model - virtual void initFull(const ModelOptions & options = default_embedded_interface_model_options); + void initFullImpl( + const ModelOptions & options = EmbeddedInterfaceModelOptions()) override; /// Initialise the materials - virtual void initMaterials(); + void + assignMaterialToElements(const ElementTypeMapArray * filter) override; + + /// Initialize the embedded shape functions + void initModel() override; /// Allows filtering of dump fields which need to be dumpes on interface mesh - virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, - const std::string & field_id, - const std::string & group_name, - const ElementKind & element_kind, - bool padding_flag); + void addDumpGroupFieldToDumper(const std::string & dumper_name, + const std::string & field_id, + const std::string & group_name, + const ElementKind & element_kind, + bool padding_flag) override; // virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, // const ElementKind & kind); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get interface mesh AKANTU_GET_MACRO(InterfaceMesh, *interface_mesh, Mesh &); /// Get associated elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(InterfaceAssociatedElements, interface_mesh->getData("associated_element"), Element); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Intersector object to build the interface mesh EmbeddedInterfaceIntersector intersector; /// Interface mesh (weak reference) Mesh * interface_mesh; /// Mesh used to create the CGAL primitives for intersections Mesh & primitive_mesh; /// Material selector for interface MaterialSelector * interface_material_selector; }; /// Material selector based on mesh data for interface elements template class InterfaceMeshDataMaterialSelector : public ElementDataMaterialSelector { public: InterfaceMeshDataMaterialSelector(const std::string & name, const EmbeddedInterfaceModel & model, UInt first_index = 1) : ElementDataMaterialSelector(model.getInterfaceMesh().getData(name), model, first_index) {} }; } // akantu #endif // __AKANTU_EMBEDDED_INTERFACE_MODEL_HH__ diff --git a/test/test_geometry/test_geometry_intersection.cc b/test/test_geometry/test_geometry_intersection.cc index 2466b614a..7b29e4962 100644 --- a/test/test_geometry/test_geometry_intersection.cc +++ b/test/test_geometry/test_geometry_intersection.cc @@ -1,127 +1,127 @@ /** * @file test_geometry_intersection.cc * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Feb 27 2015 * @date last modification: Thu Jan 14 2016 * * @brief Tests the intersection module * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh_geom_factory.hh" #include "tree_type_helper.hh" #include "geom_helper_functions.hh" #include "mesh_geom_common.hh" #include #include #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; -typedef Cartesian K; +typedef cgal::Cartesian K; typedef IntersectionTypeHelper, K>, K::Segment_3>::intersection_type result_type; -typedef Spherical SK; +typedef cgal::Spherical SK; typedef boost::variant > sk_inter_res; /*typedef CGAL::cpp11::result_of >)>::type sk_res;*/ typedef std::pair pair_type; /* -------------------------------------------------------------------------- */ int main (int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblWarning); Mesh mesh(2); mesh.read("test_geometry_triangle.msh"); MeshGeomFactory<2, _triangle_3, Triangle, K> factory(mesh); factory.constructData(); const TreeTypeHelper, K>::tree & tree = factory.getTree(); K::Point_3 a(0., 0.25, 0.), b(1., 0.25, 0.); K::Segment_3 line(a, b); K::Point_3 begin(a), intermediate(0.25, 0.25, 0.), end(0.75, 0.25, 0.); K::Segment_3 result_0(begin, intermediate), result_1(intermediate, end); std::list list_of_intersections; tree.all_intersections(line, std::back_inserter(list_of_intersections)); const result_type & intersection_0 = list_of_intersections.back(); const result_type & intersection_1 = list_of_intersections.front(); if (!intersection_0 || !intersection_1) return EXIT_FAILURE; /// *-> first is the intersection ; *->second is the primitive id if (const K::Segment_3 * segment = boost::get(&(intersection_0->first))) { if (!compareSegments(*segment, result_0)) { return EXIT_FAILURE; } } else return EXIT_FAILURE; if (const K::Segment_3 * segment = boost::get(&(intersection_1->first))) { if (!compareSegments(*segment, result_1)) { return EXIT_FAILURE; } } else return EXIT_FAILURE; SK::Sphere_3 sphere(SK::Point_3(0, 0, 0), 3.); SK::Segment_3 seg(SK::Point_3(0, 0, 0), SK::Point_3(2., 2., 2.)); SK::Line_arc_3 arc(seg); std::list s_results; CGAL::intersection(arc, sphere, std::back_inserter(s_results)); if (pair_type * pair = boost::get(&s_results.front())) { std::cout << "xi = " << to_double(pair->first.x()) << ", yi = " << to_double(pair->first.y()) << std::endl; if (!comparePoints(pair->first, SK::Circular_arc_point_3(1.0, 1.0, 1.0))) return EXIT_FAILURE; } else return EXIT_FAILURE; MeshGeomFactory<2, _triangle_3, Line_arc, SK> Sfactory(mesh); Sfactory.constructData(); finalize(); return EXIT_SUCCESS; } diff --git a/test/test_geometry/test_geometry_predicates.cc b/test/test_geometry/test_geometry_predicates.cc index 91675ae0c..2cee26181 100644 --- a/test/test_geometry/test_geometry_predicates.cc +++ b/test/test_geometry/test_geometry_predicates.cc @@ -1,88 +1,88 @@ /** * @file test_geometry_predicates.cc * * @author Lucas Frerot * * @date creation: Fri Jan 04 2013 * @date last modification: Thu Jan 14 2016 * * @brief Tests the geometry predicates * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "geom_helper_functions.hh" #include "mesh_geom_common.hh" #include /* -------------------------------------------------------------------------- */ using namespace akantu; -typedef Cartesian K; +typedef cgal::Cartesian K; typedef K::Point_3 Point; typedef K::Segment_3 Segment; int main(int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblWarning); Point a(0, 1, 0); Point b(0, 1, 1); Segment seg1(a, b); Segment seg2(b, a); if (!compareSegments(seg1, seg2)) return EXIT_FAILURE; // Testing sort + unique on list of segments std::vector > pair_list; pair_list.push_back(std::make_pair(seg1, 1)); pair_list.push_back(std::make_pair(seg2, 2)); segmentPairsLess sorter; std::sort(pair_list.begin(), pair_list.end(), sorter); std::vector >::iterator it = std::unique(pair_list.begin(), pair_list.end(), compareSegmentPairs); if (it - pair_list.begin() != 1) { std::cout << pair_list.size() << std::endl; return EXIT_FAILURE; } // Testing insertion in set std::set, segmentPairsLess> pair_set; pair_set.insert(pair_set.begin(), std::make_pair(seg1, 1)); pair_set.insert(pair_set.begin(), std::make_pair(seg2, 2)); if (pair_set.size() != 1) { std::cout << pair_set.size() << std::endl; return EXIT_FAILURE; } finalize(); return EXIT_SUCCESS; } diff --git a/test/test_geometry/test_segment_intersection_tetrahedron_4.cc b/test/test_geometry/test_segment_intersection_tetrahedron_4.cc index 77343b936..62ce9b434 100644 --- a/test/test_geometry/test_segment_intersection_tetrahedron_4.cc +++ b/test/test_geometry/test_segment_intersection_tetrahedron_4.cc @@ -1,148 +1,148 @@ /** * @file test_segment_intersection_tetrahedron_4.cc * * @author Lucas Frerot * * @date creation: Fri Feb 27 2015 * @date last modification: Thu Jan 14 2016 * * @brief Tests the intersection module with _tetrahedron_4 elements * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh_segment_intersector.hh" #include "mesh_geom_common.hh" #include /* -------------------------------------------------------------------------- */ using namespace akantu; -typedef Cartesian K; +typedef cgal::Cartesian K; typedef K::Point_3 Point; typedef K::Segment_3 Segment; /* -------------------------------------------------------------------------- */ int main (int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblError); Mesh mesh(3), interface_mesh(3, "interface_mesh"); mesh.read("test_geometry_tetrahedron.msh"); MeshSegmentIntersector<3, _tetrahedron_4> intersector(mesh, interface_mesh); intersector.constructData(); // Testing a segment going through the cube Point point(1., 1., 1.); Segment segment(CGAL::ORIGIN, point); intersector.computeIntersectionQuery(segment); std::cout << "number of seg_2 : " << interface_mesh.getNbElement(_segment_2) << std::endl; if (interface_mesh.getNbElement(_segment_2) != 2) return EXIT_FAILURE; Vector bary(2), bary1(2), bary2(2); Element test; test.type = _segment_2; test.element = 0; interface_mesh.getBarycenter(test, bary1); test.element = 1; interface_mesh.getBarycenter(test, bary2); Real first_bary[] = {1./6., 1./6., 1./6.}; Real second_bary[] = {2./3., 2./3., 2./3.}; // We don't know the order of the elements, so here we test permutations if (!( (Math::are_vector_equal(3, bary1.storage(), first_bary) && Math::are_vector_equal(3, bary2.storage(), second_bary) ) || (Math::are_vector_equal(3, bary1.storage(), second_bary) && Math::are_vector_equal(3, bary2.storage(), first_bary) ) )) return EXIT_FAILURE; // Testing a segment completely inside one element Point a(0.05, 0.05, 0.05), b(0.06, 0.06, 0.06); Segment inside_segment(a, b); intersector.computeIntersectionQuery(inside_segment); test.element = interface_mesh.getNbElement(_segment_2) - 1; interface_mesh.getBarycenter(test, bary); Real third_bary[] = {0.055, 0.055, 0.055}; if (!Math::are_vector_equal(3, bary.storage(), third_bary)) return EXIT_FAILURE; // Testing a segment whose end points are inside elements Point c(0.1, 0.1, 0.1), d(0.9, 0.9, 0.9); Segment crossing_segment(c, d); intersector.computeIntersectionQuery(crossing_segment); UInt el1 = interface_mesh.getNbElement(_segment_2) - 2; UInt el2 = el1 + 1; test.element = el1; interface_mesh.getBarycenter(test, bary1); test.element = el2; interface_mesh.getBarycenter(test, bary2); Real fourth_bary[] = {13./60., 13./60., 13./60.}; Real fifth_bary[] = {37./60., 37./60., 37./60.}; // We don't know the order of the elements, so here we test permutations if (!( (Math::are_vector_equal(3, bary1.storage(), fourth_bary) && Math::are_vector_equal(3, bary2.storage(), fifth_bary) ) || (Math::are_vector_equal(3, bary1.storage(), fifth_bary) && Math::are_vector_equal(3, bary2.storage(), fourth_bary) ) )) return EXIT_FAILURE; // Testing a segment along the edge of elements Point e(1, 0, 0), f(0, 1, 0); Segment edge_segment(e, f); UInt current_nb_elements = interface_mesh.getNbElement(_segment_2); intersector.computeIntersectionQuery(edge_segment); if (interface_mesh.getNbElement(_segment_2) != current_nb_elements + 1) return EXIT_FAILURE; test.element = interface_mesh.getNbElement(_segment_2) - 1; interface_mesh.getBarycenter(test, bary); Real sixth_bary[] = {0.5, 0.5, 0}; if (!Math::are_vector_equal(3, bary.storage(), sixth_bary)) return EXIT_FAILURE; return EXIT_SUCCESS; } diff --git a/test/test_geometry/test_segment_intersection_triangle_3.cc b/test/test_geometry/test_segment_intersection_triangle_3.cc index af6b1ccb2..787a18f70 100644 --- a/test/test_geometry/test_segment_intersection_triangle_3.cc +++ b/test/test_geometry/test_segment_intersection_triangle_3.cc @@ -1,145 +1,145 @@ /** * @file test_segment_intersection_triangle_3.cc * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Feb 27 2015 * @date last modification: Thu Jan 14 2016 * * @brief Tests the interface mesh generation * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh_segment_intersector.hh" #include "mesh_sphere_intersector.hh" #include "geom_helper_functions.hh" #include "mesh_geom_common.hh" #include /* -------------------------------------------------------------------------- */ using namespace akantu; -typedef Cartesian K; -typedef Spherical SK; +typedef cgal::Cartesian K; +typedef cgal::Spherical SK; /* -------------------------------------------------------------------------- */ int main (int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblError); Math::setTolerance(1e-10); Mesh mesh(2), interface_mesh(2, "interface_mesh"); mesh.read("test_geometry_triangle.msh"); MeshSegmentIntersector<2, _triangle_3> intersector(mesh, interface_mesh); intersector.constructData(); // Testing a segment going out of the mesh K::Point_3 a(0, 0.25, 0), b(1, 0.25, 0), c(0.25, 0, 0), d(0.25, 1, 0); K::Segment_3 h_interface(a, b), v_interface(c, d); std::list interface_list; interface_list.push_back(h_interface); interface_list.push_back(v_interface); intersector.computeIntersectionQueryList(interface_list); if (interface_mesh.getNbElement(_segment_2) != 4) return EXIT_FAILURE; Vector bary(2); Element test; test.element = 0; test.type = _segment_2; interface_mesh.getBarycenter(test, bary); Real first_bary[] = {0.125, 0.25}; if (!Math::are_vector_equal(2, bary.storage(), first_bary)) return EXIT_FAILURE; // Testing a segment completely inside an element K::Point_3 e(0.1, 0.33, 0), f(0.1, 0.67, 0); K::Segment_3 inside_segment(e, f); intersector.computeIntersectionQuery(inside_segment); test.element = interface_mesh.getNbElement(_segment_2) - 1; interface_mesh.getBarycenter(test, bary); Real second_bary[] = {0.1, 0.5}; if (!Math::are_vector_equal(2, bary.storage(), second_bary)) return EXIT_FAILURE; #if 0 - // Spherical kernel testing the addition of nodes + // cgal::Spherical kernel testing the addition of nodes std::cout << "initial mesh size = " << mesh.getNodes().size() << " nodes" << std::endl; SK::Sphere_3 sphere(SK::Point_3(0, 1, 0), 0.2*0.2); SK::Sphere_3 sphere2(SK::Point_3(1, 0, 0), 0.4999999999); MeshSphereIntersector<2, _triangle_3> intersector_sphere(mesh); intersector_sphere.constructData(); std::list sphere_list; sphere_list.push_back(sphere); sphere_list.push_back(sphere2); intersector_sphere.computeIntersectionQueryList(sphere_list); std::cout << "final mesh size = " << mesh.getNodes().size() << std::endl; const Array new_node_triangle_3 = intersector_sphere.getNewNodePerElem(); const Array & nodes = mesh.getNodes(); std::cout << "New nodes :" << std::endl; std::cout << "node 5, x=" << nodes(4,0) << ", y=" << nodes(4,1) << std::endl; std::cout << "node 6, x=" << nodes(5,0) << ", y=" << nodes(5,1) << std::endl; std::cout << "node 7, x=" << nodes(6,0) << ", y=" << nodes(6,1) << std::endl; if ( (new_node_triangle_3(0,0) != 1) || (new_node_triangle_3(1,0) != 2)){ for(UInt k=0; k != new_node_triangle_3.size(); ++k){ std::cout << new_node_triangle_3(k,0) << " new nodes in element " << k << ", node(s): " << new_node_triangle_3(k,1) << ", " << new_node_triangle_3(k,3) << ", on segment(s):" << new_node_triangle_3(k,2) << ", " << new_node_triangle_3(k,4) << std::endl; } return EXIT_FAILURE; } #endif finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc index b3176dddf..6ad371486 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_element_matrix.cc @@ -1,90 +1,87 @@ /** * @file test_embedded_element_matrix.cc * * @author Lucas Frerot * * @date creation: Wed Mar 25 2015 * @date last modification: Wed May 13 2015 * * @brief test of the class EmbeddedInterfaceModel * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ -#include "aka_common.hh" #include "embedded_interface_model.hh" +#include "sparse_matrix_aij.hh" using namespace akantu; int main(int argc, char * argv[]) { debug::setDebugLevel(dblWarning); initialize("embedded_element.dat", argc, argv); - const UInt dim = 2; + constexpr UInt dim = 2; + constexpr ElementType type = _segment_2; const Real height = 0.5; Mesh mesh(dim); mesh.read("triangle.msh"); - Array nodes_vec(2, dim, "reinforcement_nodes"); - nodes_vec.storage()[0] = 0; nodes_vec.storage()[1] = height; - nodes_vec.storage()[2] = 1; nodes_vec.storage()[3] = height; + Mesh reinforcement_mesh(dim, "reinforcement_mesh"); + auto & nodes = reinforcement_mesh.getNodes(); + nodes.push_back(Vector({0, height})); + nodes.push_back(Vector({1, height})); - Array conn_vec(1, 2, "reinforcement_connectivity"); - conn_vec.storage()[0] = 0; conn_vec.storage()[1] = 1; + reinforcement_mesh.addConnectivityType(type); + auto & connectivity = reinforcement_mesh.getConnectivity(type); + connectivity.push_back(Vector({0, 1})); Array names_vec(1, 1, "reinforcement", "reinforcement_names"); - - Mesh reinforcement_mesh(dim, "reinforcement_mesh"); - reinforcement_mesh.getNodes().copy(nodes_vec); - reinforcement_mesh.addConnectivityType(_segment_2); - reinforcement_mesh.getConnectivity(_segment_2).copy(conn_vec); - reinforcement_mesh.registerData("physical_names").alloc(1, 1, _segment_2); - reinforcement_mesh.getData("physical_names")(_segment_2).copy(names_vec); + reinforcement_mesh.registerData("physical_names").alloc(1, 1, type); + reinforcement_mesh.getData("physical_names")(type).copy(names_vec); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); - model.initFull(EmbeddedInterfaceModelOptions(_static)); + model.initFull(_analysis_method = _static); - if (model.getInterfaceMesh().getNbElement(_segment_2) != 1) + if (model.getInterfaceMesh().getNbElement(type) != 1) return EXIT_FAILURE; if (model.getInterfaceMesh().getSpatialDimension() != 2) return EXIT_FAILURE; model.assembleStiffnessMatrix(); - SparseMatrix & K = model.getStiffnessMatrix(); + SparseMatrixAIJ & K = + dynamic_cast(model.getDOFManager().getMatrix("K")); Math::setTolerance(1e-8); // Testing the assembled stiffness matrix if (!Math::are_float_equal(K(0, 0), 1. - height) || !Math::are_float_equal(K(0, 2), height - 1.) || !Math::are_float_equal(K(2, 0), height - 1.) || !Math::are_float_equal(K(2, 2), 1. - height)) return EXIT_FAILURE; - model.updateResidual(); - return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc index 26798992c..fd4603895 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc @@ -1,106 +1,103 @@ /** * @file test_embedded_interface_model.cc * * @author Lucas Frerot * * @date creation: Wed Mar 25 2015 * @date last modification: Thu Jul 09 2015 * * @brief Embedded model test based on potential energy * * @section LICENSE * * Copyright (©) 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 . * */ #include #include "aka_common.hh" #include "embedded_interface_model.hh" +#include "sparse_matrix.hh" using namespace akantu; int main (int argc, char * argv[]) { debug::setDebugLevel(dblWarning); initialize("material.dat", argc, argv); UInt dim = 2; Math::setTolerance(1e-7); // Mesh here is a 1x1 patch Mesh mesh(dim); mesh.read("embedded_mesh.msh"); Array nodes_vec(2, dim, "reinforcement_nodes"); nodes_vec.storage()[0] = 0; nodes_vec.storage()[1] = 0.5; nodes_vec.storage()[2] = 1; nodes_vec.storage()[3] = 0.5; Array conn_vec(1, 2, "reinforcement_connectivity"); conn_vec.storage()[0] = 0; conn_vec.storage()[1] = 1; Array names_vec(1, 1, "reinforcement", "reinforcement_names"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); reinforcement_mesh.getNodes().copy(nodes_vec); reinforcement_mesh.addConnectivityType(_segment_2); reinforcement_mesh.getConnectivity(_segment_2).copy(conn_vec); reinforcement_mesh.registerData("physical_names").alloc(1, 1, _segment_2); reinforcement_mesh.getData("physical_names")(_segment_2).copy(names_vec); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); - model.initFull(EmbeddedInterfaceModelOptions(_static)); + model.initFull(_analysis_method = _static); Array & nodes = mesh.getNodes(); Array & forces = model.getForce(); Array & bound = model.getBlockedDOFs(); forces(2, 0) = -250; forces(5, 0) = -500; forces(8, 0) = -250; for (UInt i = 0 ; i < mesh.getNbNodes() ; i++) { if (Math::are_float_equal(nodes(i, 0), 0.)) bound(i, 0) = true; if (Math::are_float_equal(nodes(i, 1), 0.)) bound(i, 1) = true; } model.addDumpFieldVector("displacement"); model.addDumpFieldTensor("stress"); model.setBaseNameToDumper("reinforcement", "reinforcement"); model.addDumpFieldTensorToDumper("reinforcement", "stress_embedded"); - // Assemble the global stiffness matrix - model.assembleStiffnessMatrix(); - model.updateResidual(); + model.solveStep(); - model.getStiffnessMatrix().saveMatrix("matrix_test"); + model.getDOFManager().getMatrix("K").saveMatrix("matrix_test"); - model.solveStatic<_scm_newton_raphson_tangent_not_computed, _scc_residual>(1e-7, 1); - model.updateResidual(); model.dump(); Real pot_energy = model.getEnergy("potential"); if (std::abs(pot_energy - 7.37343e-06) > 1e-5) return EXIT_FAILURE; finalize(); return 0; } diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc index 7f24996e2..76a80f295 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc @@ -1,225 +1,222 @@ /** * @file test_embedded_interface_model_prestress.cc * * @author Lucas Frerot * * @date creation: Tue Apr 28 2015 * @date last modification: Thu Oct 15 2015 * * @brief Embedded model test for prestressing (bases on stress norm) * * @section LICENSE * * Copyright (©) 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 . * */ #include "aka_common.hh" #include "embedded_interface_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; #define YG 0.483644859 #define I_eq 0.012488874 #define A_eq (1e-2 + 1. / 7. * 1.) /* -------------------------------------------------------------------------- */ struct StressSolution : public BC::Neumann::FromHigherDim { Real M; Real I; Real yg; Real pre_stress; StressSolution(UInt dim, Real M, Real I, Real yg = 0, Real pre_stress = 0) : BC::Neumann::FromHigherDim(Matrix(dim, dim)), M(M), I(I), yg(yg), pre_stress(pre_stress) {} virtual ~StressSolution() {} - void operator()(const IntegrationPoint & quad_point, - Vector & dual, - const Vector & coord, - const Vector & normals) const { + void operator()(const IntegrationPoint & /*quad_point*/, Vector & dual, + const Vector & coord, + const Vector & normals) const { UInt dim = coord.size(); if (dim < 2) AKANTU_DEBUG_ERROR("Solution not valid for 1D"); Matrix stress(dim, dim); stress.clear(); stress(0, 0) = this->stress(coord(1)); dual.mul(stress, normals); } inline Real stress(Real height) const { return -M / I * (height - yg) + pre_stress; } inline Real neutral_axis() const { return -I * pre_stress / M + yg; } }; /* -------------------------------------------------------------------------- */ int main (int argc, char * argv[]) { initialize("prestress.dat", argc, argv); debug::setDebugLevel(dblError); Math::setTolerance(1e-6); const UInt dim = 2; /* -------------------------------------------------------------------------- */ Mesh mesh(dim); mesh.read("embedded_mesh_prestress.msh"); mesh.createGroupsFromMeshData("physical_names"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); try { reinforcement_mesh.read("embedded_mesh_prestress_reinforcement.msh"); } catch(debug::Exception & e) {} reinforcement_mesh.createGroupsFromMeshData("physical_names"); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); model.initFull(EmbeddedInterfaceModelOptions(_static)); /* -------------------------------------------------------------------------- */ /* Computation of analytical residual */ /* -------------------------------------------------------------------------- */ /* * q = 1000 N/m * L = 20 m * a = 1 m */ - Real steel_area = model.getMaterial("reinforcement").getParam("area"); + Real steel_area = + model.getMaterial("reinforcement").get("area"); Real pre_stress = 1e6; Real stress_norm = 0.; StressSolution * concrete_stress = NULL, * steel_stress = NULL; Real pre_force = pre_stress * steel_area; Real pre_moment = -pre_force * (YG - 0.25); Real neutral_axis = YG - I_eq / A_eq * pre_force / pre_moment; concrete_stress = new StressSolution(dim, pre_moment, 7. * I_eq, YG, -pre_force / (7. * A_eq)); steel_stress = new StressSolution(dim, pre_moment, I_eq, YG, pre_stress - pre_force / A_eq); stress_norm = std::abs(concrete_stress->stress(1)) * (1 - neutral_axis) * 0.5 + std::abs(concrete_stress->stress(0)) * neutral_axis * 0.5 + std::abs(steel_stress->stress(0.25)) * steel_area; model.applyBC(*concrete_stress, "XBlocked"); ElementGroup & end_node = mesh.getElementGroup("EndNode"); NodeGroup & end_node_group = end_node.getNodeGroup(); NodeGroup::const_node_iterator end_node_it = end_node_group.begin(); Vector end_node_force = model.getForce().begin(dim)[*end_node_it]; end_node_force(0) += steel_stress->stress(0.25) * steel_area; Array analytical_residual(mesh.getNbNodes(), dim, "analytical_residual"); analytical_residual.copy(model.getForce()); model.getForce().clear(); delete concrete_stress; delete steel_stress; /* -------------------------------------------------------------------------- */ model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "XBlocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "YBlocked"); - // Assemble the global stiffness matrix - model.assembleStiffnessMatrix(); - - model.updateResidual(); - - if (!model.solveStatic<_scm_newton_raphson_tangent_not_computed, _scc_residual>(1e-6, 1)) + try { + model.solveStep(); + } catch (debug::Exception e) { + std::cerr << e.what() << std::endl; return EXIT_FAILURE; - - model.updateResidual(); + } /* -------------------------------------------------------------------------- */ /* Computation of FEM residual norm */ /* -------------------------------------------------------------------------- */ ElementGroup & xblocked = mesh.getElementGroup("XBlocked"); NodeGroup & boundary_nodes = xblocked.getNodeGroup(); NodeGroup::const_node_iterator nodes_it = boundary_nodes.begin(), nodes_end = boundary_nodes.end(); - Array::vector_iterator com_res = model.getResidual().begin(dim); + Array::vector_iterator com_res = model.getInternalForce().begin(dim); Array::vector_iterator ana_res = analytical_residual.begin(dim); Array::vector_iterator position = mesh.getNodes().begin(dim); Real res_sum = 0.; UInt lower_node = -1; UInt upper_node = -1; Real lower_dist = 1; Real upper_dist = 1; for (; nodes_it != nodes_end ; ++nodes_it) { UInt node_number = *nodes_it; const Vector res = com_res[node_number]; const Vector ana = ana_res[node_number]; const Vector pos = position[node_number]; if (!Math::are_float_equal(pos(1), 0.25)) { if ((std::abs(pos(1) - 0.25) < lower_dist) && (pos(1) < 0.25)) { lower_dist = std::abs(pos(1) - 0.25); lower_node = node_number; } if ((std::abs(pos(1) - 0.25) < upper_dist) && (pos(1) > 0.25)) { upper_dist = std::abs(pos(1) - 0.25); upper_node = node_number; } } for (UInt i = 0 ; i < dim ; i++) { if (!Math::are_float_equal(pos(1), 0.25)) { res_sum += std::abs(res(i)); } } } const Vector upper_res = com_res[upper_node], lower_res = com_res[lower_node]; const Vector end_node_res = com_res[*end_node_it]; Vector delta = upper_res - lower_res; delta *= lower_dist / (upper_dist + lower_dist); Vector concrete_residual = lower_res + delta; Vector steel_residual = end_node_res - concrete_residual; for (UInt i = 0 ; i < dim ; i++) { res_sum += std::abs(concrete_residual(i)); res_sum += std::abs(steel_residual(i)); } Real relative_error = std::abs(res_sum - stress_norm) / stress_norm; if (relative_error > 1e-3) return EXIT_FAILURE; finalize(); return 0; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc index 565a4cb72..0a7121dc4 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc @@ -1,314 +1,314 @@ /** * @file test_solid_mechanics_model_cube3d.cc * * @author Guillaume Anciaux * * @date creation: Wed Aug 04 2010 * @date last modification: Thu Aug 06 2015 * * @brief test of the class SolidMechanicsModel on the 3d cube * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition_functor.hh" #include "test_solid_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { constexpr Real A = 1e-1; constexpr Real E = 1.; -constexpr Real poisson = .3; -constexpr Real lambda = E * poisson / (1. + poisson) / (1. - 2. * poisson); +constexpr Real poisson = 3. / 10; +constexpr Real lambda = E * poisson / (1 + poisson) / (1 - 2 * poisson); constexpr Real mu = E / 2 / (1. + poisson); constexpr Real rho = 1.; const Real cp = std::sqrt((lambda + 2 * mu) / rho); const Real cs = std::sqrt(mu / rho); const Real c = std::sqrt(E / rho); const Vector k = {.5, 0., 0.}; const Vector psi1 = {0., 0., 1.}; const Vector psi2 = {0., 1., 0.}; const Real knorm = k.norm(); /* -------------------------------------------------------------------------- */ template struct Verification {}; /* -------------------------------------------------------------------------- */ template <> struct Verification<1> { void displacement(Vector & disp, const Vector & coord, Real current_time) { const auto & x = coord(_x); const Real omega = c * k[0]; disp(_x) = A * cos(k[0] * x - omega * current_time); } void velocity(Vector & vel, const Vector & coord, Real current_time) { const auto & x = coord(_x); const Real omega = c * k[0]; vel(_x) = omega * A * sin(k[0] * x - omega * current_time); } }; /* -------------------------------------------------------------------------- */ template <> struct Verification<2> { void displacement(Vector & disp, const Vector & X, Real current_time) { Vector kshear = {k[1], k[0]}; Vector kpush = {k[0], k[1]}; const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; disp = A * (kpush * cos(phase_p) + kshear * cos(phase_s)); } void velocity(Vector & vel, const Vector & X, Real current_time) { Vector kshear = {k[1], k[0]}; Vector kpush = {k[0], k[1]}; const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; vel = A * (kpush * omega_p * cos(phase_p) + kshear * omega_s * cos(phase_s)); } }; /* -------------------------------------------------------------------------- */ template <> struct Verification<3> { void displacement(Vector & disp, const Vector & coord, Real current_time) { const auto & X = coord; Vector kpush = k; Vector kshear1(3); Vector kshear2(3); kshear1.crossProduct(k, psi1); kshear2.crossProduct(k, psi2); const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; disp = A * (kpush * cos(phase_p) + kshear1 * cos(phase_s) + kshear2 * cos(phase_s)); } void velocity(Vector & vel, const Vector & coord, Real current_time) { const auto & X = coord; Vector kpush = k; Vector kshear1(3); Vector kshear2(3); kshear1.crossProduct(k, psi1); kshear2.crossProduct(k, psi2); const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; vel = A * (kpush * omega_p * cos(phase_p) + kshear1 * omega_s * cos(phase_s) + kshear2 * omega_s * cos(phase_s)); } }; /* -------------------------------------------------------------------------- */ template class SolutionFunctor : public BC::Dirichlet::DirichletFunctor { public: SolutionFunctor(Real current_time, SolidMechanicsModel & model) : current_time(current_time), model(model) {} public: static constexpr UInt dim = ElementClass<_type>::getSpatialDimension(); inline void operator()(UInt node, Vector & flags, Vector & primal, const Vector & coord) const { flags(0) = true; auto & vel = model.getVelocity(); auto it = vel.begin(model.getSpatialDimension()); Vector v = it[node]; Verification verif; verif.displacement(primal, coord, current_time); verif.velocity(v, coord, current_time); } private: Real current_time; SolidMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ // This fixture uses somewhat finer meshes representing bars. template class TestSMMFixtureBar : public TestSMMFixture> { using parent = TestSMMFixture>; public: void SetUp() override { this->mesh_file = "bar" + aka::to_string(this->type) + ".msh"; parent::SetUp(); getStaticParser().parse("test_solid_mechanics_model_" "dynamics_material.dat"); auto analysis_method = std::tuple_element_t<1, type_>::value; this->model->initFull(_analysis_method = analysis_method); if (this->dump_paraview) { std::stringstream base_name; base_name << "bar" << analysis_method << this->type; this->model->setBaseName(base_name.str()); this->model->addDumpFieldVector("displacement"); this->model->addDumpField("mass"); this->model->addDumpField("velocity"); this->model->addDumpField("acceleration"); this->model->addDumpFieldVector("external_force"); this->model->addDumpFieldVector("internal_force"); this->model->addDumpField("stress"); this->model->addDumpField("strain"); } time_step = this->model->getStableTimeStep() / 10.; this->model->setTimeStep(time_step); // std::cout << "timestep: " << time_step << std::endl; const auto & position = this->mesh->getNodes(); auto & displacement = this->model->getDisplacement(); auto & velocity = this->model->getVelocity(); constexpr auto dim = parent::spatial_dimension; Verification verif; for (auto && tuple : zip(make_view(position, dim), make_view(displacement, dim), make_view(velocity, dim))) { verif.displacement(std::get<1>(tuple), std::get<0>(tuple), 0.); verif.velocity(std::get<2>(tuple), std::get<0>(tuple), 0.); } if (dump_paraview) this->model->dump(); /// boundary conditions this->model->applyBC(SolutionFunctor(0., *this->model), "BC"); wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1 simulation_time = 5 / wave_velocity; max_steps = simulation_time / time_step; // 100 auto ndump = 50; dump_freq = max_steps / ndump; } protected: Real time_step; Real wave_velocity; Real simulation_time; UInt max_steps; UInt dump_freq; bool dump_paraview{false}; }; template using analysis_method_t = std::integral_constant; #ifdef AKANTU_IMPLICIT using TestAnalysisTypes = std::tuple, analysis_method_t<_explicit_lumped_mass>>; #else using TestAnalysisTypes = std::tuple>; #endif using TestTypes = gtest_list_t>; TYPED_TEST_CASE(TestSMMFixtureBar, TestTypes); /* -------------------------------------------------------------------------- */ TYPED_TEST(TestSMMFixtureBar, DynamicsExplicit) { constexpr auto dim = TestFixture::spatial_dimension; Real current_time = 0.; const auto & position = this->mesh->getNodes(); const auto & displacement = this->model->getDisplacement(); UInt nb_nodes = this->mesh->getNbNodes(); UInt nb_global_nodes = this->mesh->getNbGlobalNodes(); Real max_error{0.}; Array displacement_solution(nb_nodes, dim); Verification verif; for (UInt s = 0; s < this->max_steps; ++s, current_time += this->time_step) { if (s % this->dump_freq == 0 && this->dump_paraview) this->model->dump(); /// boundary conditions this->model->applyBC( SolutionFunctor(current_time, *this->model), "BC"); // compute the disp solution for (auto && tuple : zip(make_view(position, dim), make_view(displacement_solution, dim))) { verif.displacement(std::get<1>(tuple), std::get<0>(tuple), current_time); } // compute the error solution Real disp_error = 0.; for (auto && tuple : zip(make_view(displacement, dim), make_view(displacement_solution, dim))) { auto diff = std::get<1>(tuple) - std::get<0>(tuple); disp_error += diff.dot(diff); } this->mesh->getCommunicator().allReduce(disp_error, SynchronizerOperation::_sum); disp_error = sqrt(disp_error) / nb_global_nodes; max_error = std::max(disp_error, max_error); ASSERT_NEAR(disp_error, 0., 1e-3); this->model->solveStep(); } // std::cout << "max error: " << max_error << std::endl; } }