diff --git a/doc/dev-doc/manual/getting_started.rst b/doc/dev-doc/manual/getting_started.rst index ffd4b4bb3..4ba394ab2 100644 --- a/doc/dev-doc/manual/getting_started.rst +++ b/doc/dev-doc/manual/getting_started.rst @@ -1,302 +1,274 @@ Getting Started =============== Compiling Akantu ---------------- Akantu is a `CMake `_ project, so to configure it, you can either follow the usual way:: > cd akantu > mkdir build > cd build > ccmake .. [ Set the options that you need ] > make > make install All the Akantu options are documented in Appendix app:package-dependencies. Writing a ``main`` function --------------------------- Akantu first needs to be initialized. The memory management included in the core library handles the correct allocation and de-allocation of vectors, structures and/or objects. Moreover, in parallel computations, the initialization procedure performs the communication setup. This is achieved by the function :cpp:func:`initialize ` that is used as follows:: #include "aka_common.hh" #include "..." using namespace akantu; int main(int argc, char *argv[]) { initialize("input_file.dat", argc, argv); // your code ... } The :cpp:func:`initialize ` function takes the text inpute file and the program parameters which can be parsed by Akantu in due form (see sect:parser). Obviously it is necessary to include all files needed in main. In this manual all provided code implies the usage of ``akantu`` as namespace. Creating and Loading a Mesh --------------------------- In its current state, Akantu supports three types of meshes: Gmsh, Abaqus and Diana. Once a :cpp:class:`akantu::Mesh` object is created with a given spatial dimension, it can be filled by reading a mesh input file. The method :cpp:func:`read ` of the class :cpp:class:`Mesh ` infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified. :: UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); // Reading Gmsh files mesh.read("my_gmsh_mesh.msh"); mesh.read("my_gmsh_mesh", _miot_gmsh); The Gmsh reader adds the geometrical and physical tags as mesh data. The physical values are stored as a :cpp:type:`UInt ` data called ``tag_0``, if a string name is provided it is stored as a ``std::string`` data named ``physical_names``. The geometrical tag is stored as a :cpp:type:`UInt ` data named ``tag_1``. Using Arrays ------------ Data in Akantu can be stored in data containers implemented by the :cpp:class:`akantu::Array` class. In its most basic usage, the :cpp:class:`Array ` class implemented in \akantu is similar to the ``std::vector`` class of the Standard Template Library (STL) for C++. A simple :cpp:class:`Array ` containing a sequence of ``nb_element`` values (of a given type) can be generated with:: Array example_array(nb_element); where ``type`` usually is ``Real``, ``Int``, ``UInt`` or ``bool``. Each value is associated to an index, so that data can be accessed by typing:: auto & val = example_array(index); ``Arrays`` can also contain tuples of values for each index. In that case, the number of components per tuple must be specified at the :cpp:class:`Array ` creation. For example, if we want to create an :cpp:class:`Array ` to store the coordinates (sequences of three values) of ten nodes, the appropriate code is the following:: UInt nb_nodes = 10; UInt spatial_dimension = 3; Array position(nb_nodes, spatial_dimension); In this case the :math:`x` position of the eighth node number will be given by ``position(7, 0)`` (in C++, numbering starts at 0 and not 1). If the number of components for the sequences is not specified, the default value of 1 is used. Here is a list of some basic operations that can be performed on :cpp:class:`Array `: - ``resize(size)`` change the size of the :cpp:class:`Array `. - ``clear()`` set all entries of the :cpp:class:`Array ` to zero. - ``set(t)`` set all entries of the :cpp:class:`Array ` to ``t``. - ``copy(const Array & other)`` copy another :cpp:class:`Array ` into the current one. The two :cpp:class:`Array ` should have the same number of components. - ``push_back(tuple)`` append a tuple with the correct number of components at the end of the :cpp:class:`Array `. - ``erase(i)`` erase the value at the i-th position. - ``find(value)`` search ``value`` in the current :cpp:class:`Array `. Return position index of the first occurence or -1 if not found. - ``storage()`` Return the address of the allocated memory of the :cpp:class:`Array `. Array iterators ------------------- It is very common in Akantu to loop over arrays to perform a specific treatment. This ranges from geometric calculation on nodal quantities to tensor algebra (in constitutive laws for example). The :cpp:class:`Array ` object has the possibility to request iterators in order to make the writing of loops easier and enhance readability. For instance, a loop over the nodal coordinates can be performed like:: // accessing the nodal coordinates Array // with spatial_dimension components const auto & nodes = mesh.getNodes(); for (const auto & coords : make_view(nodes, spatial_dimension)) { // do what you need .... } In that example, each ``coords`` is a ``Vector`` containing geometrical array of size ``spatial_dimension`` and the iteration is conveniently performed by the :cpp:class:`Array ` iterator. The :cpp:class:`Array ` object is intensively used to store second order tensor values. In that case, it should be specified that the returned object type is a matrix when constructing the iterator. This is done when calling the :cpp:func:`make_view `. For instance, assuming that we have a :cpp:class:`Array ` storing stresses, we can loop over the stored tensors by:: for (const auto & stress : make_view(stresses, spatial_dimension, spatial_dimension)) { // stress is of type `const Matrix&` } In that last example, the :cpp:class:`Matrix ` objects are ``spatial_dimension`` :math:`\times` ``spatial_dimension`` matrices. The light objects :cpp:class:`Matrix ` and :cpp:class:`Vector ` can be used and combined to do most common linear algebra. If the number of component is 1, it is possible to use :cpp:func:`make_view ` to this effect. In general, a mesh consists of several kinds of elements. Consequently, the amount of data to be stored can differ for each element type. The straightforward example is the connectivity array, namely the sequences of nodes belonging to each element (linear triangular elements have fewer nodes than, say, rectangular quadratic elements etc.). A particular data structure called :cpp:class:`ElementTypeMapArray ` is provided to easily manage this kind of data. It consists of a group of ``Arrays``, each associated to an element type. The following code can retrieve the ``ElementTypeMapArray`` which stores the connectivity arrays for a mesh:: const ElementTypeMapArray & connectivities = mesh.getConnectivities(); Then, the specific array associated to a given element type can be obtained by:: const Array & connectivity_triangle = connectivities(_triangle_3); where the first order 3-node triangular element was used in the presented piece of code. Vector & Matrix ``````````````` The :cpp:class:`Array ` iterators as presented in the previous section can be shaped as :cpp:class:`Vector ` or :cpp:class:`Matrix `. This objects represent 1st and 2nd order tensors. As such they come with some functionalities that we will present a bit more into detail in this here. ``Vector`` ''''''''''''' - Accessors: - ``v(i)`` gives the ``i`` -th component of the vector ``v`` - ``v[i]`` gives the ``i`` -th component of the vector ``v`` - ``v.size()`` gives the number of component - Level 1: (results are scalars) - ``v.norm()`` returns the geometrical norm (:math:`L_2`) - - ``v.norm()`` returns the :math:`L_N` norm defined as :math:`\left(\sum_i - |v(i)|^N\right)^{1/N}`. N can take any positive integer value. - There are also some particular values for the most commonly used - norms, ``L_1`` for the Manhattan norm, ``L_2`` for the geometrical - norm and ``L_inf`` for the norm infinity. + - ``v.norm()`` returns the :math:`L_N` norm defined as :math:`\left(\sum_i|v(i)|^N\right)^{1/N}`. N can take any positive integer value. + There are also some particular values for the most commonly used norms, ``L_1`` for the Manhattan norm, ``L_2`` for the geometrical norm and ``L_inf`` for the norm infinity. - ``v.dot(x)`` return the dot product of ``v`` and ``x`` - ``v.distance(x)`` return the geometrical norm of :math:`v - x` - Level 2: (results are vectors) - ``v += s``, ``v -= s``, ``v *= s``, ``v /= s`` those are element-wise operators that sum, substract, multiply or divide all the component of ``v`` by the scalar ``s`` - ``v += x``, ``v -= x`` sums or substracts the vector ``x`` to/from ``v`` - - ``v.mul(A, x, alpha)`` stores the result of :math:`\alpha \boldsymbol{A} - \vec{x}` in ``v``, :math:`\alpha` is equal to 1 by default - - ``v.solve(A, b)`` stores the result of the resolution of the system - :math:`\boldsymbol{A} \vec{x} = \vec{b}` in ``v`` + - ``v.mul(A, x, alpha)`` stores the result of :math:`\alpha \boldsymbol{A} \vec{x}` in ``v``, :math:`\alpha` is equal to 1 by default + - ``v.solve(A, b)`` stores the result of the resolution of the system :math:`\boldsymbol{A} \vec{x} = \vec{b}` in ``v`` - ``v.crossProduct(v1, v2)`` computes the cross product of ``v1`` and ``v2`` and stores the result in ``v`` ``Matrix`` ''''''''''''' - Accessors: - ``A(i, j)`` gives the component :math:`A_{ij}` of the matrix ``A`` - ``A(i)`` gives the :math:`i^{th}` column of the matrix as a ``Vector`` - - ``A[k]`` gives the :math:`k^{th}` component of the matrix, matrices are - stored in a column major way, which means that to access :math:`A_{ij}`, - :math:`k = i + j M` + - ``A[k]`` gives the :math:`k^{th}` component of the matrix, matrices are stored in a column major way, which means that to access :math:`A_{ij}`, :math:`k = i + j M` - ``A.rows()`` gives the number of rows of ``A`` (:math:`M`) - ``A.cols()`` gives the number of columns of ``A`` (:math:`N`) - - ``A.size()`` gives the number of component in the matrix (:math:`M \times - N`) + - ``A.size()`` gives the number of component in the matrix (:math:`M \times N`) - Level 1: (results are scalars) - ``A.norm()`` is equivalent to ``A.norm()`` - ``A.norm()`` returns the :math:`L_N` norm defined as :math:`\left(\sum_i\sum_j |A(i,j)|^N\right)^{1/N}`. N can take any positive integer value. There are also some particular values for the most commonly used norms, ``L_1`` for the Manhattan norm, ``L_2`` for the geometrical norm and ``L_inf`` for the norm infinity. - ``A.trace()`` return the trace of ``A`` - ``A.det()`` return the determinant of ``A`` - ``A.doubleDot(B)`` return the double dot product of ``A`` and ``B``, :math:`\mat{A}:\mat{B}` - Level 3: (results are matrices) - ``A.eye(s)``, ``Matrix::eye(s)`` fills/creates a matrix with the :math:`s\mat{I}` with :math:`\mat{I}` the identity matrix - ``A.inverse(B)`` stores :math:`\mat{B}^{-1}` in ``A`` - ``A.transpose()`` returns :math:`\mat{A}^{t}` - ``A.outerProduct(v1, v2)`` stores :math:`\vec{v_1} \vec{v_2}^{t}` in ``A`` - ``C.mul(A, B, alpha)``: stores the result of the product of ``A`` and code{B} time the scalar ``alpha`` in ``C``. ``t_A`` and ``t_B`` are boolean defining if ``A`` and ``B`` should be transposed or not. - +----------+----------+--------------+ - |``t_A`` |``t_B`` |result | - | | | | - +----------+----------+--------------+ - |false |false |:math:`\mat{C}| - | | |= \alpha | - | | |\mat{A} | - | | |\mat{B}` | - | | | | - +----------+----------+--------------+ - |false |true |:math:`\mat{C}| - | | |= \alpha | - | | |\mat{A} | - | | |\mat{B}^t` | - | | | | - +----------+----------+--------------+ - |true |false |:math:`\mat{C}| - | | |= \alpha | - | | |\mat{A}^t | - | | |\mat{B}` | - | | | | - +----------+----------+--------------+ - |true |true |:math:`\mat{C}| - | | |= \alpha | - | | |\mat{A}^t | - | | |\mat{B}^t` | - +----------+----------+--------------+ - - - ``A.eigs(d, V)`` this method computes the eigenvalues and eigenvectors of - ``A`` and store the results in ``d`` and ``V`` such that :math:`d(i) = - \lambda_i` and :math:`V(i) = \vec{v_i}` with :math:`\mat{A}\vec{v_i} = - \lambda_i\vec{v_i}` and :math:`\lambda_1 > ... > \lambda_i > ... > - \lambda_N` + +----------+----------+---------------------------------------------+ + |``t_A`` |``t_B`` |result | + +----------+----------+---------------------------------------------+ + |false |false |:math:`\mat{C} = \alpha \mat{A} \mat{B}` | + +----------+----------+---------------------------------------------+ + |false |true |:math:`\mat{C} = \alpha \mat{A} \mat{B}^t` | + +----------+----------+---------------------------------------------+ + |true |false |:math:`\mat{C} = \alpha \mat{A}^t \mat{B}` | + +----------+----------+---------------------------------------------+ + |true |true |:math:`\mat{C} = \alpha \mat{A}^t \mat{B}^t` | + +----------+----------+---------------------------------------------+ + + - ``A.eigs(d, V)`` this method computes the eigenvalues and eigenvectors of ``A`` and store the results in ``d`` and ``V`` such that :math:`d(i) = \lambda_i` and :math:`V(i) = \vec{v_i}` with :math:`\mat{A}\vec{v_i} = \lambda_i\vec{v_i}` and :math:`\lambda_1 > ... > \lambda_i > ... > \lambda_N` diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index be17e197d..ab91f4436 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,642 +1,645 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Constants */ /* -------------------------------------------------------------------------- */ namespace { [[gnu::unused]] constexpr UInt _all_dimensions{ std::numeric_limits::max()}; #ifdef AKANTU_NDEBUG [[gnu::unused]] constexpr Real REAL_INIT_VALUE{0.}; #else [[gnu::unused]] constexpr Real REAL_INIT_VALUE{ std::numeric_limits::quiet_NaN()}; #endif } // namespace /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; using MemoryID = UInt; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "aka_enum_macros.hh" /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of /// events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) // clang-format on /// enum ModelType defines which type of physics is solved -AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) -AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) -AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) +AKANTU_CLASS_ENUM_DECLARE_WITH_IO(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { model, solid_mechanics_model, solid_mechanics_model_cohesive, heat_transfer_model, structural_mechanics_model, embedded_model, }; #endif /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_NON_LINEAR_SOLVER_TYPES \ (linear) \ (newton_raphson) \ (newton_raphson_modified) \ (lumped) \ (gmres) \ (bfgs) \ (cg) \ (auto) // clang-format on -AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) -AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType, - AKANTU_NON_LINEAR_SOLVER_TYPES) -AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType, - AKANTU_NON_LINEAR_SOLVER_TYPES) +AKANTU_CLASS_ENUM_DECLARE_WITH_IO(NonLinearSolverType, + AKANTU_NON_LINEAR_SOLVER_TYPES) #else /// Type of non linear resolution available in akantu enum class NonLinearSolverType { _linear, ///< No non linear convergence loop _newton_raphson, ///< Regular Newton-Raphson _newton_raphson_modified, ///< Newton-Raphson with initial tangent _lumped, ///< Case of lumped mass or equivalent matrix - _gmres, + _gmres, ///< _bfgs, _cg, _auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_TIME_STEP_SOLVER_TYPE \ (static) \ (dynamic) \ (dynamic_lumped) \ (not_defined) // clang-format on -AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) -AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType, - AKANTU_TIME_STEP_SOLVER_TYPE) -AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) +AKANTU_CLASS_ENUM_DECLARE_WITH_IO(TimeStepSolverType, + AKANTU_TIME_STEP_SOLVER_TYPE) #else /// Type of time stepping solver enum class TimeStepSolverType { _static, ///< Static solution _dynamic, ///< Dynamic solver _dynamic_lumped, ///< Dynamic solver with lumped mass _not_defined, ///< For not defined cases }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_INTEGRATION_SCHEME_TYPE \ (pseudo_time) \ (forward_euler) \ (trapezoidal_rule_1) \ (backward_euler) \ (central_difference) \ (fox_goodwin) \ (trapezoidal_rule_2) \ (linear_acceleration) \ (newmark_beta) \ (generalized_trapezoidal) // clang-format on -AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) -AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType, - AKANTU_INTEGRATION_SCHEME_TYPE) -AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType, - AKANTU_INTEGRATION_SCHEME_TYPE) +AKANTU_CLASS_ENUM_DECLARE_WITH_IO(IntegrationSchemeType, + AKANTU_INTEGRATION_SCHEME_TYPE) #else /// Type of integration scheme enum class IntegrationSchemeType { _pseudo_time, ///< Pseudo Time _forward_euler, ///< GeneralizedTrapezoidal(0) _trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _backward_euler, ///< GeneralizedTrapezoidal(1) _central_difference, ///< NewmarkBeta(0, 1/2) _fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ (residual) \ (solution) \ (residual_mass_wgh) // clang-format on -AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria, - AKANTU_SOLVE_CONVERGENCE_CRITERIA) -AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria, - AKANTU_SOLVE_CONVERGENCE_CRITERIA) -AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria, - AKANTU_SOLVE_CONVERGENCE_CRITERIA) +AKANTU_CLASS_ENUM_DECLARE_WITH_IO(SolveConvergenceCriteria, + AKANTU_SOLVE_CONVERGENCE_CRITERIA) #else /// enum SolveConvergenceCriteria different convergence criteria enum class SolveConvergenceCriteria { _residual, ///< Use residual to test the convergence _solution, ///< Use solution to test the convergence _residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; #endif /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum MatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; +#if !defined(DOXYGEN) +// clang-format off +#define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ + (residual) \ + (solution) \ + (residual_mass_wgh) +#define AKANTU_CLUSTERING_STRATEGY \ + (facets) \ + (nodes) +// clang-format on +AKANTU_CLASS_ENUM_DECLARE_WITH_IO(ClusteringStrategy, + AKANTU_CLUSTERING_STRATEGY) +#else +/// @enum ClusteringStrategy type of strategy to use to detect +/// connected elements while clustering a mesh +enum class ClusteringStrategy { _facets, _nodes }; +#endif + /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_SYNCHRONIZATION_TAG \ (whatever) \ (update) \ (ask_nodes) \ (size) \ (smm_mass) \ (smm_for_gradu) \ (smm_boundary) \ (smm_uv) \ (smm_res) \ (smm_init_mat) \ (smm_stress) \ (smmc_facets) \ (smmc_facets_conn) \ (smmc_facets_stress) \ (smmc_damage) \ (giu_global_conn) \ (ce_groups) \ (gm_clusters) \ (htm_temperature) \ (htm_gradient_temperature) \ (htm_phi) \ (htm_gradient_phi) \ (mnl_for_average) \ (mnl_weight) \ (nh_criterion) \ (test) \ (user_1) \ (user_2) \ (material_id) \ (for_dump) \ (cf_nodal) \ (cf_incr) \ (solver_solution) // clang-format on -AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) -AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) +AKANTU_CLASS_ENUM_DECLARE_WITH_IO(SynchronizationTag, + AKANTU_SYNCHRONIZATION_TAG) #else /// @enum SynchronizationTag type of synchronizations enum class SynchronizationTag { //--- Generic tags --- _whatever, _update, _ask_nodes, _size, //--- SolidMechanicsModel tags --- _smm_mass, ///< synchronization of the SolidMechanicsModel.mass _smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _smm_uv, ///< synchronization of the nodal velocities and displacement _smm_res, ///< synchronization of the nodal residual _smm_init_mat, ///< synchronization of the data to initialize materials _smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _smmc_facets, ///< synchronization of facet data to setup facet synch _smmc_facets_conn, ///< synchronization of facet global connectivity _smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _htm_temperature, ///< synchronization of the nodal temperature _htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _htm_phi, ///< synchronization of the nodal level set value phi _htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _mnl_for_average, ///< synchronization of data to average in non local /// material _mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _nh_criterion, // --- General tags --- _test, ///< Test tag _user_1, ///< tag for user simulations _user_2, ///< tag for user simulations _material_id, ///< synchronization of the material ids _for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _cf_nodal, ///< synchronization of disp, velo, and current position _cf_incr, ///< synchronization of increment // --- Solver tags --- _solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; #endif /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } std::ostream & operator<<(std::ostream & stream, NodeFlag flag); } // namespace akantu AKANTU_ENUM_HASH(GhostType) namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT ' ' #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_DEREF_PTR(name, ptr) \ inline decltype(auto) get##name() const { \ if (not ptr) { \ AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ } \ return (*ptr); \ } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ struct TensorTrait {}; struct TensorProxyTrait {}; } // namespace akantu /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ namespace aka { /* ------------------------------------------------------------------------ */ template using is_tensor = std::is_base_of; template using is_tensor_proxy = std::is_base_of; /* ------------------------------------------------------------------------ */ template using is_scalar = std::is_arithmetic; /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> bool is_of_type(T && t) { return ( dynamic_cast>::value, std::add_const_t, R>>>(&t) != nullptr); } /* -------------------------------------------------------------------------- */ template bool is_of_type(std::unique_ptr & t) { return ( dynamic_cast::value, std::add_const_t, R>>>( t.get()) != nullptr); } /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { static_assert( disjunction< std::is_base_of, std::decay_t>, // down-cast std::is_base_of, std::decay_t> // up-cast >::value, "Type T and R are not valid for a as_type conversion"); return dynamic_cast>::value, std::add_const_t, R>>>(t); } /* -------------------------------------------------------------------------- */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { return &as_type(*t); } /* -------------------------------------------------------------------------- */ template decltype(auto) as_type(const std::shared_ptr & t) { return std::dynamic_pointer_cast(t); } } // namespace aka #include "aka_common_inline_impl.hh" #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(); #define AKANTU_CURRENT_FUNCTION \ (std::string(__func__) + "():" + std::to_string(__LINE__)) } // namespace akantu /* -------------------------------------------------------------------------- */ #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/common/aka_enum_macros.hh b/src/common/aka_enum_macros.hh index e0e82be7e..4f9eb19f6 100644 --- a/src/common/aka_enum_macros.hh +++ b/src/common/aka_enum_macros.hh @@ -1,132 +1,137 @@ /** * @file aka_enum_macros.hh * * @author Nicolas Richart * * @date creation Wed Oct 31 2018 * * @brief A Documented file. * * * 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 #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_ENUM_MACROS_HH__ #define __AKANTU_AKA_ENUM_MACROS_HH__ #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))) #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 #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_CLASS_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, prefix) \ } \ AKANTU_ENUM_HASH(type_name) \ namespace std { \ inline string to_string(const ::akantu::type_name & type) { \ using namespace akantu; \ static unordered_map<::akantu::type_name, 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, prefix, list))}; \ return convert.at(type); \ } \ } \ namespace akantu { \ inline std::ostream & operator<<(std::ostream & stream, \ const type_name & type) { \ stream << std::to_string(type); \ return stream; \ } #define AKANTU_ENUM_INPUT_STREAM_(type_name, list, prefix) \ 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, prefix, list))}; \ try { \ type = convert.at(str); \ } catch (std::out_of_range &) { \ std::ostringstream values; \ std::for_each(convert.begin(), convert.end(), [&values](auto && pair) { \ static bool first = true; \ if (not first) \ values << ", "; \ values << "\"" << pair.first << "\""; \ first = false; \ }); \ AKANTU_EXCEPTION("The value " << str << " is not a valid " \ << BOOST_PP_STRINGIZE(type_name) \ << " valid values are " << values.str()); \ } \ return stream; \ } #define AKANTU_CLASS_ENUM_OUTPUT_STREAM(type_name, list) \ AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, type_name::_) #define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ AKANTU_ENUM_OUTPUT_STREAM_(type_name, list, ) #define AKANTU_CLASS_ENUM_INPUT_STREAM(type_name, list) \ AKANTU_ENUM_INPUT_STREAM_(type_name, list, type_name::_) #define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ AKANTU_ENUM_INPUT_STREAM_(type_name, list, ) +#define AKANTU_CLASS_ENUM_DECLARE_WITH_IO(type_name, list) \ + AKANTU_CLASS_ENUM_DECLARE(type_name, list) \ + AKANTU_CLASS_ENUM_OUTPUT_STREAM(type_name, list) \ + AKANTU_CLASS_ENUM_INPUT_STREAM(type_name, list) + #endif /* __AKANTU_AKA_ENUM_MACROS_HH__ */ diff --git a/src/mesh/group_manager.cc b/src/mesh/group_manager.cc index 6884b902e..d1c941ec8 100644 --- a/src/mesh/group_manager.cc +++ b/src/mesh/group_manager.cc @@ -1,1038 +1,1110 @@ /** * @file group_manager.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Feb 20 2018 * * @brief Stores information about ElementGroup and NodeGroup * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "group_manager.hh" #include "aka_csr.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -GroupManager::GroupManager(const Mesh & mesh, const ID & id, +GroupManager::GroupManager(Mesh & mesh, const ID & id, const MemoryID & mem_id) : id(id), memory_id(mem_id), mesh(mesh) { AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ GroupManager::~GroupManager() = default; /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::createNodeGroup(const std::string & group_name, bool replace_group) { AKANTU_DEBUG_IN(); auto it = node_groups.find(group_name); if (it != node_groups.end()) { if (replace_group) { it->second.reset(); } else { AKANTU_EXCEPTION( "Trying to create a node group that already exists:" << group_name); } } std::stringstream sstr; sstr << this->id << ":" << group_name << "_node_group"; auto && ptr = std::make_unique(group_name, mesh, sstr.str(), memory_id); auto & node_group = *ptr; // \todo insert_or_assign in c++17 if (it != node_groups.end()) { it->second = std::move(ptr); } else { node_groups[group_name] = std::move(ptr); } AKANTU_DEBUG_OUT(); return node_group; } /* -------------------------------------------------------------------------- */ template NodeGroup & GroupManager::createFilteredNodeGroup(const std::string & group_name, const NodeGroup & source_node_group, T & filter) { AKANTU_DEBUG_IN(); NodeGroup & node_group = this->createNodeGroup(group_name); node_group.append(source_node_group); if (T::type == FilterFunctor::_node_filter_functor) { node_group.applyNodeFilter(filter); } else { AKANTU_ERROR("ElementFilter cannot be applied to NodeGroup yet." << " Needs to be implemented."); } AKANTU_DEBUG_OUT(); return node_group; } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, bool replace_group) { AKANTU_DEBUG_IN(); auto it = element_groups.find(group_name); if (it != element_groups.end()) { if (replace_group) { it->second.reset(); } else { AKANTU_EXCEPTION("Trying to create a element group that already exists:" << group_name); } } NodeGroup & new_node_group = createNodeGroup(group_name + "_nodes", replace_group); auto && ptr = std::make_unique( group_name, mesh, new_node_group, dimension, this->id + ":" + group_name + "_element_group", memory_id); auto & element_group = *ptr; if (it != element_groups.end()) { it->second = std::move(ptr); } else { element_groups[group_name] = std::move(ptr); } AKANTU_DEBUG_OUT(); return element_group; } /* -------------------------------------------------------------------------- */ void GroupManager::destroyElementGroup(const std::string & group_name, bool destroy_node_group) { AKANTU_DEBUG_IN(); auto eit = element_groups.find(group_name); if (eit != element_groups.end()) { if (destroy_node_group) destroyNodeGroup(eit->second->getNodeGroup().getName()); element_groups.erase(eit); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::destroyNodeGroup(const std::string & group_name) { AKANTU_DEBUG_IN(); auto nit = node_groups.find(group_name); if (nit != node_groups.end()) { node_groups.erase(nit); } AKANTU_DEBUG_OUT(); } // /* -------------------------------------------------------------------------- // */ void GroupManager::destroyAllElementGroups(bool destroy_node_groups) { // AKANTU_DEBUG_IN(); // if (destroy_node_groups) // for (auto && data : element_groups) { // destroyNodeGroup(std::get<1>(data)->getNodeGroup().getName()); // } // element_groups.clear(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group) { AKANTU_DEBUG_IN(); if (element_groups.find(group_name) != element_groups.end()) AKANTU_EXCEPTION( "Trying to create a element group that already exists:" << group_name); auto && ptr = std::make_unique( group_name, mesh, node_group, dimension, id + ":" + group_name + "_element_group", memory_id); auto & element_group = *ptr; element_groups[group_name] = std::move(ptr); AKANTU_DEBUG_OUT(); return element_group; } /* -------------------------------------------------------------------------- */ template ElementGroup & GroupManager::createFilteredElementGroup( const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter) { AKANTU_DEBUG_IN(); if (T::type == FilterFunctor::_node_filter_functor) { auto & filtered_node_group = this->createFilteredNodeGroup( group_name + "_nodes", node_group, filter); AKANTU_DEBUG_OUT(); return this->createElementGroup(group_name, dimension, filtered_node_group); } else if (T::type == FilterFunctor::_element_filter_functor) { AKANTU_ERROR( "Cannot handle an ElementFilter yet. Needs to be implemented."); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ class ClusterSynchronizer : public DataAccessor { using DistantIDs = std::set>; public: ClusterSynchronizer(GroupManager & group_manager, UInt element_dimension, std::string cluster_name_prefix, ElementTypeMapArray & element_to_fragment, const ElementSynchronizer & element_synchronizer, UInt nb_cluster) : group_manager(group_manager), element_dimension(element_dimension), cluster_name_prefix(std::move(cluster_name_prefix)), element_to_fragment(element_to_fragment), element_synchronizer(element_synchronizer), nb_cluster(nb_cluster) {} UInt synchronize() { Communicator & comm = Communicator::getStaticCommunicator(); UInt rank = comm.whoAmI(); UInt nb_proc = comm.getNbProc(); /// find starting index to renumber local clusters Array nb_cluster_per_proc(nb_proc); nb_cluster_per_proc(rank) = nb_cluster; comm.allGather(nb_cluster_per_proc); starting_index = std::accumulate(nb_cluster_per_proc.begin(), nb_cluster_per_proc.begin() + rank, 0); UInt global_nb_fragment = std::accumulate(nb_cluster_per_proc.begin() + rank, nb_cluster_per_proc.end(), starting_index); /// create the local to distant cluster pairs with neighbors element_synchronizer.synchronizeOnce(*this, SynchronizationTag::_gm_clusters); /// count total number of pairs Array nb_pairs(nb_proc); // This is potentially a bug for more than // 2**31 pairs, but due to a all gatherv after // it must be int to match MPI interfaces nb_pairs(rank) = distant_ids.size(); comm.allGather(nb_pairs); UInt total_nb_pairs = std::accumulate(nb_pairs.begin(), nb_pairs.end(), 0); /// generate pairs global array UInt local_pair_index = std::accumulate(nb_pairs.storage(), nb_pairs.storage() + rank, 0); Array total_pairs(total_nb_pairs, 2); for (auto & ids : distant_ids) { total_pairs(local_pair_index, 0) = ids.first; total_pairs(local_pair_index, 1) = ids.second; ++local_pair_index; } /// communicate pairs to all processors nb_pairs *= 2; comm.allGatherV(total_pairs, nb_pairs); /// renumber clusters /// generate fragment list std::vector> global_clusters; UInt total_nb_cluster = 0; Array is_fragment_in_cluster(global_nb_fragment, 1, false); std::queue fragment_check_list; while (total_pairs.size() != 0) { /// create a new cluster ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; UInt first_fragment = total_pairs(0, 0); UInt second_fragment = total_pairs(0, 1); total_pairs.erase(0); fragment_check_list.push(first_fragment); fragment_check_list.push(second_fragment); while (!fragment_check_list.empty()) { UInt current_fragment = fragment_check_list.front(); UInt * total_pairs_end = total_pairs.storage() + total_pairs.size() * 2; UInt * fragment_found = std::find(total_pairs.storage(), total_pairs_end, current_fragment); if (fragment_found != total_pairs_end) { UInt position = fragment_found - total_pairs.storage(); UInt pair = position / 2; UInt other_index = (position + 1) % 2; fragment_check_list.push(total_pairs(pair, other_index)); total_pairs.erase(pair); } else { fragment_check_list.pop(); current_cluster.insert(current_fragment); is_fragment_in_cluster(current_fragment) = true; } } } /// add to FragmentToCluster all local fragments for (UInt c = 0; c < global_nb_fragment; ++c) { if (!is_fragment_in_cluster(c)) { ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; current_cluster.insert(c); } } /// reorganize element groups to match global clusters for (UInt c = 0; c < global_clusters.size(); ++c) { /// create new element group corresponding to current cluster std::stringstream sstr; sstr << cluster_name_prefix << "_" << c; ElementGroup & cluster = group_manager.createElementGroup(sstr.str(), element_dimension, true); auto it = global_clusters[c].begin(); auto end = global_clusters[c].end(); /// append to current element group all fragments that belong to /// the same cluster if they exist for (; it != end; ++it) { Int local_index = *it - starting_index; if (local_index < 0 || local_index >= Int(nb_cluster)) continue; std::stringstream tmp_sstr; tmp_sstr << "tmp_" << cluster_name_prefix << "_" << local_index; AKANTU_DEBUG_ASSERT(group_manager.elementGroupExists(tmp_sstr.str()), "Temporary fragment \"" << tmp_sstr.str() << "\" not found"); cluster.append(group_manager.getElementGroup(tmp_sstr.str())); group_manager.destroyElementGroup(tmp_sstr.str(), true); } } return total_nb_cluster; } private: /// functions for parallel communications inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_gm_clusters) return elements.size() * sizeof(UInt); return 0; } inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override { if (tag != SynchronizationTag::_gm_clusters) return; Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { const Element & el = *el_it; /// for each element pack its global cluster index buffer << element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; } } inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override { if (tag != SynchronizationTag::_gm_clusters) return; Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { UInt distant_cluster; buffer >> distant_cluster; const Element & el = *el_it; UInt local_cluster = element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; distant_ids.insert(std::make_pair(local_cluster, distant_cluster)); } } private: GroupManager & group_manager; UInt element_dimension; std::string cluster_name_prefix; ElementTypeMapArray & element_to_fragment; const ElementSynchronizer & element_synchronizer; UInt nb_cluster; DistantIDs distant_ids; UInt starting_index; }; /* -------------------------------------------------------------------------- */ /// \todo this function doesn't work in 1D UInt GroupManager::createBoundaryGroupFromGeometry() { UInt spatial_dimension = mesh.getSpatialDimension(); return createClusters(spatial_dimension - 1, "boundary"); } +/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- */ +class ClusteringStrategyImplementation { +public: + ClusteringStrategyImplementation(Mesh & mesh, + ElementTypeMapArray & seen_elements) + : mesh(mesh), seen_elements(seen_elements) {} + + virtual ~ClusteringStrategyImplementation() = default; + + virtual void queueConnectedElements(const Element & element, + std::queue & queue) = 0; + +protected: + void queueElement(const Element & element, std::queue & queue) { + auto & seen = seen_elements(element); + if (not seen) { + seen = true; + queue.push(element); + } + } + +protected: + Mesh & mesh; + +private: + ElementTypeMapArray & seen_elements; +}; + +/* -------------------------------------------------------------------------- */ +class ClusteringStrategyFacets : public ClusteringStrategyImplementation { +public: + ClusteringStrategyFacets(Mesh & mesh, + ElementTypeMapArray & seen_elements, + UInt element_dimension, Mesh * mesh_facets = nullptr) + : ClusteringStrategyImplementation(mesh, seen_elements), + element_dimension(element_dimension), mesh_facets(mesh_facets) { + + if (not this->mesh_facets and this->element_dimension > 0) { + MeshAccessor mesh_accessor(this->mesh); + this->mesh_facets = std::make_unique( + this->mesh.getSpatialDimension(), mesh_accessor.getNodesSharedPtr(), + "mesh_facets_for_clusters"); + + this->mesh_facets->defineMeshParent(mesh); + + MeshUtils::buildAllFacets(this->mesh, *this->mesh_facets, + this->element_dimension, + this->element_dimension - 1); + } + } + +public: + void queueConnectedElements(const Element & element, + std::queue & queue) override { + const Array & element_to_facet = + this->mesh_facets->getSubelementToElement(element.type, + element.ghost_type); + + UInt nb_facet_per_element = element_to_facet.getNbComponent(); + + for (UInt f = 0; f < nb_facet_per_element; ++f) { + const Element & facet = element_to_facet(element.element, f); + + if (facet == ElementNull) + continue; + + const std::vector & connected_elements = + this->mesh_facets->getElementToSubelement( + facet.type, facet.ghost_type)(facet.element); + + for (UInt elem = 0; elem < connected_elements.size(); ++elem) { + const Element & check_el = connected_elements[elem]; + + // check if this element has to be skipped + if (check_el == ElementNull || check_el == element) + continue; + + this->queueElement(check_el, queue); + } + } + } + +private: + UInt element_dimension; + std::unique_ptr mesh_facets; +}; + +/* -------------------------------------------------------------------------- */ +class ClusteringStrategyNodes : public ClusteringStrategyImplementation { +public: + ClusteringStrategyNodes(Mesh & mesh, + ElementTypeMapArray & seen_elements, + UInt element_dimension) + : ClusteringStrategyImplementation(mesh, seen_elements) { + MeshUtils::buildNode2Elements(this->mesh, node_to_elements, + element_dimension); + } + +public: + void queueConnectedElements(const Element & element, + std::queue & queue) override { + auto conn = const_cast(this->mesh).getConnectivity(element); + for (auto node : conn) { + for (auto element : node_to_elements.getRow(node)) { + this->queueElement(element, queue); + } + } + } + +protected: + CSR node_to_elements; +}; + /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix, + const ClusteringStrategy & clustering_strategy_type, const GroupManager::ClusteringFilter & filter) { return createClusters(element_dimension, std::move(cluster_name_prefix), - filter, mesh_facets); + clustering_strategy_type, filter, mesh_facets); } /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, std::string cluster_name_prefix, + const ClusteringStrategy & clustering_strategy_type, const GroupManager::ClusteringFilter & filter) { std::unique_ptr mesh_facets; if (!mesh_facets && element_dimension > 0) { MeshAccessor mesh_accessor(const_cast(mesh)); mesh_facets = std::make_unique(mesh.getSpatialDimension(), mesh_accessor.getNodesSharedPtr(), "mesh_facets_for_clusters"); mesh_facets->defineMeshParent(mesh); MeshUtils::buildAllFacets(mesh, *mesh_facets, element_dimension, element_dimension - 1); } return createClusters(element_dimension, std::move(cluster_name_prefix), - filter, *mesh_facets); + clustering_strategy_type, filter, *mesh_facets); } /* -------------------------------------------------------------------------- */ //// \todo if needed element list construction can be optimized by //// templating the filter class -UInt GroupManager::createClusters(UInt element_dimension, - const std::string & cluster_name_prefix, - const GroupManager::ClusteringFilter & filter, - Mesh & mesh_facets) { +UInt GroupManager::createClusters( + UInt element_dimension, const std::string & cluster_name_prefix, + const ClusteringStrategy & clustering_strategy_type, + const GroupManager::ClusteringFilter & filter, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); - UInt nb_proc = mesh.getCommunicator().getNbProc(); - std::string tmp_cluster_name_prefix = cluster_name_prefix; + auto nb_proc = mesh.getCommunicator().getNbProc(); + auto tmp_cluster_name_prefix = cluster_name_prefix; - ElementTypeMapArray * element_to_fragment = nullptr; + std::unique_ptr> element_to_fragment = nullptr; if (nb_proc > 1 && mesh.isDistributed()) { - element_to_fragment = - new ElementTypeMapArray("element_to_fragment", id, memory_id); - + element_to_fragment = std::make_unique>( + "element_to_fragment", id, memory_id); element_to_fragment->initialize( mesh, _nb_component = 1, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); - // mesh.initElementTypeMapArray(*element_to_fragment, 1, element_dimension, - // false, _ek_not_defined, true); tmp_cluster_name_prefix = "tmp_" + tmp_cluster_name_prefix; } ElementTypeMapArray seen_elements("seen_elements", id, memory_id); seen_elements.initialize(mesh, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); - // mesh.initElementTypeMapArray(seen_elements, 1, element_dimension, false, - // _ek_not_defined, true); for (auto ghost_type : ghost_types) { Element el; el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { el.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & seen_elements_array = seen_elements(type, ghost_type); for (UInt e = 0; e < nb_element; ++e) { el.element = e; if (!filter(el)) seen_elements_array(e) = true; } } } - Array checked_node(mesh.getNbNodes(), 1, false); + std::unique_ptr clustering_strategy; + switch (clustering_strategy_type) { + case ClusteringStrategy::_facets: + clustering_strategy = std::make_unique( + mesh, seen_elements, element_dimension, &mesh_facets); + break; + case ClusteringStrategy::_nodes: + clustering_strategy = std::make_unique( + mesh, seen_elements, element_dimension); + + break; + default: + AKANTU_EXCEPTION(clustering_strategy_type + << " is not a recognized clustering strategy"); + } UInt nb_cluster = 0; for (auto ghost_type : ghost_types) { Element uns_el; uns_el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { uns_el.type = type; Array & seen_elements_vec = seen_elements(uns_el.type, uns_el.ghost_type); for (UInt e = 0; e < seen_elements_vec.size(); ++e) { // skip elements that have been already seen if (seen_elements_vec(e) == true) continue; // set current element uns_el.element = e; seen_elements_vec(e) = true; /// create a new cluster std::stringstream sstr; sstr << tmp_cluster_name_prefix << "_" << nb_cluster; ElementGroup & cluster = createElementGroup(sstr.str(), element_dimension, true); ++nb_cluster; // point element are cluster by themself if (element_dimension == 0) { - cluster.add(uns_el); - - UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(uns_el.type); - Vector connect = - mesh.getConnectivity(uns_el.type, uns_el.ghost_type) - .begin(nb_nodes_per_element)[uns_el.element]; - for (UInt n = 0; n < nb_nodes_per_element; ++n) { - /// add element's nodes to the cluster - UInt node = connect[n]; - if (!checked_node(node)) { - cluster.addNode(node); - checked_node(node) = true; - } - } - + cluster.add(uns_el, true, false); continue; } std::queue element_to_add; element_to_add.push(uns_el); /// keep looping until current cluster is complete (no more /// connected elements) while (!element_to_add.empty()) { /// take first element and erase it in the queue Element el = element_to_add.front(); element_to_add.pop(); /// if parallel, store cluster index per element if (nb_proc > 1 && mesh.isDistributed()) (*element_to_fragment)(el.type, el.ghost_type)(el.element) = nb_cluster - 1; /// add current element to the cluster cluster.add(el); - const Array & element_to_facet = - mesh_facets.getSubelementToElement(el.type, el.ghost_type); - - UInt nb_facet_per_element = element_to_facet.getNbComponent(); - - for (UInt f = 0; f < nb_facet_per_element; ++f) { - const Element & facet = element_to_facet(el.element, f); - - if (facet == ElementNull) - continue; - - const std::vector & connected_elements = - mesh_facets.getElementToSubelement( - facet.type, facet.ghost_type)(facet.element); - - for (UInt elem = 0; elem < connected_elements.size(); ++elem) { - const Element & check_el = connected_elements[elem]; - - // check if this element has to be skipped - if (check_el == ElementNull || check_el == el) - continue; - - Array & seen_elements_vec_current = - seen_elements(check_el.type, check_el.ghost_type); - - if (seen_elements_vec_current(check_el.element) == false) { - seen_elements_vec_current(check_el.element) = true; - element_to_add.push(check_el); - } - } - } - - UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); - Vector connect = mesh.getConnectivity(el.type, el.ghost_type) - .begin(nb_nodes_per_element)[el.element]; - for (UInt n = 0; n < nb_nodes_per_element; ++n) { - /// add element's nodes to the cluster - UInt node = connect[n]; - if (!checked_node(node)) { - cluster.addNode(node, false); - checked_node(node) = true; - } - } + // connected component based on strategy + clustering_strategy->queueConnectedElements(el, element_to_add); } + + cluster.getNodeGroup().optimize(); } } } if (nb_proc > 1 && mesh.isDistributed()) { ClusterSynchronizer cluster_synchronizer( *this, element_dimension, cluster_name_prefix, *element_to_fragment, this->mesh.getElementSynchronizer(), nb_cluster); nb_cluster = cluster_synchronizer.synchronize(); - delete element_to_fragment; } if (mesh.isDistributed()) this->synchronizeGroupNames(); AKANTU_DEBUG_OUT(); return nb_cluster; } /* -------------------------------------------------------------------------- */ template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name) { std::set group_names; const auto & datas = mesh.getData(dataset_name); std::map group_dim; for (auto ghost_type : ghost_types) { for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements (" << type << ":" << ghost_type << ") in the map from MeshData (" << dataset.size() << ") " << dataset_name << " and in the mesh (" << nb_element << ")!"); for (UInt e(0); e < nb_element; ++e) { std::stringstream sstr; sstr << dataset(e); std::string gname = sstr.str(); group_names.insert(gname); auto it = group_dim.find(gname); if (it == group_dim.end()) { group_dim[gname] = mesh.getSpatialDimension(type); } else { it->second = std::max(it->second, mesh.getSpatialDimension(type)); } } } } auto git = group_names.begin(); auto gend = group_names.end(); for (; git != gend; ++git) createElementGroup(*git, group_dim[*git]); if (mesh.isDistributed()) this->synchronizeGroupNames(); Element el; for (auto ghost_type : ghost_types) { el.ghost_type = ghost_type; for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { el.type = type; const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements in the map from " "MeshData and in the mesh!"); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(el.type); Array::const_iterator> cit = mesh.getConnectivity(type, ghost_type).begin(nb_nodes_per_element); for (UInt e(0); e < nb_element; ++e, ++cit) { el.element = e; std::stringstream sstr; sstr << dataset(e); ElementGroup & group = getElementGroup(sstr.str()); group.add(el, false, false); const Vector & connect = *cit; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; group.addNode(node, false); } } } } git = group_names.begin(); for (; git != gend; ++git) { getElementGroup(*git).optimize(); } } template void GroupManager::createGroupsFromMeshData( const std::string & dataset_name); template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name); /* -------------------------------------------------------------------------- */ void GroupManager::createElementGroupFromNodeGroup( const std::string & name, const std::string & node_group_name, UInt dimension) { NodeGroup & node_group = getNodeGroup(node_group_name); ElementGroup & group = createElementGroup(name, dimension, node_group); group.fillFromNodeGroup(); } /* -------------------------------------------------------------------------- */ void GroupManager::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "GroupManager [" << std::endl; std::set node_group_seen; for (auto & group : iterateElementGroups()) { group.printself(stream, indent + 1); node_group_seen.insert(group.getNodeGroup().getName()); } for (auto & group : iterateNodeGroups()) { if (node_group_seen.find(group.getName()) == node_group_seen.end()) group.printself(stream, indent + 1); } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ UInt GroupManager::getNbElementGroups(UInt dimension) const { if (dimension == _all_dimensions) return element_groups.size(); auto it = element_groups.begin(); auto end = element_groups.end(); UInt count = 0; for (; it != end; ++it) count += (it->second->getDimension() == dimension); return count; } /* -------------------------------------------------------------------------- */ void GroupManager::checkAndAddGroups(DynamicCommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); UInt nb_node_group; buffer >> nb_node_group; AKANTU_DEBUG_INFO("Received " << nb_node_group << " node group names"); for (UInt ng = 0; ng < nb_node_group; ++ng) { std::string node_group_name; buffer >> node_group_name; if (node_groups.find(node_group_name) == node_groups.end()) { this->createNodeGroup(node_group_name); } AKANTU_DEBUG_INFO("Received node goup name: " << node_group_name); } UInt nb_element_group; buffer >> nb_element_group; AKANTU_DEBUG_INFO("Received " << nb_element_group << " element group names"); for (UInt eg = 0; eg < nb_element_group; ++eg) { std::string element_group_name; buffer >> element_group_name; std::string node_group_name; buffer >> node_group_name; UInt dim; buffer >> dim; AKANTU_DEBUG_INFO("Received element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with node group " << node_group_name); NodeGroup & node_group = *node_groups[node_group_name]; if (element_groups.find(element_group_name) == element_groups.end()) { this->createElementGroup(element_group_name, dim, node_group); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::fillBufferWithGroupNames( DynamicCommunicationBuffer & comm_buffer) const { AKANTU_DEBUG_IN(); // packing node group names; UInt nb_groups = this->node_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " node group names"); auto nnames_it = node_groups.begin(); auto nnames_end = node_groups.end(); for (; nnames_it != nnames_end; ++nnames_it) { std::string node_group_name = nnames_it->first; comm_buffer << node_group_name; AKANTU_DEBUG_INFO("Sending node goupe name: " << node_group_name); } // packing element group names with there associated node group name nb_groups = this->element_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " element group names"); auto gnames_it = this->element_groups.begin(); auto gnames_end = this->element_groups.end(); for (; gnames_it != gnames_end; ++gnames_it) { ElementGroup & element_group = *(gnames_it->second); std::string element_group_name = gnames_it->first; std::string node_group_name = element_group.getNodeGroup().getName(); UInt dim = element_group.getDimension(); comm_buffer << element_group_name; comm_buffer << node_group_name; comm_buffer << dim; AKANTU_DEBUG_INFO("Sending element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with the node group " << node_group_name); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::synchronizeGroupNames() { AKANTU_DEBUG_IN(); const Communicator & comm = mesh.getCommunicator(); Int nb_proc = comm.getNbProc(); Int my_rank = comm.whoAmI(); if (nb_proc == 1) return; if (my_rank == 0) { for (Int p = 1; p < nb_proc; ++p) { DynamicCommunicationBuffer recv_buffer; auto tag = Tag::genTag(p, 0, Tag::_ELEMENT_GROUP); comm.receive(recv_buffer, p, tag); AKANTU_DEBUG_INFO("Got " << printMemorySize(recv_buffer.size()) << " from proc " << p << " " << tag); this->checkAndAddGroups(recv_buffer); } DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); AKANTU_DEBUG_INFO("Initiating broadcast with " << printMemorySize(comm_buffer.size())); comm.broadcast(comm_buffer); } else { DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); auto tag = Tag::genTag(my_rank, 0, Tag::_ELEMENT_GROUP); AKANTU_DEBUG_INFO("Sending " << printMemorySize(comm_buffer.size()) << " to proc " << 0 << " " << tag); comm.send(comm_buffer, 0, tag); DynamicCommunicationBuffer recv_buffer; comm.broadcast(recv_buffer); AKANTU_DEBUG_INFO("Receiving broadcast with " << printMemorySize(recv_buffer.size())); this->checkAndAddGroups(recv_buffer); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const ElementGroup & GroupManager::getElementGroup(const std::string & name) const { auto it = element_groups.find(name); if (it == element_groups.end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::getElementGroup(const std::string & name) { auto it = element_groups.find(name); if (it == element_groups.end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const NodeGroup & GroupManager::getNodeGroup(const std::string & name) const { auto it = node_groups.find(name); if (it == node_groups.end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::getNodeGroup(const std::string & name) { auto it = node_groups.find(name); if (it == node_groups.end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ template void GroupManager::renameGroup(GroupsType & groups, const std::string & name, const std::string & new_name) { auto it = groups.find(name); if (it == groups.end()) { AKANTU_EXCEPTION("There are no group named " << name << " associated to the group manager: " << id); } auto && group_ptr = std::move(it->second); group_ptr->name = new_name; groups.erase(it); groups[new_name] = std::move(group_ptr); } /* -------------------------------------------------------------------------- */ void GroupManager::renameElementGroup(const std::string & name, const std::string & new_name) { renameGroup(element_groups, name, new_name); } /* -------------------------------------------------------------------------- */ void GroupManager::renameNodeGroup(const std::string & name, const std::string & new_name) { renameGroup(node_groups, name, new_name); } /* -------------------------------------------------------------------------- */ void GroupManager::copyElementGroup(const std::string & name, const std::string & new_name) { const auto & grp = getElementGroup(name); auto & new_grp = createElementGroup(new_name, grp.getDimension()); new_grp.getElements().copy(grp.getElements()); } /* -------------------------------------------------------------------------- */ void GroupManager::copyNodeGroup(const std::string & name, const std::string & new_name) { const auto & grp = getNodeGroup(name); auto & new_grp = createNodeGroup(new_name); new_grp.getNodes().copy(grp.getNodes()); } } // namespace akantu diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index 27bf2d6f1..e591273e2 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,351 +1,356 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Feb 07 2018 * * @brief Stores information relevent to the notion of element and nodes * groups. * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GROUP_MANAGER_HH__ #define __AKANTU_GROUP_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_iterators.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementGroup; class NodeGroup; class Mesh; class Element; class ElementSynchronizer; template class CommunicationBufferTemplated; namespace dumpers { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: using ElementGroups = std::map>; using NodeGroups = std::map>; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - GroupManager(const Mesh & mesh, const ID & id = "group_manager", + GroupManager(Mesh & mesh, const ID & id = "group_manager", const MemoryID & memory_id = 0); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ public: using node_group_iterator = NodeGroups::iterator; using element_group_iterator = ElementGroups::iterator; using const_node_group_iterator = NodeGroups::const_iterator; using const_element_group_iterator = ElementGroups::const_iterator; #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, \ param_in, param_out) \ [[deprecated( \ "use iterate(Element|Node)Groups or " \ "(element|node)GroupExists")]] inline BOOST_PP_CAT(BOOST_PP_CAT(const_, \ group_type), \ _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ }; \ \ [[deprecated("use iterate(Element|Node)Groups or " \ "(element|node)GroupExists")]] inline BOOST_PP_CAT(group_type, \ _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ } #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \ AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION( \ group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY()) AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name); public: decltype(auto) iterateNodeGroups() { return make_dereference_adaptor(make_values_adaptor(node_groups)); } decltype(auto) iterateNodeGroups() const { return make_dereference_adaptor(make_values_adaptor(node_groups)); } decltype(auto) iterateElementGroups() { return make_dereference_adaptor(make_values_adaptor(element_groups)); } decltype(auto) iterateElementGroups() const { return make_dereference_adaptor(make_values_adaptor(element_groups)); } /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* ------------------------------------------------------------------------ */ public: class ClusteringFilter { public: virtual bool operator()(const Element &) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create an element group and the associated node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension = _all_dimensions, bool replace_group = false); /* ------------------------------------------------------------------------ */ /// renames an element group void renameElementGroup(const std::string & name, const std::string & new_name); /// renames a node group void renameNodeGroup(const std::string & name, const std::string & new_name); /// copy an existing element group void copyElementGroup(const std::string & name, const std::string & new_name); /// copy an existing node group void copyNodeGroup(const std::string & name, const std::string & new_name); /* ------------------------------------------------------------------------ */ /// create a node group from another node group but filtered template NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// create an element group from another element group but filtered template ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy a node group void destroyNodeGroup(const std::string & group_name); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); // /// destroy all element groups and the associated node groups // void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// create groups based on values stored in a given mesh data template void createGroupsFromMeshData(const std::string & dataset_name); /// create boundaries group by a clustering algorithm \todo extend to parallel UInt createBoundaryGroupFromGeometry(); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix = "cluster", + const ClusteringStrategy & clustering_strategy_type = + ClusteringStrategy::_facets, const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix = "cluster", + const ClusteringStrategy & clustering_strategy_type = + ClusteringStrategy::_facets, const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix, + const ClusteringStrategy & clustering_strategy_type, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); /// register an elemental field to the given group name (overloading for /// ElementalPartionField) template class dump_type> std::shared_ptr createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// ElementalField) template class ret_type, template class, bool> class dump_type> std::shared_ptr createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template class dump_type> std::shared_ptr createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); template class ftype> std::shared_ptr createNodalField(const ftype * field, const std::string & group_name, UInt padding_size = 0); template class ftype> std::shared_ptr createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size); protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated & comm_buffer) const; /// take a buffer and create the missing groups localy void checkAndAddGroups(CommunicationBufferTemplated & buffer); /// register an elemental field to the given group name template inline std::shared_ptr createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem); /// register an elemental field to the given group name template inline std::shared_ptr createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: // AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; UInt getNbNodeGroups() { return node_groups.size(); }; bool elementGroupExists(const std::string & name) { return element_groups.find(name) != element_groups.end(); } bool nodeGroupExists(const std::string & name) { return node_groups.find(name) != node_groups.end(); } private: template void renameGroup(GroupsType & groups, const std::string & name, const std::string & new_name); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id to create element and node groups ID id; /// memory_id to create element and node groups MemoryID memory_id; /// list of the node groups managed NodeGroups node_groups; /// list of the element groups managed ElementGroups element_groups; /// Mesh to which the element belongs - const Mesh & mesh; + Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const GroupManager & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AKANTU_GROUP_MANAGER_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc index 0e640ce26..5548c714d 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc @@ -1,542 +1,543 @@ /** * @file fragment_manager.cc * * @author Aurelia Isabel Cuba Ramos * @author Marco Vocialta * * @date creation: Thu Jan 23 2014 * @date last modification: Tue Feb 20 2018 * * @brief Group manager to handle fragments * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fragment_manager.hh" #include "aka_iterators.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "material_cohesive.hh" #include "mesh_iterators.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ FragmentManager::FragmentManager(SolidMechanicsModelCohesive & model, bool dump_data, const ID & id, const MemoryID & memory_id) : GroupManager(model.getMesh(), id, memory_id), model(model), mass_center(0, model.getSpatialDimension(), "mass_center"), mass(0, model.getSpatialDimension(), "mass"), velocity(0, model.getSpatialDimension(), "velocity"), inertia_moments(0, model.getSpatialDimension(), "inertia_moments"), principal_directions( 0, model.getSpatialDimension() * model.getSpatialDimension(), "principal_directions"), quad_coordinates("quad_coordinates", id), mass_density("mass_density", id), nb_elements_per_fragment(0, 1, "nb_elements_per_fragment"), dump_data(dump_data) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); /// compute quadrature points' coordinates quad_coordinates.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(quad_coordinates, spatial_dimension, // spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getMesh().getNodes(), quad_coordinates); /// store mass density per quadrature point mass_density.initialize(mesh, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(mass_density, 1, spatial_dimension, // _not_ghost); storeMassDensityPerIntegrationPoint(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ class CohesiveElementFilter : public GroupManager::ClusteringFilter { public: CohesiveElementFilter(const SolidMechanicsModelCohesive & model, const Real max_damage = 1.) : model(model), is_unbroken(max_damage) {} bool operator()(const Element & el) const override { if (Mesh::getKind(el.type) == _ek_regular) return true; const Array & mat_indexes = model.getMaterialByElement(el.type, el.ghost_type); const Array & mat_loc_num = model.getMaterialLocalNumbering(el.type, el.ghost_type); const auto & mat = static_cast( model.getMaterial(mat_indexes(el.element))); UInt el_index = mat_loc_num(el.element); UInt nb_quad_per_element = model.getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(el.type, el.ghost_type); const Array & damage_array = mat.getDamage(el.type, el.ghost_type); AKANTU_DEBUG_ASSERT(nb_quad_per_element * el_index < damage_array.size(), "This quadrature point is out of range"); const Real * element_damage = damage_array.storage() + nb_quad_per_element * el_index; UInt unbroken_quads = std::count_if( element_damage, element_damage + nb_quad_per_element, is_unbroken); if (unbroken_quads > 0) return true; return false; } private: struct IsUnbrokenFunctor { IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {} bool operator()(const Real & x) { return x < max_damage; } const Real max_damage; }; const SolidMechanicsModelCohesive & model; const IsUnbrokenFunctor is_unbroken; }; /* -------------------------------------------------------------------------- */ void FragmentManager::buildFragments(Real damage_limit) { AKANTU_DEBUG_IN(); - ElementSynchronizer * cohesive_synchronizer = + auto * cohesive_synchronizer = const_cast(&model.getCohesiveSynchronizer()); if (cohesive_synchronizer) { cohesive_synchronizer->synchronize(model, SynchronizationTag::_smmc_damage); } auto & mesh_facets = const_cast(mesh.getMeshFacets()); - UInt spatial_dimension = model.getSpatialDimension(); + auto spatial_dimension = model.getSpatialDimension(); std::string fragment_prefix("fragment"); /// generate fragments global_nb_fragment = createClusters(spatial_dimension, mesh_facets, fragment_prefix, + ClusteringStrategy::_facets, CohesiveElementFilter(model, damage_limit)); nb_fragment = getNbElementGroups(spatial_dimension); fragment_index.resize(nb_fragment); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { auto name = std::get<0>(data).getName(); /// get fragment index std::string fragment_index_string = name.substr(fragment_prefix.size() + 1); std::get<1>(data) = std::stoul(fragment_index_string); } /// compute fragments' mass computeMass(); if (dump_data) { createDumpDataArray(fragment_index, "fragments", true); createDumpDataArray(mass, "fragments mass"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeMass() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// create a unit field per quadrature point, since to compute mass /// it's neccessary to integrate only density ElementTypeMapArray unit_field("unit_field", id); unit_field.initialize(model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); integrateFieldOnFragments(unit_field, mass); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeCenterOfMass() { AKANTU_DEBUG_IN(); /// integrate position multiplied by density integrateFieldOnFragments(quad_coordinates, mass_center); /// divide it by the fragments' mass Real * mass_storage = mass.storage(); Real * mass_center_storage = mass_center.storage(); UInt total_components = mass_center.size() * mass_center.getNbComponent(); for (UInt i = 0; i < total_components; ++i) mass_center_storage[i] /= mass_storage[i]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeVelocity() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// compute velocity per quadrature point ElementTypeMapArray velocity_field("velocity_field", id); velocity_field.initialize( model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getVelocity(), velocity_field); /// integrate on fragments integrateFieldOnFragments(velocity_field, velocity); /// divide it by the fragments' mass Real * mass_storage = mass.storage(); Real * velocity_storage = velocity.storage(); UInt total_components = velocity.size() * velocity.getNbComponent(); for (UInt i = 0; i < total_components; ++i) velocity_storage[i] /= mass_storage[i]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Given the distance @f$ \mathbf{r} @f$ between a quadrature point * and its center of mass, the moment of inertia is computed as \f[ * I_\mathrm{CM} = \mathrm{tr}(\mathbf{r}\mathbf{r}^\mathrm{T}) * \mathbf{I} - \mathbf{r}\mathbf{r}^\mathrm{T} \f] for more * information check Wikipedia * (http://en.wikipedia.org/wiki/Moment_of_inertia#Identities_for_a_skew-symmetric_matrix) * */ void FragmentManager::computeInertiaMoments() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); computeCenterOfMass(); /// compute local coordinates products with respect to the center of match ElementTypeMapArray moments_coords("moments_coords", id); moments_coords.initialize(model.getFEEngine(), _nb_component = spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); /// loop over fragments for (auto && data : zip(iterateElementGroups(), make_view(mass_center, spatial_dimension))) { const auto & el_list = std::get<0>(data).getElements(); auto & mass_center = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); auto & moments_coords_array = moments_coords(type); const auto & quad_coordinates_array = quad_coordinates(type); const auto & el_list_array = el_list(type); auto moments_begin = moments_coords_array.begin(spatial_dimension, spatial_dimension); auto quad_coordinates_begin = quad_coordinates_array.begin(spatial_dimension); Vector relative_coords(spatial_dimension); for (UInt el = 0; el < el_list_array.size(); ++el) { UInt global_el = el_list_array(el); /// loop over quadrature points for (UInt q = 0; q < nb_quad_per_element; ++q) { UInt global_q = global_el * nb_quad_per_element + q; Matrix moments_matrix = moments_begin[global_q]; const Vector & quad_coord_vector = quad_coordinates_begin[global_q]; /// to understand this read the documentation written just /// before this function relative_coords = quad_coord_vector; relative_coords -= mass_center; moments_matrix.outerProduct(relative_coords, relative_coords); Real trace = moments_matrix.trace(); moments_matrix *= -1.; moments_matrix += Matrix::eye(spatial_dimension, trace); } } } } /// integrate moments Array integrated_moments(global_nb_fragment, spatial_dimension * spatial_dimension); integrateFieldOnFragments(moments_coords, integrated_moments); /// compute and store principal moments inertia_moments.resize(global_nb_fragment); principal_directions.resize(global_nb_fragment); auto integrated_moments_it = integrated_moments.begin(spatial_dimension, spatial_dimension); auto inertia_moments_it = inertia_moments.begin(spatial_dimension); auto principal_directions_it = principal_directions.begin(spatial_dimension, spatial_dimension); for (UInt frag = 0; frag < global_nb_fragment; ++frag, ++integrated_moments_it, ++inertia_moments_it, ++principal_directions_it) { integrated_moments_it->eig(*inertia_moments_it, *principal_directions_it); } if (dump_data) createDumpDataArray(inertia_moments, "moments of inertia"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeAllData() { AKANTU_DEBUG_IN(); buildFragments(); computeVelocity(); computeInertiaMoments(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::storeMassDensityPerIntegrationPoint() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : mesh.elementTypes(_spatial_dimension = spatial_dimension, _element_kind = _ek_regular)) { Array & mass_density_array = mass_density(type); UInt nb_element = mesh.getNbElement(type); UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); mass_density_array.resize(nb_element * nb_quad_per_element); const Array & mat_indexes = model.getMaterialByElement(type); Real * mass_density_it = mass_density_array.storage(); /// store mass_density for each element and quadrature point for (UInt el = 0; el < nb_element; ++el) { Material & mat = model.getMaterial(mat_indexes(el)); for (UInt q = 0; q < nb_quad_per_element; ++q, ++mass_density_it) *mass_density_it = mat.getRho(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::integrateFieldOnFragments( ElementTypeMapArray & field, Array & output) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); UInt nb_component = output.getNbComponent(); /// integration part output.resize(global_nb_fragment); output.clear(); auto output_begin = output.begin(nb_component); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & el_list = std::get<0>(data).getElements(); auto fragment_index = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); const Array & density_array = mass_density(type); Array & field_array = field(type); const Array & elements = el_list(type); /// generate array to be integrated by filtering fragment's elements Array integration_array(elements.size() * nb_quad_per_element, nb_component); auto field_array_begin = field_array.begin_reinterpret( nb_quad_per_element, nb_component, field_array.size() / nb_quad_per_element); auto density_array_begin = density_array.begin_reinterpret( nb_quad_per_element, density_array.size() / nb_quad_per_element); for (auto && data : enumerate(make_view( integration_array, nb_quad_per_element, nb_component))) { UInt global_el = elements(std::get<0>(data)); auto & int_array = std::get<1>(data); int_array = field_array_begin[global_el]; /// multiply field by density const Vector & density_vector = density_array_begin[global_el]; for (UInt i = 0; i < nb_quad_per_element; ++i) { for (UInt j = 0; j < nb_component; ++j) { int_array(i, j) *= density_vector(i); } } } /// integrate the field over the fragment Array integrated_array(elements.size(), nb_component); model.getFEEngine().integrate(integration_array, integrated_array, nb_component, type, _not_ghost, elements); /// sum over all elements and store the result Vector output_tmp(output_begin[fragment_index]); output_tmp += std::accumulate(integrated_array.begin(nb_component), integrated_array.end(nb_component), Vector(nb_component)); } } /// sum output over all processors const Communicator & comm = mesh.getCommunicator(); comm.allReduce(output, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeNbElementsPerFragment() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); nb_elements_per_fragment.resize(global_nb_fragment); nb_elements_per_fragment.clear(); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & el_list = std::get<0>(data).getElements(); auto fragment_index = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_element = el_list(type).size(); nb_elements_per_fragment(fragment_index) += nb_element; } } /// sum values over all processors const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_elements_per_fragment, SynchronizerOperation::_sum); if (dump_data) createDumpDataArray(nb_elements_per_fragment, "elements per fragment"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FragmentManager::createDumpDataArray(Array & data, std::string name, bool fragment_index_output) { AKANTU_DEBUG_IN(); if (data.size() == 0) return; auto & mesh_not_const = const_cast(mesh); auto && spatial_dimension = mesh.getSpatialDimension(); auto && nb_component = data.getNbComponent(); auto && data_begin = data.begin(nb_component); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & fragment = std::get<0>(data); auto fragment_idx = std::get<1>(data); /// loop over cluster types for (auto & type : fragment.elementTypes(spatial_dimension)) { /// init mesh data auto & mesh_data = mesh_not_const.getDataPointer( name, type, _not_ghost, nb_component); auto mesh_data_begin = mesh_data.begin(nb_component); /// fill mesh data for (const auto & elem : fragment.getElements(type)) { Vector md_tmp = mesh_data_begin[elem]; if (fragment_index_output) { md_tmp(0) = fragment_idx; } else { md_tmp = data_begin[fragment_idx]; } } } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh index 23960075a..447b56471 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh @@ -1,265 +1,265 @@ /** * @file material_cohesive_linear_inline_impl.hh * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Apr 22 2015 * @date last modification: Wed Feb 21 2018 * * @brief Inline functions of the MaterialCohesiveLinear * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template inline Real MaterialCohesiveLinear::computeEffectiveNorm( const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_traction) const { normal_traction.mul(stress, normal); Real normal_contrib = normal_traction.dot(normal); /// in 3D tangential components must be summed Real tangent_contrib = 0; if (dim == 2) { Real tangent_contrib_tmp = normal_traction.dot(tangent); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } else if (dim == 3) { for (UInt s = 0; s < dim - 1; ++s) { const Vector tangent_v(tangent.storage() + s * dim, dim); Real tangent_contrib_tmp = normal_traction.dot(tangent_v); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } } tangent_contrib = std::sqrt(tangent_contrib); normal_contrib = std::max(Real(0.), normal_contrib); return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv); } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTractionOnQuad( Vector & traction, Vector & opening, const Vector & normal, Real & delta_max, const Real & delta_c, const Vector & insertion_stress, const Real & sigma_c, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_traction, Vector & contact_opening) { /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm / delta_c < -Math::getTolerance(); // penetration = normal_opening_norm < 0.; - if (this->contact_after_breaking == false && + if (not this->contact_after_breaking && Math::are_float_equal(damage, 1.)) penetration = false; if (penetration) { /// use penalty coefficient in case of penetration contact_traction = normal_opening; contact_traction *= this->penalty; contact_opening = normal_opening; /// don't consider penetration contribution for delta opening = tangential_opening; normal_opening.clear(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction.clear(); contact_opening.clear(); } delta = std::sqrt(delta); /// update maximum displacement and damage delta_max = std::max(delta_max, delta); damage = std::min(delta_max / delta_c, Real(1.)); /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ if (Math::are_float_equal(damage, 1.)) traction.clear(); else if (Math::are_float_equal(damage, 0.)) { if (penetration) traction.clear(); else traction = insertion_stress; } else { traction = tangential_opening; traction *= this->beta2_kappa; traction += normal_opening; AKANTU_DEBUG_ASSERT(delta_max != 0., "Division by zero, tolerance might be too low"); traction *= sigma_c / delta_max * (1. - damage); } } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTangentTractionOnQuad( Matrix & tangent, Real & delta_max, const Real & delta_c, const Real & sigma_c, Vector & opening, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_opening) { /** * During the update of the residual the interpenetrations are * stored in the array "contact_opening", therefore, in the case * of penetration, in the array "opening" there are only the * tangential components. */ opening += contact_opening; /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm < 0.0; if (this->contact_after_breaking == false && Math::are_float_equal(damage, 1.)) penetration = false; Real derivative = 0; // derivative = d(t/delta)/ddelta Real t = 0; Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); if (penetration) { /// stiffness in compression given by the penalty parameter tangent += n_outer_n; tangent *= penalty; opening = tangential_opening; normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; } else { delta += normal_opening_norm * normal_opening_norm; } delta = std::sqrt(delta); /** * Delta has to be different from 0 to have finite values of * tangential stiffness. At the element insertion, delta = * 0. Therefore, a fictictious value is defined, for the * evaluation of the first value of K. */ if (delta < Math::getTolerance()) delta = delta_c / 1000.; if (delta >= delta_max) { if (delta <= delta_c) { derivative = -sigma_c / (delta * delta); t = sigma_c * (1 - delta / delta_c); } else { derivative = 0.; t = 0.; } } else if (delta < delta_max) { Real tmax = sigma_c * (1 - delta_max / delta_c); t = tmax / delta_max * delta; } /// computation of the derivative of the constitutive law (dT/ddelta) Matrix I(spatial_dimension, spatial_dimension); I.eye(this->beta2_kappa); Matrix nn(n_outer_n); nn *= (1. - this->beta2_kappa); nn += I; nn *= t / delta; Vector t_tilde(normal_opening); t_tilde *= (1. - this->beta2_kappa2); Vector mm(opening); mm *= this->beta2_kappa2; t_tilde += mm; Vector t_hat(normal_opening); t_hat += this->beta2_kappa * tangential_opening; Matrix prov(spatial_dimension, spatial_dimension); prov.outerProduct(t_hat, t_tilde); prov *= derivative / delta; prov += nn; Matrix prov_t = prov.transpose(); tangent += prov_t; } /* -------------------------------------------------------------------------- */ } // namespace akantu /* -------------------------------------------------------------------------- */ #endif //__AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH__