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