diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh
index 29de7f949..95abb457a 100644
--- a/src/common/aka_common.hh
+++ b/src/common/aka_common.hh
@@ -1,698 +1,714 @@
/**
* @file aka_common.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <boost/preprocessor.hpp>
#include <limits>
#include <list>
#include <memory>
#include <string>
#include <type_traits>
#include <unordered_map>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Constants */
/* -------------------------------------------------------------------------- */
namespace {
[[gnu::unused]] constexpr UInt _all_dimensions{
std::numeric_limits<UInt>::max()};
#ifdef AKANTU_NDEBUG
[[gnu::unused]] constexpr Real REAL_INIT_VALUE{0.};
#else
[[gnu::unused]] constexpr Real REAL_INIT_VALUE{
std::numeric_limits<Real>::quiet_NaN()};
#endif
} // namespace
/* -------------------------------------------------------------------------- */
/* Common types */
/* -------------------------------------------------------------------------- */
using ID = std::string;
} // 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) \
(contact_mechanics_model) \
(coupler_solid_contact) \
(coupler_solid_cohesive_contact) \
(phase_field_model) \
(coupler_solid_phasefield)
// clang-format on
/// enum ModelType defines which type of physics is solved
AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES)
AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES)
#else
enum class ModelType {
model,
solid_mechanics_model,
solid_mechanics_model_cohesive,
heat_transfer_model,
structural_mechanics_model,
embedded_model,
};
#endif
/// enum AnalysisMethod type of solving method used to solve the equation of
/// motion
enum AnalysisMethod {
_static = 0,
_implicit_dynamic = 1,
_explicit_lumped_mass = 2,
_explicit_lumped_capacity = 2,
_explicit_consistent_mass = 3,
_explicit_contact = 4,
_implicit_contact = 5
};
/// enum DOFSupportType defines which kind of dof that can exists
enum DOFSupportType { _dst_nodal, _dst_generic };
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_NON_LINEAR_SOLVER_TYPES \
(linear) \
(newton_raphson) \
(newton_raphson_modified) \
(lumped) \
(gmres) \
(bfgs) \
(cg) \
(newton_raphson_contact) \
(auto)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType,
AKANTU_NON_LINEAR_SOLVER_TYPES)
AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType,
AKANTU_NON_LINEAR_SOLVER_TYPES)
#else
/// Type of non linear resolution available in akantu
enum class NonLinearSolverType {
_linear, ///< No non linear convergence loop
_newton_raphson, ///< Regular Newton-Raphson
_newton_raphson_modified, ///< Newton-Raphson with initial tangent
_lumped, ///< Case of lumped mass or equivalent matrix
_gmres,
_bfgs,
_cg,
- _newton_raphson_contact, ///< Regular Newton-Raphson modified
- /// for contact problem
- _auto, ///< This will take a default value that make sense in case of
- /// model::getNewSolver
+ _newton_raphson_contact, ///< Regular Newton-Raphson modified
+ /// for contact problem
+ _auto, ///< This will take a default value that make sense in case of
+ /// model::getNewSolver
};
#endif
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_TIME_STEP_SOLVER_TYPE \
(static) \
(dynamic) \
(dynamic_lumped) \
(not_defined)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType,
AKANTU_TIME_STEP_SOLVER_TYPE)
AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE)
#else
/// Type of time stepping solver
enum class TimeStepSolverType {
_static, ///< Static solution
_dynamic, ///< Dynamic solver
_dynamic_lumped, ///< Dynamic solver with lumped mass
_not_defined, ///< For not defined cases
};
#endif
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_INTEGRATION_SCHEME_TYPE \
(pseudo_time) \
(forward_euler) \
(trapezoidal_rule_1) \
(backward_euler) \
(central_difference) \
(fox_goodwin) \
(trapezoidal_rule_2) \
(linear_acceleration) \
(newmark_beta) \
(generalized_trapezoidal)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType,
AKANTU_INTEGRATION_SCHEME_TYPE)
AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType,
AKANTU_INTEGRATION_SCHEME_TYPE)
#else
/// Type of integration scheme
enum class IntegrationSchemeType {
_pseudo_time, ///< Pseudo Time
_forward_euler, ///< GeneralizedTrapezoidal(0)
_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2)
_backward_euler, ///< GeneralizedTrapezoidal(1)
_central_difference, ///< NewmarkBeta(0, 1/2)
_fox_goodwin, ///< NewmarkBeta(1/6, 1/2)
_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2)
_linear_acceleration, ///< NewmarkBeta(1/3, 1/2)
_newmark_beta, ///< generic NewmarkBeta with user defined
/// alpha and beta
_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user
/// defined alpha
};
#endif
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_SOLVE_CONVERGENCE_CRITERIA \
(residual) \
(solution) \
(residual_mass_wgh)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria,
AKANTU_SOLVE_CONVERGENCE_CRITERIA)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria,
AKANTU_SOLVE_CONVERGENCE_CRITERIA)
AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria,
AKANTU_SOLVE_CONVERGENCE_CRITERIA)
#else
/// enum SolveConvergenceCriteria different convergence criteria
enum class SolveConvergenceCriteria {
_residual, ///< Use residual to test the convergence
_solution, ///< Use solution to test the convergence
_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to
///< testb
};
#endif
/// enum CohesiveMethod type of insertion of cohesive elements
enum CohesiveMethod { _intrinsic, _extrinsic };
/// @enum MatrixType type of sparse matrix used
enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined };
/// @enum Type of contact detection
-enum DetectionType { _explicit, _implicit};
+enum DetectionType { _explicit, _implicit };
+
+#if !defined(DOXYGEN)
+// clang-format off
+#define AKANTU_CONTACT_STATE \
+ (no_contact) \
+ (stick) \
+ (slip)
+// clang-format on
+AKANTU_CLASS_ENUM_DECLARE(ContactState,
+ AKANTU_CONTACT_STATE)
+AKANTU_CLASS_ENUM_OUTPUT_STREAM(ContactState,
+ AKANTU_CONTACT_STATE)
+AKANTU_CLASS_ENUM_INPUT_STREAM(ContactState,
+ AKANTU_CONTACT_STATE)
+#else
/// @enum no contact or stick or slip state
enum class ContactState {
_no_contact = 0,
_stick = 1,
_slip = 2,
};
-
+#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) \
(ce_insertion_order) \
(gm_clusters) \
(htm_temperature) \
(htm_gradient_temperature) \
(htm_phi) \
(htm_gradient_phi) \
(pfm_damage) \
(pfm_driving) \
(pfm_history) \
(pfm_energy) \
(csp_damage) \
(csp_strain) \
(mnl_for_average) \
(mnl_weight) \
(nh_criterion) \
(test) \
(user_1) \
(user_2) \
(material_id) \
(for_dump) \
(cf_nodal) \
(cf_incr) \
(solver_solution)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG)
#else
/// @enum SynchronizationTag type of synchronizations
enum class SynchronizationTag {
//--- Generic tags ---
_whatever,
_update,
_ask_nodes,
_size,
//--- SolidMechanicsModel tags ---
_smm_mass, ///< synchronization of the SolidMechanicsModel.mass
_smm_for_gradu, ///< synchronization of the
/// SolidMechanicsModel.displacement
_smm_boundary, ///< synchronization of the boundary, forces, velocities
/// and displacement
_smm_uv, ///< synchronization of the nodal velocities and displacement
_smm_res, ///< synchronization of the nodal residual
_smm_init_mat, ///< synchronization of the data to initialize materials
_smm_stress, ///< synchronization of the stresses to compute the
///< internal
/// forces
_smmc_facets, ///< synchronization of facet data to setup facet synch
_smmc_facets_conn, ///< synchronization of facet global connectivity
_smmc_facets_stress, ///< synchronization of facets' stress to setup
///< facet
/// synch
_smmc_damage, ///< synchronization of damage
// --- GlobalIdsUpdater tags ---
_giu_global_conn, ///< synchronization of global connectivities
// --- CohesiveElementInserter tags ---
_ce_groups, ///< synchronization of cohesive element insertion depending
/// on facet groups
_ce_insertion_order, ///< synchronization of the order of insertion of
/// cohesive elements
// --- 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
// --- PhaseFieldModel tags ---
- _pfm_damage, ///< synchronization of the nodal damage
- _pfm_driving, ///< synchronization of the driving forces to
- /// compute the internal
- _pfm_history, ///< synchronization of the damage history to
- /// compute the internal
- _pfm_energy, ///< synchronization of the damage energy
- /// density to compute the internal
+ _pfm_damage, ///< synchronization of the nodal damage
+ _pfm_driving, ///< synchronization of the driving forces to
+ /// compute the internal
+ _pfm_history, ///< synchronization of the damage history to
+ /// compute the internal
+ _pfm_energy, ///< synchronization of the damage energy
+ /// density to compute the internal
// --- CouplerSolidPhaseField tags ---
- _csp_damage, ///< synchronization of the damage from phase
- /// model to solid model
- _csp_strain, ///< synchronization of the strain from solid
- /// model to phase model
-
+ _csp_damage, ///< synchronization of the damage from phase
+ /// model to solid model
+ _csp_strain, ///< synchronization of the strain from solid
+ /// model to phase model
+
// --- LevelSet tags ---
_htm_phi, ///< synchronization of the nodal level set value phi
_htm_gradient_phi, ///< synchronization of the element gradient phi
_pfm_damage,
_pfm_gradient_damage,
-
+
//--- Material non local ---
_mnl_for_average, ///< synchronization of data to average in non local
/// material
_mnl_weight, ///< synchronization of data for the weight computations
// --- NeighborhoodSynchronization tags ---
_nh_criterion,
// --- General tags ---
_test, ///< Test tag
_user_1, ///< tag for user simulations
_user_2, ///< tag for user simulations
_material_id, ///< synchronization of the material ids
_for_dump, ///< everything that needs to be synch before dump
// --- Contact & Friction ---
_cf_nodal, ///< synchronization of disp, velo, and current position
_cf_incr, ///< synchronization of increment
// --- Solver tags ---
_solver_solution ///< synchronization of the solution obained with the
/// PETSc solver
};
#endif
/// @enum GhostType type of ghost
enum GhostType {
_not_ghost = 0,
_ghost = 1,
_casper // not used but a real cute ghost
};
/// Define the flag that can be set to a node
enum class NodeFlag : std::uint8_t {
_normal = 0x00,
_distributed = 0x01,
_master = 0x03,
_slave = 0x05,
_pure_ghost = 0x09,
_shared_mask = 0x0F,
_periodic = 0x10,
_periodic_master = 0x30,
_periodic_slave = 0x50,
_periodic_mask = 0xF0,
_local_master_mask = 0xCC, // ~(_master & _periodic_mask)
};
inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) {
using under = std::underlying_type_t<NodeFlag>;
return NodeFlag(under(a) & under(b));
}
inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) {
using under = std::underlying_type_t<NodeFlag>;
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<NodeFlag>;
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<GhostType_def>;
namespace {
constexpr ghost_type_t ghost_types{_casper};
}
/// 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 const auto & get##name() const { \
if (not(ptr)) { \
AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \
} \
return (*(ptr)); \
}
#define AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(name, ptr) \
inline auto & get##name() { \
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<type> & get##name(const support & el_type, \
GhostType ghost_type = _not_ghost) \
con { /* NOLINT */ \
return variable(el_type, ghost_type); \
} // NOLINT
#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 <typename T> std::string printMemorySize(UInt size);
/* -------------------------------------------------------------------------- */
struct TensorTrait {};
struct TensorProxyTrait {};
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* Type traits */
/* -------------------------------------------------------------------------- */
namespace aka {
/* ------------------------------------------------------------------------ */
template <typename T> using is_tensor = std::is_base_of<akantu::TensorTrait, T>;
template <typename T>
using is_tensor_proxy = std::is_base_of<akantu::TensorProxyTrait, T>;
/* ------------------------------------------------------------------------ */
template <typename T> using is_scalar = std::is_arithmetic<T>;
/* ------------------------------------------------------------------------ */
template <typename R, typename T,
std::enable_if_t<std::is_reference<T>::value> * = nullptr>
bool is_of_type(T && t) {
return (
dynamic_cast<std::add_pointer_t<
std::conditional_t<std::is_const<std::remove_reference_t<T>>::value,
std::add_const_t<R>, R>>>(&t) != nullptr);
}
/* -------------------------------------------------------------------------- */
template <typename R, typename T> bool is_of_type(std::unique_ptr<T> & t) {
return (
dynamic_cast<std::add_pointer_t<
std::conditional_t<std::is_const<T>::value, std::add_const_t<R>, R>>>(
t.get()) != nullptr);
}
/* ------------------------------------------------------------------------ */
template <typename R, typename T,
std::enable_if_t<std::is_reference<T>::value> * = nullptr>
decltype(auto) as_type(T && t) {
static_assert(
disjunction<
std::is_base_of<std::decay_t<T>, std::decay_t<R>>, // down-cast
std::is_base_of<std::decay_t<R>, std::decay_t<T>> // up-cast
>::value,
"Type T and R are not valid for a as_type conversion");
return dynamic_cast<std::add_lvalue_reference_t<
std::conditional_t<std::is_const<std::remove_reference_t<T>>::value,
std::add_const_t<R>, R>>>(t);
}
/* -------------------------------------------------------------------------- */
template <typename R, typename T,
std::enable_if_t<std::is_pointer<T>::value> * = nullptr>
decltype(auto) as_type(T && t) {
return &as_type<R>(*t);
}
/* -------------------------------------------------------------------------- */
template <typename R, typename T>
decltype(auto) as_type(const std::shared_ptr<T> & t) {
return std::dynamic_pointer_cast<R>(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 <typename a, typename b> struct hash<std::pair<a, b>> {
hash() = default;
size_t operator()(const std::pair<a, b> & p) const {
size_t seed = ah(p.first);
return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) +
(seed >> 2);
}
private:
const hash<a> ah{};
const hash<b> bh{};
};
} // namespace std
#endif // AKANTU_COMMON_HH_
diff --git a/src/common/aka_math.hh b/src/common/aka_math.hh
index 7e1121d91..a1df46c8e 100644
--- a/src/common/aka_math.hh
+++ b/src/common/aka_math.hh
@@ -1,283 +1,284 @@
/**
* @file aka_math.hh
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Marion Estelle Chambart <marion.chambart@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Aug 04 2010
* @date last modification: Mon Sep 11 2017
*
* @brief mathematical operations
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
#include <utility>
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_AKA_MATH_H_
#define AKANTU_AKA_MATH_H_
namespace akantu {
/* -------------------------------------------------------------------------- */
namespace Math {
/// tolerance for functions that need one
extern Real tolerance; // NOLINT
/* ------------------------------------------------------------------------ */
/* Matrix algebra */
/* ------------------------------------------------------------------------ */
/// @f$ y = A*x @f$
void matrix_vector(UInt m, UInt n, const Array<Real> & A,
const Array<Real> & x, Array<Real> & y, Real alpha = 1.);
/// @f$ y = A*x @f$
inline void matrix_vector(UInt m, UInt n, Real * A, Real * x, Real * y,
Real alpha = 1.);
/// @f$ y = A^t*x @f$
inline void matrixt_vector(UInt m, UInt n, Real * A, Real * x, Real * y,
Real alpha = 1.);
/// @f$ C = A*B @f$
void matrix_matrix(UInt m, UInt n, UInt k, const Array<Real> & A,
const Array<Real> & B, Array<Real> & C, Real alpha = 1.);
/// @f$ C = A*B^t @f$
void matrix_matrixt(UInt m, UInt n, UInt k, const Array<Real> & A,
const Array<Real> & B, Array<Real> & C, Real alpha = 1.);
/// @f$ C = A*B @f$
inline void matrix_matrix(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
/// @f$ C = A^t*B @f$
inline void matrixt_matrix(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
/// @f$ C = A*B^t @f$
inline void matrix_matrixt(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
/// @f$ C = A^t*B^t @f$
inline void matrixt_matrixt(UInt m, UInt n, UInt k, Real * A, Real * B,
Real * C, Real alpha = 1.);
template <bool tr_A, bool tr_B>
inline void matMul(UInt m, UInt n, UInt k, Real alpha, Real * A, Real * B,
Real beta, Real * C);
template <bool tr_A>
inline void matVectMul(UInt m, UInt n, Real alpha, Real * A, Real * x,
Real beta, Real * y);
inline void aXplusY(UInt n, Real alpha, Real * x, Real * y);
inline void matrix33_eigenvalues(Real * A, Real * Adiag);
inline void matrix22_eigenvalues(Real * A, Real * Adiag);
template <UInt dim> inline void eigenvalues(Real * A, Real * d);
/// solve @f$ A x = \Lambda x @f$ and return d and V such as @f$ A V[i:] =
/// d[i] V[i:]@f$
template <typename T> void matrixEig(UInt n, T * A, T * d, T * V = nullptr);
/// determinent of a 2x2 matrix
Real det2(const Real * mat);
/// determinent of a 3x3 matrix
Real det3(const Real * mat);
/// determinent of a nxn matrix
template <UInt n> Real det(const Real * mat);
/// determinent of a nxn matrix
template <typename T> T det(UInt n, const T * A);
/// inverse a nxn matrix
template <UInt n> inline void inv(const Real * A, Real * inv);
/// inverse a nxn matrix
template <typename T> inline void inv(UInt n, const T * A, T * inv);
/// inverse a 3x3 matrix
inline void inv3(const Real * mat, Real * inv);
/// inverse a 2x2 matrix
inline void inv2(const Real * mat, Real * inv);
/// solve A x = b using a LU factorization
template <typename T>
inline void solve(UInt n, const T * A, T * x, const T * b);
/// return the double dot product between 2 tensors in 2d
inline Real matrixDoubleDot22(Real * A, Real * B);
/// return the double dot product between 2 tensors in 3d
inline Real matrixDoubleDot33(Real * A, Real * B);
/// extension of the double dot product to two 2nd order tensor in dimension n
inline Real matrixDoubleDot(UInt n, Real * A, Real * B);
/* ------------------------------------------------------------------------ */
/* Array algebra */
/* ------------------------------------------------------------------------ */
/// vector cross product
inline void vectorProduct3(const Real * v1, const Real * v2, Real * res);
/// normalize a vector
inline void normalize2(Real * v);
/// normalize a vector
inline void normalize3(Real * v);
/// return norm of a 2-vector
inline Real norm2(const Real * v);
/// return norm of a 3-vector
inline Real norm3(const Real * v);
/// return norm of a vector
inline Real norm(UInt n, const Real * v);
/// return the dot product between 2 vectors in 2d
inline Real vectorDot2(const Real * v1, const Real * v2);
/// return the dot product between 2 vectors in 3d
inline Real vectorDot3(const Real * v1, const Real * v2);
/// return the dot product between 2 vectors
inline Real vectorDot(Real * v1, Real * v2, UInt n);
/* ------------------------------------------------------------------------ */
/* Geometry */
/* ------------------------------------------------------------------------ */
/// compute normal a normal to a vector
inline void normal2(const Real * vec, Real * normal);
/// compute normal a normal to a vector
inline void normal3(const Real * vec1, const Real * vec2, Real * normal);
/// compute the tangents to an array of normal vectors
void compute_tangents(const Array<Real> & normals, Array<Real> & tangents);
/// distance in 2D between x and y
inline Real distance_2d(const Real * x, const Real * y);
/// distance in 3D between x and y
inline Real distance_3d(const Real * x, const Real * y);
/// radius of the in-circle of a triangle in 2d space
- static inline Real triangle_inradius_2d(const Vector<Real> & coord1, const Vecotr<Real> & coord2,
- const Vector<Real> & coord3);
+ static inline Real triangle_inradius(const Vector<Real> & coord1,
+ const Vector<Real> & coord2,
+ const Vector<Real> & coord3);
/// radius of the in-circle of a tetrahedron
inline Real tetrahedron_inradius(const Real * coord1, const Real * coord2,
const Real * coord3, const Real * coord4);
/// volume of a tetrahedron
inline Real tetrahedron_volume(const Real * coord1, const Real * coord2,
const Real * coord3, const Real * coord4);
/// compute the barycenter of n points
inline void barycenter(const Real * coord, UInt nb_points,
UInt spatial_dimension, Real * barycenter);
/// vector between x and y
inline void vector_2d(const Real * x, const Real * y, Real * res);
/// vector pointing from x to y in 3 spatial dimension
inline void vector_3d(const Real * x, const Real * y, Real * res);
/// test if two scalar are equal within a given tolerance
inline bool are_float_equal(Real x, Real y);
/// test if two vectors are equal within a given tolerance
inline bool are_vector_equal(UInt n, Real * x, Real * y);
#ifdef isnan
#error \
"You probably included <math.h> which is incompatible with aka_math please use\
<cmath> or add a \"#undef isnan\" before akantu includes"
#endif
/// test if a real is a NaN
inline bool isnan(Real x);
/// test if the line x and y intersects each other
inline bool intersects(Real x_min, Real x_max, Real y_min, Real y_max);
/// test if a is in the range [x_min, x_max]
inline bool is_in_range(Real a, Real x_min, Real x_max);
inline Real getTolerance() { return Math::tolerance; }
inline void setTolerance(Real tol) { Math::tolerance = tol; }
template <UInt p, typename T> inline T pow(T x);
template <class T1, class T2,
std::enable_if_t<std::is_integral<T1>::value and
std::is_integral<T2>::value> * = nullptr>
inline Real kronecker(T1 i, T2 j) {
return static_cast<Real>(i == j);
}
/// reduce all the values of an array, the summation is done in place and the
/// array is modified
Real reduce(Array<Real> & array);
class NewtonRaphson {
public:
NewtonRaphson(Real tolerance, Real max_iteration)
: tolerance(tolerance), max_iteration(max_iteration) {}
template <class Functor> Real solve(const Functor & funct, Real x_0);
private:
Real tolerance;
Real max_iteration;
};
struct NewtonRaphsonFunctor {
explicit NewtonRaphsonFunctor(const std::string & name) : name(name) {}
virtual ~NewtonRaphsonFunctor() = default;
NewtonRaphsonFunctor(const NewtonRaphsonFunctor & other) = default;
NewtonRaphsonFunctor(NewtonRaphsonFunctor && other) noexcept = default;
NewtonRaphsonFunctor &
operator=(const NewtonRaphsonFunctor & other) = default;
NewtonRaphsonFunctor &
operator=(NewtonRaphsonFunctor && other) noexcept = default;
virtual Real f(Real x) const = 0;
virtual Real f_prime(Real x) const = 0;
std::string name;
};
} // namespace Math
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "aka_math_tmpl.hh"
#endif /* AKANTU_AKA_MATH_H_ */
diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh
index 335b7228f..61a75f471 100644
--- a/src/common/aka_math_tmpl.hh
+++ b/src/common/aka_math_tmpl.hh
@@ -1,838 +1,839 @@
/**
* @file aka_math_tmpl.hh
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Mathilde Radiguet <mathilde.radiguet@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Aug 04 2010
* @date last modification: Tue Feb 20 2018
*
* @brief Implementation of the inline functions of the math toolkit
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_blas_lapack.hh"
#include "aka_math.hh"
+#include "aka_types.hh"
/* -------------------------------------------------------------------------- */
#include <cmath>
#include <typeinfo>
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace Math {
/* ------------------------------------------------------------------------ */
inline void matrix_vector(UInt im, UInt in,
Real * A, // NOLINT(readability-non-const-parameter)
Real * x, // NOLINT(readability-non-const-parameter)
Real * y, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// y = alpha*op(A)*x + beta*y
char tran_A = 'N';
int incx = 1;
int incy = 1;
double beta = 0.;
int m = im;
int n = in;
aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy);
#else
std::fill_n(y, im, 0.);
for (UInt i = 0; i < im; ++i) {
for (UInt j = 0; j < in; ++j) {
y[i] += A[i + j * im] * x[j];
}
y[i] *= alpha;
}
#endif
}
/* ------------------------------------------------------------------------ */
inline void
matrixt_vector(UInt im, UInt in,
Real * A, // NOLINT(readability-non-const-parameter)
Real * x, // NOLINT(readability-non-const-parameter)
Real * y, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// y = alpha*op(A)*x + beta*y
char tran_A = 'T';
int incx = 1;
int incy = 1;
double beta = 0.;
int m = im;
int n = in;
aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy);
#else
std::fill_n(y, in, 0.);
for (UInt i = 0; i < im; ++i) {
for (UInt j = 0; j < in; ++j) {
y[j] += A[j * im + i] * x[i];
}
y[i] *= alpha;
}
#endif
}
/* ------------------------------------------------------------------------ */
inline void matrix_matrix(UInt im, UInt in, UInt ik,
Real * A, // NOLINT(readability-non-const-parameter)
Real * B, // NOLINT(readability-non-const-parameter)
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'N';
char trans_b = 'N';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &k, &beta, C,
&m);
#else
std::fill_n(C, im * in, 0.);
for (UInt j = 0; j < in; ++j) {
UInt _jb = j * ik;
UInt _jc = j * im;
for (UInt i = 0; i < im; ++i) {
for (UInt l = 0; l < ik; ++l) {
UInt _la = l * im;
C[i + _jc] += A[i + _la] * B[l + _jb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* ------------------------------------------------------------------------ */
inline void
matrixt_matrix(UInt im, UInt in, UInt ik,
Real * A, // NOLINT(readability-non-const-parameter)
Real * B, // NOLINT(readability-non-const-parameter)
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'T';
char trans_b = 'N';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &k, &beta, C,
&m);
#else
std::fill_n(C, im * in, 0.);
for (UInt j = 0; j < in; ++j) {
UInt _jc = j * im;
UInt _jb = j * ik;
for (UInt i = 0; i < im; ++i) {
UInt _ia = i * ik;
for (UInt l = 0; l < ik; ++l) {
C[i + _jc] += A[l + _ia] * B[l + _jb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* ------------------------------------------------------------------------ */
inline void
matrix_matrixt(UInt im, UInt in, UInt ik,
Real * A, // NOLINT(readability-non-const-parameter)
Real * B, // NOLINT(readability-non-const-parameter)
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'N';
char trans_b = 'T';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &n, &beta, C,
&m);
#else
std::fill_n(C, im * in, 0.);
for (UInt j = 0; j < in; ++j) {
UInt _jc = j * im;
for (UInt i = 0; i < im; ++i) {
for (UInt l = 0; l < ik; ++l) {
UInt _la = l * im;
UInt _lb = l * in;
C[i + _jc] += A[i + _la] * B[j + _lb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* ------------------------------------------------------------------------ */
inline void
matrixt_matrixt(UInt im, UInt in, UInt ik,
Real * A, // NOLINT(readability-non-const-parameter)
Real * B, // NOLINT(readability-non-const-parameter)
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
char trans_a = 'T';
char trans_b = 'T';
double beta = 0.;
int m = im, n = in, k = ik;
aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &n, &beta, C,
&m);
#else
std::fill_n(C, im * in, 0.);
for (UInt j = 0; j < in; ++j) {
UInt _jc = j * im;
for (UInt i = 0; i < im; ++i) {
UInt _ia = i * ik;
for (UInt l = 0; l < ik; ++l) {
UInt _lb = l * in;
C[i + _jc] += A[l + _ia] * B[j + _lb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* ------------------------------------------------------------------------ */
inline void aXplusY(UInt n, Real alpha,
Real * x, // NOLINT(readability-non-const-parameter)
Real * y) {
#ifdef AKANTU_USE_BLAS
/// y := alpha x + y
int incx = 1, incy = 1;
aka_axpy(&n, &alpha, x, &incx, y, &incy);
#else
for (UInt i = 0; i < n; ++i) {
*(y++) += alpha * *(x++);
}
#endif
}
/* ------------------------------------------------------------------------ */
inline Real vectorDot(Real * v1, // NOLINT(readability-non-const-parameter)
Real * v2, // NOLINT(readability-non-const-parameter)
UInt in) {
#ifdef AKANTU_USE_BLAS
/// d := v1 . v2
int incx = 1, incy = 1, n = in;
Real d = aka_dot(&n, v1, &incx, v2, &incy);
#else
Real d = 0;
for (UInt i = 0; i < in; ++i) {
d += v1[i] * v2[i];
}
#endif
return d;
}
/* ------------------------------------------------------------------------ */
template <bool tr_A, bool tr_B>
inline void matMul(UInt m, UInt n, UInt k, Real alpha,
Real * A, // NOLINT(readability-non-const-parameter)
Real * B, // NOLINT(readability-non-const-parameter)
Real /*beta*/, Real * C) {
if (tr_A) {
if (tr_B) {
matrixt_matrixt(m, n, k, A, B, C, alpha);
} else {
matrixt_matrix(m, n, k, A, B, C, alpha);
}
} else {
if (tr_B) {
matrix_matrixt(m, n, k, A, B, C, alpha);
} else {
matrix_matrix(m, n, k, A, B, C, alpha);
}
}
}
/* ------------------------------------------------------------------------ */
template <bool tr_A>
inline void matVectMul(UInt m, UInt n, Real alpha,
Real * A, // NOLINT(readability-non-const-parameter)
Real * x, // NOLINT(readability-non-const-parameter)
Real /*beta*/, Real * y) {
if (tr_A) {
matrixt_vector(m, n, A, x, y, alpha);
} else {
matrix_vector(m, n, A, x, y, alpha);
}
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline void matrixEig(UInt n,
T * A, // NOLINT(readability-non-const-parameter)
T * d, T * V) {
// Matrix A is row major, so the lapack function in fortran will
// process A^t. Asking for the left eigenvectors of A^t will give the
// transposed right eigenvectors of A so in the C++ code the right
// eigenvectors.
char jobvr{'N'};
if (V != nullptr) {
jobvr = 'V'; // compute left eigenvectors
}
char jobvl{'N'}; // compute right eigenvectors
auto * di = new T[n]; // imaginary part of the eigenvalues
int info;
int N = n;
T wkopt;
int lwork = -1;
// query and allocate the optimal workspace
aka_geev<T>(&jobvl, &jobvr, &N, A, &N, d, di, nullptr, &N, V, &N, &wkopt,
&lwork, &info);
lwork = int(wkopt);
auto * work = new T[lwork];
// solve the eigenproblem
aka_geev<T>(&jobvl, &jobvr, &N, A, &N, d, di, nullptr, &N, V, &N, work,
&lwork, &info);
AKANTU_DEBUG_ASSERT(
info == 0,
"Problem computing eigenvalues/vectors. DGEEV exited with the value "
<< info);
delete[] work;
delete[] di; // I hope for you that there was no complex eigenvalues !!!
}
/* ------------------------------------------------------------------------ */
inline void
matrix22_eigenvalues(Real * A, // NOLINT(readability-non-const-parameter)
Real * Adiag) {
/// d = determinant of Matrix A
Real d = det2(A);
/// b = trace of Matrix A
Real b = A[0] + A[3];
Real c = std::sqrt(b * b - 4 * d);
Adiag[0] = .5 * (b + c);
Adiag[1] = .5 * (b - c);
}
/* ------------------------------------------------------------------------ */
inline void
matrix33_eigenvalues(Real * A, // NOLINT(readability-non-const-parameter)
Real * Adiag) {
matrixEig(3, A, Adiag);
}
/* ------------------------------------------------------------------------ */
template <UInt dim>
inline void eigenvalues(Real * A, // NOLINT(readability-non-const-parameter)
Real * d) {
if (dim == 1) {
d[0] = A[0];
} else if (dim == 2) {
matrix22_eigenvalues(A, d);
}
// else if(dim == 3) { matrix33_eigenvalues(A, d); }
else {
matrixEig(dim, A, d);
}
}
/* ------------------------------------------------------------------------ */
inline Real det2(const Real * mat) {
return mat[0] * mat[3] - mat[1] * mat[2];
}
/* ------------------------------------------------------------------------ */
inline Real det3(const Real * mat) {
return mat[0] * (mat[4] * mat[8] - mat[7] * mat[5]) -
mat[3] * (mat[1] * mat[8] - mat[7] * mat[2]) +
mat[6] * (mat[1] * mat[5] - mat[4] * mat[2]);
}
/* ------------------------------------------------------------------------ */
template <UInt n> inline Real det(const Real * mat) {
if (n == 1) {
return *mat;
}
if (n == 2) {
return det2(mat);
}
if (n == 3) {
return det3(mat);
}
return det(n, mat);
}
/* ------------------------------------------------------------------------ */
template <typename T> inline T det(UInt n, const T * A) {
int N = n;
int info;
auto * ipiv = new int[N + 1];
auto * LU = new T[N * N];
std::copy(A, A + N * N, LU);
// LU factorization of A
aka_getrf(&N, &N, LU, &N, ipiv, &info);
if (info > 0) {
AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info
<< " )");
}
// det(A) = det(L) * det(U) = 1 * det(U) = product_i U_{ii}
T det = 1.;
for (int i = 0; i < N; ++i) {
det *= (2 * (ipiv[i] == i) - 1) * LU[i * n + i];
}
delete[] ipiv;
delete[] LU;
return det;
}
/* ------------------------------------------------------------------------ */
inline void normal2(const Real * vec, Real * normal) {
normal[0] = vec[1];
normal[1] = -vec[0];
normalize2(normal);
}
/* ------------------------------------------------------------------------ */
inline void normal3(const Real * vec1, const Real * vec2, Real * normal) {
vectorProduct3(vec1, vec2, normal);
normalize3(normal);
}
/* ------------------------------------------------------------------------ */
inline void normalize2(Real * vec) {
Real norm = norm2(vec);
vec[0] /= norm;
vec[1] /= norm;
}
/* ------------------------------------------------------------------------ */
inline void normalize3(Real * vec) {
Real norm = norm3(vec);
vec[0] /= norm;
vec[1] /= norm;
vec[2] /= norm;
}
/* ------------------------------------------------------------------------ */
inline Real norm2(const Real * vec) {
return sqrt(vec[0] * vec[0] + vec[1] * vec[1]);
}
/* ------------------------------------------------------------------------ */
inline Real norm3(const Real * vec) {
return sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
}
/* ------------------------------------------------------------------------ */
inline Real norm(UInt n, const Real * vec) {
Real norm = 0.;
for (UInt i = 0; i < n; ++i) {
norm += vec[i] * vec[i];
}
return sqrt(norm);
}
/* ------------------------------------------------------------------------ */
inline void inv2(const Real * mat, Real * inv) {
Real det_mat = det2(mat);
inv[0] = mat[3] / det_mat;
inv[1] = -mat[1] / det_mat;
inv[2] = -mat[2] / det_mat;
inv[3] = mat[0] / det_mat;
}
/* ------------------------------------------------------------------------ */
inline void inv3(const Real * mat, Real * inv) {
Real det_mat = det3(mat);
inv[0] = (mat[4] * mat[8] - mat[7] * mat[5]) / det_mat;
inv[1] = (mat[2] * mat[7] - mat[8] * mat[1]) / det_mat;
inv[2] = (mat[1] * mat[5] - mat[4] * mat[2]) / det_mat;
inv[3] = (mat[5] * mat[6] - mat[8] * mat[3]) / det_mat;
inv[4] = (mat[0] * mat[8] - mat[6] * mat[2]) / det_mat;
inv[5] = (mat[2] * mat[3] - mat[5] * mat[0]) / det_mat;
inv[6] = (mat[3] * mat[7] - mat[6] * mat[4]) / det_mat;
inv[7] = (mat[1] * mat[6] - mat[7] * mat[0]) / det_mat;
inv[8] = (mat[0] * mat[4] - mat[3] * mat[1]) / det_mat;
}
/* ------------------------------------------------------------------------ */
template <UInt n> inline void inv(const Real * A, Real * Ainv) {
if (n == 1) {
*Ainv = 1. / *A;
} else if (n == 2) {
inv2(A, Ainv);
} else if (n == 3) {
inv3(A, Ainv);
} else {
inv(n, A, Ainv);
}
}
/* ------------------------------------------------------------------------ */
template <typename T> inline void inv(UInt n, const T * A, T * invA) {
int N = n;
int info;
auto * ipiv = new int[N + 1];
int lwork = N * N;
auto * work = new T[lwork];
std::copy(A, A + n * n, invA);
aka_getrf(&N, &N, invA, &N, ipiv, &info);
if (info > 0) {
AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info
<< " )");
}
aka_getri(&N, invA, &N, ipiv, work, &lwork, &info);
if (info != 0) {
AKANTU_ERROR("Cannot invert the matrix (info: " << info << " )");
}
delete[] ipiv;
delete[] work;
}
/* ------------------------------------------------------------------------ */
template <typename T>
inline void solve(UInt n, const T * A, T * x, const T * b) {
int N = n;
int info;
auto * ipiv = new int[N];
auto * lu_A = new T[N * N];
std::copy(A, A + N * N, lu_A);
aka_getrf(&N, &N, lu_A, &N, ipiv, &info);
if (info > 0) {
AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info
<< " )");
}
char trans = 'N';
int nrhs = 1;
std::copy(b, b + N, x);
aka_getrs(&trans, &N, &nrhs, lu_A, &N, ipiv, x, &N, &info);
if (info != 0) {
AKANTU_ERROR("Cannot solve the system (info: " << info << " )");
}
delete[] ipiv;
delete[] lu_A;
}
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
inline Real
matrixDoubleDot22(Real * A, // NOLINT(readability-non-const-parameter)
Real * B // NOLINT(readability-non-const-parameter)
) {
Real d;
d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3];
return d;
}
/* ------------------------------------------------------------------------ */
inline Real
matrixDoubleDot33(Real * A, // NOLINT(readability-non-const-parameter)
Real * B // NOLINT(readability-non-const-parameter)
) {
Real d;
d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3] + A[4] * B[4] +
A[5] * B[5] + A[6] * B[6] + A[7] * B[7] + A[8] * B[8];
return d;
}
/* ------------------------------------------------------------------------ */
inline Real
matrixDoubleDot(UInt n,
Real * A, // NOLINT(readability-non-const-parameter)
Real * B // NOLINT(readability-non-const-parameter)
) {
Real d = 0.;
for (UInt i = 0; i < n; ++i) {
for (UInt j = 0; j < n; ++j) {
d += A[i * n + j] * B[i * n + j];
}
}
return d;
}
/* ------------------------------------------------------------------------ */
inline void vectorProduct3(const Real * v1, const Real * v2, Real * res) {
res[0] = v1[1] * v2[2] - v1[2] * v2[1];
res[1] = v1[2] * v2[0] - v1[0] * v2[2];
res[2] = v1[0] * v2[1] - v1[1] * v2[0];
}
/* ------------------------------------------------------------------------ */
inline Real vectorDot2(const Real * v1, const Real * v2) {
return (v1[0] * v2[0] + v1[1] * v2[1]);
}
/* ------------------------------------------------------------------------ */
inline Real vectorDot3(const Real * v1, const Real * v2) {
return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);
}
/* ------------------------------------------------------------------------ */
inline Real distance_2d(const Real * x, const Real * y) {
return std::sqrt((y[0] - x[0]) * (y[0] - x[0]) +
(y[1] - x[1]) * (y[1] - x[1]));
}
/* ------------------------------------------------------------------------ */
inline Real triangle_inradius(const Vector<Real> & coord1, const Vector<Real> & coord2,
const Vector<Real> & coord3) {
/**
* @f{eqnarray*}{
* r &=& A / s \\
* A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)}
* \\ s &=& \frac{a + b + c}{2}
* @f}
*/
Real a, b, c;
a = coord1.distance(coord2);
b = coord2.distance(coord3);
c = coord1.distance(coord3);
Real s;
s = (a + b + c) * 0.5;
return std::sqrt((s - a) * (s - b) * (s - c) / s);
}
/* ------------------------------------------------------------------------ */
inline Real distance_3d(const Real * x, const Real * y) {
return std::sqrt((y[0] - x[0]) * (y[0] - x[0]) +
(y[1] - x[1]) * (y[1] - x[1]) +
(y[2] - x[2]) * (y[2] - x[2]));
}
/* ------------------------------------------------------------------------ */
inline Real tetrahedron_volume(const Real * coord1, const Real * coord2,
const Real * coord3, const Real * coord4) {
Real xx[9];
xx[0] = coord2[0];
xx[1] = coord2[1];
xx[2] = coord2[2];
xx[3] = coord3[0];
xx[4] = coord3[1];
xx[5] = coord3[2];
xx[6] = coord4[0];
xx[7] = coord4[1];
xx[8] = coord4[2];
auto vol = det3(xx);
xx[0] = coord1[0];
xx[1] = coord1[1];
xx[2] = coord1[2];
xx[3] = coord3[0];
xx[4] = coord3[1];
xx[5] = coord3[2];
xx[6] = coord4[0];
xx[7] = coord4[1];
xx[8] = coord4[2];
vol -= det3(xx);
xx[0] = coord1[0];
xx[1] = coord1[1];
xx[2] = coord1[2];
xx[3] = coord2[0];
xx[4] = coord2[1];
xx[5] = coord2[2];
xx[6] = coord4[0];
xx[7] = coord4[1];
xx[8] = coord4[2];
vol += det3(xx);
xx[0] = coord1[0];
xx[1] = coord1[1];
xx[2] = coord1[2];
xx[3] = coord2[0];
xx[4] = coord2[1];
xx[5] = coord2[2];
xx[6] = coord3[0];
xx[7] = coord3[1];
xx[8] = coord3[2];
vol -= det3(xx);
vol /= 6;
return vol;
}
/* ------------------------------------------------------------------------ */
inline Real tetrahedron_inradius(const Real * coord1, const Real * coord2,
const Real * coord3, const Real * coord4) {
auto l12 = distance_3d(coord1, coord2);
auto l13 = distance_3d(coord1, coord3);
auto l14 = distance_3d(coord1, coord4);
auto l23 = distance_3d(coord2, coord3);
auto l24 = distance_3d(coord2, coord4);
auto l34 = distance_3d(coord3, coord4);
auto s1 = (l12 + l23 + l13) * 0.5;
s1 = std::sqrt(s1 * (s1 - l12) * (s1 - l23) * (s1 - l13));
auto s2 = (l12 + l24 + l14) * 0.5;
s2 = std::sqrt(s2 * (s2 - l12) * (s2 - l24) * (s2 - l14));
auto s3 = (l23 + l34 + l24) * 0.5;
s3 = std::sqrt(s3 * (s3 - l23) * (s3 - l34) * (s3 - l24));
auto s4 = (l13 + l34 + l14) * 0.5;
s4 = std::sqrt(s4 * (s4 - l13) * (s4 - l34) * (s4 - l14));
auto volume = tetrahedron_volume(coord1, coord2, coord3, coord4);
return 3 * volume / (s1 + s2 + s3 + s4);
}
/* ------------------------------------------------------------------------ */
inline void barycenter(const Real * coord, UInt nb_points,
UInt spatial_dimension, Real * barycenter) {
std::fill_n(barycenter, spatial_dimension, 0.);
for (UInt n = 0; n < nb_points; ++n) {
UInt offset = n * spatial_dimension;
for (UInt i = 0; i < spatial_dimension; ++i) {
barycenter[i] += coord[offset + i] / (Real)nb_points;
}
}
}
/* ------------------------------------------------------------------------ */
inline void vector_2d(const Real * x, const Real * y, Real * res) {
res[0] = y[0] - x[0];
res[1] = y[1] - x[1];
}
/* ------------------------------------------------------------------------ */
inline void vector_3d(const Real * x, const Real * y, Real * res) {
res[0] = y[0] - x[0];
res[1] = y[1] - x[1];
res[2] = y[2] - x[2];
}
/* ------------------------------------------------------------------------ */
/// Combined absolute and relative tolerance test proposed in
/// Real-time collision detection by C. Ericson (2004)
inline bool are_float_equal(const Real x, const Real y) {
Real abs_max = std::max(std::abs(x), std::abs(y));
abs_max = std::max(abs_max, Real(1.));
return std::abs(x - y) <= (tolerance * abs_max);
}
/* ------------------------------------------------------------------------ */
inline bool isnan(Real x) {
#if defined(__INTEL_COMPILER)
#pragma warning(push)
#pragma warning(disable : 1572)
#endif // defined(__INTEL_COMPILER)
// x = x return false means x = quiet_NaN
return !(x == x);
#if defined(__INTEL_COMPILER)
#pragma warning(pop)
#endif // defined(__INTEL_COMPILER)
}
/* ------------------------------------------------------------------------ */
inline bool are_vector_equal(UInt n, Real * x, Real * y) {
bool test = true;
for (UInt i = 0; i < n; ++i) {
test &= are_float_equal(x[i], y[i]);
}
return test;
}
/* ------------------------------------------------------------------------ */
inline bool intersects(Real x_min, Real x_max, Real y_min, Real y_max) {
return not((x_max < y_min) or (x_min > y_max));
}
/* ------------------------------------------------------------------------ */
inline bool is_in_range(Real a, Real x_min, Real x_max) {
return ((a >= x_min) and (a <= x_max));
}
/* ------------------------------------------------------------------------ */
template <UInt p, typename T> inline T pow(T x) {
return (pow<p - 1, T>(x) * x);
}
template <> inline UInt pow<0, UInt>(__attribute__((unused)) UInt x) {
return (1);
}
template <> inline Real pow<0, Real>(__attribute__((unused)) Real x) {
return (1.);
}
/* ------------------------------------------------------------------------ */
template <class Functor>
Real NewtonRaphson::solve(const Functor & funct, Real x_0) {
Real x = x_0;
Real f_x = funct.f(x);
UInt iter = 0;
while (std::abs(f_x) > this->tolerance && iter < this->max_iteration) {
x -= f_x / funct.f_prime(x);
f_x = funct.f(x);
iter++;
}
AKANTU_DEBUG_ASSERT(iter < this->max_iteration,
"Newton Raphson ("
<< funct.name << ") solve did not converge in "
<< this->max_iteration << " iterations (tolerance: "
<< this->tolerance << ")");
return x;
}
} // namespace Math
} // namespace akantu
diff --git a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh
index eef8044ca..327c12940 100644
--- a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh
+++ b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh
@@ -1,116 +1,116 @@
/**
* @file element_class_segment_2_inline_impl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jul 16 2010
* @date last modification: Wed Oct 11 2017
*
* @brief Specialization of the element_class class for the type _segment_2
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*
* @verbatim
q
--x--------|--------x---> x
-1 0 1
@endverbatim
*
* @f{eqnarray*}{
* w_1(x) &=& 1/2(1 - x) \\
* w_2(x) &=& 1/2(1 + x)
* @f}
*
* @f{eqnarray*}{
* x_{q} &=& 0
* @f}
*/
/* -------------------------------------------------------------------------- */
#include "element_class.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_segment_2, _gt_segment_2,
_itp_lagrange_segment_2, _ek_regular, 1,
_git_segment, 1);
/* -------------------------------------------------------------------------- */
template <>
template <class vector_type>
inline void InterpolationElement<_itp_lagrange_segment_2>::computeShapes(
const vector_type & natural_coords, vector_type & N) {
/// natural coordinate
Real c = natural_coords(0);
/// shape functions
N(0) = 0.5 * (1 - c);
N(1) = 0.5 * (1 + c);
}
/* -------------------------------------------------------------------------- */
template <>
template <class vector_type, class matrix_type>
inline void InterpolationElement<_itp_lagrange_segment_2>::computeDNDS(
__attribute__((unused)) const vector_type & natural_coords,
matrix_type & dnds) {
/// dN1/de
dnds(0, 0) = -.5;
/// dN2/de
dnds(0, 1) = .5;
}
/* -------------------------------------------------------------------------- */
template<>
template <class vector_type, class matrix_type>
inline void InterpolationElement<_itp_lagrange_segment_2>::computeD2NDS2(
const vector_type & /*natural_coords*/, matrix_type & d2nds2) {
d2nds2.zero();
}
/* -------------------------------------------------------------------------- */
template <>
inline void
InterpolationElement<_itp_lagrange_segment_2>::computeSpecialJacobian(
const Matrix<Real> & dxds, Real & jac) {
jac = dxds.norm<L_2>();
}
/* -------------------------------------------------------------------------- */
template <>
inline Real
GeometricalElement<_gt_segment_2>::getInradius(const Matrix<Real> & coord) {
-Vector<Real> a(coord(0))
-Vector<Real> b(coord(1))
-return a.distance(b);
+ Vector<Real> a(coord(0));
+ Vector<Real> b(coord(1));
+ return a.distance(b);
}
// /* --------------------------------------------------------------------------
// */
// template<> inline bool ElementClass<_segment_2>::contains(const Vector<Real>
// & natural_coords) {
// if (natural_coords(0) < -1.) return false;
// if (natural_coords(0) > 1.) return false;
// return true;
// }
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/io/dumper/dumper_compute.hh b/src/io/dumper/dumper_compute.hh
index e84b97bb5..e43209183 100644
--- a/src/io/dumper/dumper_compute.hh
+++ b/src/io/dumper/dumper_compute.hh
@@ -1,287 +1,393 @@
/**
* @file dumper_compute.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date creation: Tue Sep 02 2014
* @date last modification: Sun Dec 03 2017
*
* @brief Field that map a function to another field
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_COMPUTE_HH_
#define AKANTU_DUMPER_COMPUTE_HH_
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "dumper_field.hh"
#include "dumper_iohelper.hh"
#include "dumper_type_traits.hh"
+/* -------------------------------------------------------------------------- */
+#include <aka_iterators.hh>
+/* -------------------------------------------------------------------------- */
#include <io_helper.hh>
-
+/* -------------------------------------------------------------------------- */
+#include <type_traits>
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
-
+ /* ------------------------------------------------------------------------ */
class ComputeFunctorInterface {
public:
virtual ~ComputeFunctorInterface() = default;
virtual UInt getDim() = 0;
virtual UInt getNbComponent(UInt old_nb_comp) = 0;
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
template <typename return_type>
class ComputeFunctorOutput : public ComputeFunctorInterface {
public:
ComputeFunctorOutput() = default;
~ComputeFunctorOutput() override = default;
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
template <typename input_type, typename return_type>
class ComputeFunctor : public ComputeFunctorOutput<return_type> {
public:
ComputeFunctor() = default;
~ComputeFunctor() override = default;
- virtual return_type func(const input_type & d, Element global_index) = 0;
+ virtual return_type func(const input_type & /*d*/, Element /*global_index*/) {
+ AKANTU_TO_IMPLEMENT();
+ }
+ virtual return_type func(const input_type & /*d*/) { AKANTU_TO_IMPLEMENT(); }
};
- /* --------------------------------------------------------------------------
- */
- template <typename SubFieldCompute, typename _return_type>
+ /* ------------------------------------------------------------------------ */
+ template <class EnumType>
+ class ComputeUIntFromEnum
+ : public ComputeFunctor<Vector<EnumType>, Vector<UInt>> {
+ public:
+ ComputeUIntFromEnum() = default;
+
+ inline Vector<UInt> func(const Vector<EnumType> & in) override {
+ Vector<UInt> out(in.size());
+ for (auto && data : zip(in, out)) {
+ std::get<1>(data) =
+ static_cast<std::underlying_type_t<EnumType>>(std::get<0>(data));
+ }
+
+ return out;
+ }
+
+ UInt getDim() override { return 1; };
+ UInt getNbComponent(UInt old_nb_comp) override { return old_nb_comp; };
+ };
+
+ /* ------------------------------------------------------------------------ */
+ template <typename SubFieldCompute, typename _return_type,
+ class support_type_ = typename SubFieldCompute::support_type>
class FieldCompute : public Field {
- /* ------------------------------------------------------------------------
- */
- /* Typedefs */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Typedefs */
+ /* ---------------------------------------------------------------------- */
public:
+ using return_type = _return_type;
+ using support_type = support_type_;
+
+ private:
using sub_iterator = typename SubFieldCompute::iterator;
using sub_types = typename SubFieldCompute::types;
using sub_return_type = typename sub_types::return_type;
+ using data_type = typename sub_types::data_type;
+ using functor_type = ComputeFunctor<sub_return_type, return_type>;
+
+ public:
+ class iterator {
+ public:
+ iterator(const sub_iterator & it, functor_type & func)
+ : it(it), func(func) {}
+
+ bool operator!=(const iterator & it) const { return it.it != this->it; }
+
+ iterator operator++() {
+ ++this->it;
+ return *this;
+ }
+
+ return_type operator*() { return func.func(*it); }
+
+ protected:
+ sub_iterator it;
+ functor_type & func;
+ };
+
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
+ public:
+ FieldCompute(SubFieldCompute & cont,
+ std::unique_ptr<ComputeFunctorInterface> func)
+ : sub_field(aka::as_type<SubFieldCompute>(cont.shared_from_this())),
+ func(aka::as_type<functor_type>(func.release())) {
+ this->checkHomogeneity();
+ };
+
+ void registerToDumper(const std::string & id,
+ iohelper::Dumper & dumper) override {
+ dumper.addNodeDataField(id, *this);
+ }
+
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
+ public:
+ iterator begin() { return iterator(sub_field->begin(), *func); }
+ iterator end() { return iterator(sub_field->end(), *func); }
+
+ UInt getDim() { return func->getDim(); }
+
+ UInt size() {
+ throw;
+ // return Functor::size();
+ return 0;
+ }
+
+ void checkHomogeneity() override { this->homogeneous = true; };
+
+ iohelper::DataType getDataType() {
+ return iohelper::getDataType<data_type>();
+ }
+
+ /// for connection to a FieldCompute
+ inline std::shared_ptr<Field> connect(FieldComputeProxy & proxy) override;
+
+ /// for connection to a FieldCompute
+ std::unique_ptr<ComputeFunctorInterface>
+ connect(HomogenizerProxy & proxy) override;
+
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
+ public:
+ std::shared_ptr<SubFieldCompute> sub_field;
+ std::unique_ptr<functor_type> func;
+ };
+
+ /* ------------------------------------------------------------------------ */
+ template <typename SubFieldCompute, typename _return_type>
+ class FieldCompute<SubFieldCompute, _return_type, Element> : public Field {
+ /* ---------------------------------------------------------------------- */
+ /* Typedefs */
+ /* ---------------------------------------------------------------------- */
+ public:
using return_type = _return_type;
+ using support_type = Element;
+
+ using sub_iterator = typename SubFieldCompute::iterator;
+ using sub_types = typename SubFieldCompute::types;
+ using sub_return_type = typename sub_types::return_type;
using data_type = typename sub_types::data_type;
+ using functor_type = ComputeFunctor<sub_return_type, return_type>;
using types =
TypeTraits<data_type, return_type, ElementTypeMapArray<data_type>>;
+ public:
class iterator {
public:
- iterator(const sub_iterator & it,
- ComputeFunctor<sub_return_type, return_type> & func)
+ iterator(const sub_iterator & it, functor_type & func)
: it(it), func(func) {}
bool operator!=(const iterator & it) const { return it.it != this->it; }
iterator operator++() {
++this->it;
return *this;
}
UInt currentGlobalIndex() { return this->it.currentGlobalIndex(); }
return_type operator*() { return func.func(*it, it.getCurrentElement()); }
Element getCurrentElement() { return this->it.getCurrentElement(); }
UInt element_type() { return this->it.element_type(); }
protected:
sub_iterator it;
- ComputeFunctor<sub_return_type, return_type> & func;
+ functor_type & func;
};
- /* ------------------------------------------------------------------------
- */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
public:
FieldCompute(SubFieldCompute & cont,
std::unique_ptr<ComputeFunctorInterface> func)
: sub_field(aka::as_type<SubFieldCompute>(cont.shared_from_this())),
- func(aka::as_type<ComputeFunctor<sub_return_type, return_type>>(
- func.release())) {
+ func(aka::as_type<functor_type>(func.release())) {
this->checkHomogeneity();
};
~FieldCompute() override = default;
void registerToDumper(const std::string & id,
iohelper::Dumper & dumper) override {
dumper.addElemDataField(id, *this);
}
- /* ------------------------------------------------------------------------
- */
- /* Class Members */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
public:
iterator begin() { return iterator(sub_field->begin(), *func); }
iterator end() { return iterator(sub_field->end(), *func); }
UInt getDim() { return func->getDim(); }
UInt size() {
throw;
// return Functor::size();
return 0;
}
void checkHomogeneity() override { this->homogeneous = true; };
iohelper::DataType getDataType() {
return iohelper::getDataType<data_type>();
}
/// get the number of components of the hosted field
ElementTypeMap<UInt>
getNbComponents(UInt dim = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_not_defined) override {
ElementTypeMap<UInt> nb_components;
const auto & old_nb_components =
this->sub_field->getNbComponents(dim, ghost_type, kind);
for (auto type : old_nb_components.elementTypes(dim, ghost_type, kind)) {
UInt nb_comp = old_nb_components(type, ghost_type);
nb_components(type, ghost_type) = func->getNbComponent(nb_comp);
}
return nb_components;
};
/// for connection to a FieldCompute
inline std::shared_ptr<Field> connect(FieldComputeProxy & proxy) override;
/// for connection to a FieldCompute
std::unique_ptr<ComputeFunctorInterface>
connect(HomogenizerProxy & proxy) override;
- /* ------------------------------------------------------------------------
- */
- /* Class Members */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
public:
std::shared_ptr<SubFieldCompute> sub_field;
- std::unique_ptr<ComputeFunctor<sub_return_type, return_type>> func;
+ std::unique_ptr<functor_type> func;
};
- /* --------------------------------------------------------------------------
- */
-
- /* --------------------------------------------------------------------------
- */
-
+ /* ------------------------------------------------------------------------ */
class FieldComputeProxy {
- /* ------------------------------------------------------------------------
- */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
public:
FieldComputeProxy(std::unique_ptr<ComputeFunctorInterface> func)
: func(std::move(func)){};
inline static std::shared_ptr<Field>
createFieldCompute(std::shared_ptr<Field> & field,
std::unique_ptr<ComputeFunctorInterface> func) {
FieldComputeProxy compute_proxy(std::move(func));
return field->connect(compute_proxy);
}
template <typename T> std::shared_ptr<Field> connectToField(T * ptr) {
if (aka::is_of_type<ComputeFunctorOutput<Vector<Real>>>(func)) {
return this->connectToFunctor<Vector<Real>>(ptr);
}
if (aka::is_of_type<ComputeFunctorOutput<Vector<UInt>>>(func)) {
return this->connectToFunctor<Vector<UInt>>(ptr);
}
if (aka::is_of_type<ComputeFunctorOutput<Matrix<UInt>>>(func)) {
return this->connectToFunctor<Matrix<UInt>>(ptr);
}
if (aka::is_of_type<ComputeFunctorOutput<Matrix<Real>>>(func)) {
return this->connectToFunctor<Matrix<Real>>(ptr);
}
throw;
}
template <typename output, typename T>
std::shared_ptr<Field> connectToFunctor(T * ptr) {
return std::make_shared<FieldCompute<T, output>>(*ptr, std::move(func));
}
template <typename output, typename SubFieldCompute, typename return_type1,
typename return_type2>
std::shared_ptr<Field>
connectToFunctor(FieldCompute<FieldCompute<SubFieldCompute, return_type1>,
return_type2> * /*ptr*/) {
throw; // return new FieldCompute<T,output>(*ptr,func);
return nullptr;
}
template <typename output, typename SubFieldCompute, typename return_type1,
typename return_type2, typename return_type3,
typename return_type4>
std::shared_ptr<Field> connectToFunctor(
FieldCompute<FieldCompute<FieldCompute<FieldCompute<SubFieldCompute,
return_type1>,
return_type2>,
return_type3>,
return_type4> * /*ptr*/) {
throw; // return new FieldCompute<T,output>(*ptr,func);
return nullptr;
}
- /* ------------------------------------------------------------------------
- */
- /* Class Members */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
public:
std::unique_ptr<ComputeFunctorInterface> func;
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
/// for connection to a FieldCompute
+ template <typename SubFieldCompute, typename return_type,
+ typename support_type_>
+ inline std::shared_ptr<Field>
+ FieldCompute<SubFieldCompute, return_type, support_type_>::connect(
+ FieldComputeProxy & proxy) {
+ return proxy.connectToField(this);
+ }
+
template <typename SubFieldCompute, typename return_type>
inline std::shared_ptr<Field>
- FieldCompute<SubFieldCompute, return_type>::connect(
+ FieldCompute<SubFieldCompute, return_type, Element>::connect(
FieldComputeProxy & proxy) {
return proxy.connectToField(this);
}
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
} // namespace dumpers
} // namespace akantu
#endif /* AKANTU_DUMPER_COMPUTE_HH_ */
diff --git a/src/io/dumper/dumper_elemental_field.hh b/src/io/dumper/dumper_elemental_field.hh
index 409bad555..7a8e739c4 100644
--- a/src/io/dumper/dumper_elemental_field.hh
+++ b/src/io/dumper/dumper_elemental_field.hh
@@ -1,75 +1,75 @@
/**
* @file dumper_elemental_field.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Sun Dec 03 2017
*
* @brief description of elemental fields
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_ELEMENTAL_FIELD_HH_
#define AKANTU_DUMPER_ELEMENTAL_FIELD_HH_
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
#include "dumper_field.hh"
#include "dumper_generic_elemental_field.hh"
#ifdef AKANTU_IGFEM
#include "dumper_igfem_elemental_field.hh"
#endif
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
/* -------------------------------------------------------------------------- */
template <typename T, template <class> class ret = Vector,
bool filtered = false>
class ElementalField
: public GenericElementalField<SingleType<T, ret, filtered>,
elemental_field_iterator> {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
using types = SingleType<T, ret, filtered>;
using field_type = typename types::field_type;
using iterator = elemental_field_iterator<types>;
-
+ using support_type = Element;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ElementalField(const field_type & field,
UInt spatial_dimension = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined)
: GenericElementalField<types, elemental_field_iterator>(
field, spatial_dimension, ghost_type, element_kind) {}
};
/* -------------------------------------------------------------------------- */
} // namespace dumpers
} // namespace akantu
#endif /* AKANTU_DUMPER_ELEMENTAL_FIELD_HH_ */
diff --git a/src/io/dumper/dumper_generic_elemental_field.hh b/src/io/dumper/dumper_generic_elemental_field.hh
index 26fbf377e..7a2a2ffa3 100644
--- a/src/io/dumper/dumper_generic_elemental_field.hh
+++ b/src/io/dumper/dumper_generic_elemental_field.hh
@@ -1,220 +1,221 @@
/**
* @file dumper_generic_elemental_field.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Sep 02 2014
* @date last modification: Wed Nov 08 2017
*
* @brief Generic interface for elemental fields
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH_
#define AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH_
/* -------------------------------------------------------------------------- */
#include "dumper_element_iterator.hh"
#include "dumper_field.hh"
#include "dumper_homogenizing_field.hh"
#include "element_type_map_filter.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
/* -------------------------------------------------------------------------- */
template <class _types, template <class> class iterator_type>
class GenericElementalField : public Field {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
// check dumper_type_traits.hh for additional information over these types
using types = _types;
using data_type = typename types::data_type;
using it_type = typename types::it_type;
using field_type = typename types::field_type;
using array_type = typename types::array_type;
using array_iterator = typename types::array_iterator;
using field_type_iterator = typename field_type::type_iterator;
using iterator = iterator_type<types>;
+ using support_type = Element;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
GenericElementalField(const field_type & field,
UInt spatial_dimension = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined)
: field(field), spatial_dimension(spatial_dimension),
ghost_type(ghost_type), element_kind(element_kind) {
this->checkHomogeneity();
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// get the number of components of the hosted field
ElementTypeMap<UInt>
getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_not_defined) override {
return this->field.getNbComponents(dim, ghost_type, kind);
};
/// return the size of the contained data: i.e. the number of elements ?
virtual UInt size() {
checkHomogeneity();
return this->nb_total_element;
}
/// return the iohelper datatype to be dumped
iohelper::DataType getDataType() {
return iohelper::getDataType<data_type>();
}
protected:
/// return the number of entries per element
UInt getNbDataPerElem(ElementType type,
GhostType ghost_type = _not_ghost) const {
if (!nb_data_per_elem.exists(type, ghost_type)) {
return field(type, ghost_type).getNbComponent();
}
return nb_data_per_elem(type, this->ghost_type);
}
/// check if the same quantity of data for all element types
void checkHomogeneity() override;
public:
void registerToDumper(const std::string & id,
iohelper::Dumper & dumper) override {
dumper.addElemDataField(id, *this);
};
/// for connection to a FieldCompute
inline std::shared_ptr<Field> connect(FieldComputeProxy & proxy) override {
return proxy.connectToField(this);
}
/// for connection to a Homogenizer
inline std::unique_ptr<ComputeFunctorInterface>
connect(HomogenizerProxy & proxy) override {
return proxy.connectToField(this);
};
virtual iterator begin() {
/// type iterators on the elemental field
auto types = this->field.elementTypes(this->spatial_dimension,
this->ghost_type, this->element_kind);
auto tit = types.begin();
auto end = types.end();
/// skip all types without data
for (; tit != end and this->field(*tit, this->ghost_type).empty();
++tit) {
}
auto type = *tit;
if (tit == end) {
return this->end();
}
/// getting information for the field of the given type
const auto & vect = this->field(type, this->ghost_type);
UInt nb_data_per_elem = this->getNbDataPerElem(type);
/// define element-wise iterator
auto view = make_view(vect, nb_data_per_elem);
auto it = view.begin();
auto it_end = view.end();
/// define data iterator
iterator rit =
iterator(this->field, tit, end, it, it_end, this->ghost_type);
rit.setNbDataPerElem(this->nb_data_per_elem);
return rit;
}
virtual iterator end() {
auto types = this->field.elementTypes(this->spatial_dimension,
this->ghost_type, this->element_kind);
auto tit = types.begin();
auto end = types.end();
auto type = *tit;
for (; tit != end; ++tit) {
type = *tit;
}
const array_type & vect = this->field(type, this->ghost_type);
UInt nb_data = this->getNbDataPerElem(type);
auto it = make_view(vect, nb_data).end();
auto rit = iterator(this->field, end, end, it, it, this->ghost_type);
rit.setNbDataPerElem(this->nb_data_per_elem);
return rit;
}
virtual UInt getDim() {
if (this->homogeneous) {
auto tit = this->field
.elementTypes(this->spatial_dimension, this->ghost_type,
this->element_kind)
.begin();
return this->getNbDataPerElem(*tit);
}
throw;
return 0;
}
void setNbDataPerElem(const ElementTypeMap<UInt> & nb_data) override {
nb_data_per_elem = nb_data;
}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// the ElementTypeMapArray embedded in the field
const field_type & field;
/// total number of elements
UInt nb_total_element;
/// the spatial dimension of the problem
UInt spatial_dimension;
/// whether this is a ghost field or not (for type selection)
GhostType ghost_type;
/// The element kind to operate on
ElementKind element_kind;
/// The number of data per element type
ElementTypeMap<UInt> nb_data_per_elem;
};
} // namespace dumpers
} // namespace akantu
/* -------------------------------------------------------------------------- */
#include "dumper_generic_elemental_field_tmpl.hh"
/* -------------------------------------------------------------------------- */
#endif /* AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH_ */
diff --git a/src/io/dumper/dumper_homogenizing_field.hh b/src/io/dumper/dumper_homogenizing_field.hh
index 96395e2b1..fe8b52990 100644
--- a/src/io/dumper/dumper_homogenizing_field.hh
+++ b/src/io/dumper/dumper_homogenizing_field.hh
@@ -1,204 +1,192 @@
/**
* @file dumper_homogenizing_field.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Sep 02 2014
* @date last modification: Wed Nov 08 2017
*
* @brief description of field homogenizing field
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_HOMOGENIZING_FIELD_HH_
#define AKANTU_DUMPER_HOMOGENIZING_FIELD_HH_
/* -------------------------------------------------------------------------- */
#include "dumper_compute.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
-/* -------------------------------------------------------------------------- */
-
-template <typename type>
-inline type
-typeConverter(const type & input,
- [[gnu::unused]] Vector<typename type::value_type> & res,
- [[gnu::unused]] UInt nb_data) {
- throw;
- return input;
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename type>
-inline Matrix<type> typeConverter(const Matrix<type> & input,
- Vector<type> & res, UInt nb_data) {
-
- Matrix<type> tmp(res.storage(), input.rows(), nb_data / input.rows());
- Matrix<type> tmp2(tmp, true);
- return tmp2;
-}
-
-/* -------------------------------------------------------------------------- */
+ /* ------------------------------------------------------------------------ */
+ template <typename type>
+ inline type
+ typeConverter(const type & input,
+ [[gnu::unused]] Vector<typename type::value_type> & res,
+ [[gnu::unused]] UInt nb_data) {
+ throw;
+ return input;
+ }
-template <typename type>
-inline Vector<type> typeConverter(const Vector<type> & /*unused*/,
- Vector<type> & res, UInt /*unused*/) {
- return res;
-}
+ /* ------------------------------------------------------------------------ */
+ template <typename type>
+ inline Matrix<type> typeConverter(const Matrix<type> & input,
+ Vector<type> & res, UInt nb_data) {
-/* -------------------------------------------------------------------------- */
+ Matrix<type> tmp(res.storage(), input.rows(), nb_data / input.rows());
+ Matrix<type> tmp2(tmp, true);
+ return tmp2;
+ }
-template <typename type>
-class AvgHomogenizingFunctor : public ComputeFunctor<type, type> {
- /* ------------------------------------------------------------------------ */
- /* Typedefs */
/* ------------------------------------------------------------------------ */
- using value_type = typename type::value_type;
+ template <typename type>
+ inline Vector<type> typeConverter(const Vector<type> & /*unused*/,
+ Vector<type> & res, UInt /*unused*/) {
+ return res;
+ }
/* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------ */
-public:
- AvgHomogenizingFunctor(ElementTypeMap<UInt> & nb_datas) {
-
- auto types = nb_datas.elementTypes();
- auto tit = types.begin();
- auto end = types.end();
+ template <typename type>
+ class AvgHomogenizingFunctor : public ComputeFunctor<type, type> {
+ /* ---------------------------------------------------------------------- */
+ /* Typedefs */
+ /* ---------------------------------------------------------------------- */
+ private:
+ using value_type = typename type::value_type;
+
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
+ public:
+ AvgHomogenizingFunctor(ElementTypeMap<UInt> & nb_datas) {
+
+ auto types = nb_datas.elementTypes();
+ auto tit = types.begin();
+ auto end = types.end();
+
+ nb_data = nb_datas(*tit);
+
+ for (; tit != end; ++tit) {
+ if (nb_data != nb_datas(*tit)) {
+ throw;
+ }
+ }
+ }
- nb_data = nb_datas(*tit);
+ /* ---------------------------------------------------------------------- */
+ /* Methods */
+ /* ---------------------------------------------------------------------- */
+ public:
+ type func(const type & d, Element /*global_index*/) override {
+ Vector<value_type> res(this->nb_data);
- for (; tit != end; ++tit) {
- if (nb_data != nb_datas(*tit)) {
+ if (d.size() % this->nb_data) {
throw;
}
- }
- }
+ UInt nb_to_average = d.size() / this->nb_data;
- /* ------------------------------------------------------------------------ */
- /* Methods */
- /* ------------------------------------------------------------------------ */
+ value_type * ptr = d.storage();
+ for (UInt i = 0; i < nb_to_average; ++i) {
+ Vector<value_type> tmp(ptr, this->nb_data);
+ res += tmp;
+ ptr += this->nb_data;
+ }
+ res /= nb_to_average;
+ return typeConverter(d, res, this->nb_data);
+ };
- type func(const type & d, Element /*global_index*/) override {
- Vector<value_type> res(this->nb_data);
+ UInt getDim() override { return nb_data; };
+ UInt getNbComponent(UInt /*old_nb_comp*/) override { throw; };
- if (d.size() % this->nb_data) {
- throw;
- }
- UInt nb_to_average = d.size() / this->nb_data;
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
- value_type * ptr = d.storage();
- for (UInt i = 0; i < nb_to_average; ++i) {
- Vector<value_type> tmp(ptr, this->nb_data);
- res += tmp;
- ptr += this->nb_data;
- }
- res /= nb_to_average;
- return typeConverter(d, res, this->nb_data);
+ /// The size of data: i.e. the size of the vector to be returned
+ UInt nb_data;
};
-
- UInt getDim() override { return nb_data; };
- UInt getNbComponent(UInt /*old_nb_comp*/) override { throw; };
-
- /* ------------------------------------------------------------------------ */
- /* Class Members */
/* ------------------------------------------------------------------------ */
- /// The size of data: i.e. the size of the vector to be returned
- UInt nb_data;
-};
-/* -------------------------------------------------------------------------- */
-
-/* -------------------------------------------------------------------------- */
-
-class HomogenizerProxy {
- /* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
-public:
- HomogenizerProxy() = default;
-
-public:
- inline static std::unique_ptr<ComputeFunctorInterface>
- createHomogenizer(Field & field);
+ class HomogenizerProxy {
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
+ public:
+ HomogenizerProxy() = default;
+
+ public:
+ inline static std::unique_ptr<ComputeFunctorInterface>
+ createHomogenizer(Field & field);
+
+ template <typename T>
+ inline std::unique_ptr<ComputeFunctorInterface> connectToField(T * field) {
+ ElementTypeMap<UInt> nb_components = field->getNbComponents();
+
+ using ret_type = typename T::types::return_type;
+ return this->instantiateHomogenizer<ret_type>(nb_components);
+ }
- template <typename T>
- inline std::unique_ptr<ComputeFunctorInterface> connectToField(T * field) {
- ElementTypeMap<UInt> nb_components = field->getNbComponents();
+ template <typename ret_type>
+ inline std::unique_ptr<ComputeFunctorInterface>
+ instantiateHomogenizer(ElementTypeMap<UInt> & nb_components);
+ };
- using ret_type = typename T::types::return_type;
- return this->instantiateHomogenizer<ret_type>(nb_components);
- }
+ /* ------------------------------------------------------------------------ */
template <typename ret_type>
inline std::unique_ptr<ComputeFunctorInterface>
- instantiateHomogenizer(ElementTypeMap<UInt> & nb_components);
-};
-
-/* -------------------------------------------------------------------------- */
-
-template <typename ret_type>
-inline std::unique_ptr<ComputeFunctorInterface>
-HomogenizerProxy::instantiateHomogenizer(ElementTypeMap<UInt> & nb_components) {
- using Homogenizer = dumpers::AvgHomogenizingFunctor<ret_type>;
- return std::make_unique<Homogenizer>(nb_components);
-}
-
-template <>
-inline std::unique_ptr<ComputeFunctorInterface>
-HomogenizerProxy::instantiateHomogenizer<Vector<iohelper::ElemType>>([
- [gnu::unused]] ElementTypeMap<UInt> & nb_components) {
- throw;
- return nullptr;
-}
+ HomogenizerProxy::instantiateHomogenizer(
+ ElementTypeMap<UInt> & nb_components) {
+ using Homogenizer = dumpers::AvgHomogenizingFunctor<ret_type>;
+ return std::make_unique<Homogenizer>(nb_components);
+ }
-/* -------------------------------------------------------------------------- */
-/// for connection to a FieldCompute
-template <typename SubFieldCompute, typename return_type>
-inline std::unique_ptr<ComputeFunctorInterface>
-FieldCompute<SubFieldCompute, return_type>::connect(HomogenizerProxy & proxy) {
- return proxy.connectToField(this);
-}
+ template <>
+ inline std::unique_ptr<ComputeFunctorInterface>
+ HomogenizerProxy::instantiateHomogenizer<Vector<iohelper::ElemType>>(
+ [[gnu::unused]] ElementTypeMap<UInt> & nb_components) {
+ throw;
+ return nullptr;
+ }
-/* -------------------------------------------------------------------------- */
-inline std::unique_ptr<ComputeFunctorInterface>
-HomogenizerProxy::createHomogenizer(Field & field) {
- HomogenizerProxy homogenizer_proxy;
- return field.connect(homogenizer_proxy);
-}
+ /* ------------------------------------------------------------------------ */
+ /// for connection to a FieldCompute
+ template <typename SubFieldCompute, typename return_type,
+ typename support_type_>
+ inline std::unique_ptr<ComputeFunctorInterface>
+ FieldCompute<SubFieldCompute, return_type, support_type_>::connect(
+ HomogenizerProxy & proxy) {
+ return proxy.connectToField(this);
+ }
-/* -------------------------------------------------------------------------- */
-// inline ComputeFunctorInterface & createHomogenizer(Field & field){
-// HomogenizerProxy::createHomogenizer(field);
-// throw;
-// ComputeFunctorInterface * ptr = NULL;
-// return *ptr;
-// }
-
-// /* --------------------------------------------------------------------------
-// */
+ /* ------------------------------------------------------------------------ */
+ inline std::unique_ptr<ComputeFunctorInterface>
+ HomogenizerProxy::createHomogenizer(Field & field) {
+ HomogenizerProxy homogenizer_proxy;
+ return field.connect(homogenizer_proxy);
+ }
} // namespace dumpers
} // namespace akantu
#endif /* AKANTU_DUMPER_HOMOGENIZING_FIELD_HH_ */
diff --git a/src/io/dumper/dumper_internal_material_field.hh b/src/io/dumper/dumper_internal_material_field.hh
index 8bad2ffef..ba233717e 100644
--- a/src/io/dumper/dumper_internal_material_field.hh
+++ b/src/io/dumper/dumper_internal_material_field.hh
@@ -1,71 +1,72 @@
/**
* @file dumper_internal_material_field.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Nov 08 2017
*
* @brief description of material internal field
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_INTERNAL_MATERIAL_FIELD_HH_
#define AKANTU_DUMPER_INTERNAL_MATERIAL_FIELD_HH_
/* -------------------------------------------------------------------------- */
#include "dumper_quadrature_point_iterator.hh"
#ifdef AKANTU_IGFEM
#include "dumper_igfem_material_internal_field.hh"
#endif
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
/* -------------------------------------------------------------------------- */
template <typename T, bool filtered = false>
class InternalMaterialField
: public GenericElementalField<SingleType<T, Vector, filtered>,
quadrature_point_iterator> {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
using types = SingleType<T, Vector, filtered>;
using parent = GenericElementalField<types, quadrature_point_iterator>;
using field_type = typename types::field_type;
+ using support_type = Element;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
InternalMaterialField(const field_type & field,
UInt spatial_dimension = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined)
: parent(field, spatial_dimension, ghost_type, element_kind) {}
};
} // namespace dumpers
} // namespace akantu
#endif /* AKANTU_DUMPER_INTERNAL_MATERIAL_FIELD_HH_ */
diff --git a/src/io/dumper/dumper_material_padders.hh b/src/io/dumper/dumper_material_padders.hh
index 7ee640f0d..d3b02a5f0 100644
--- a/src/io/dumper/dumper_material_padders.hh
+++ b/src/io/dumper/dumper_material_padders.hh
@@ -1,319 +1,305 @@
/**
* @file dumper_material_padders.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Sep 02 2014
* @date last modification: Wed Nov 29 2017
*
* @brief Material padders for plane stress/ plane strain
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_MATERIAL_PADDERS_HH_
#define AKANTU_DUMPER_MATERIAL_PADDERS_HH_
/* -------------------------------------------------------------------------- */
#include "dumper_padding_helper.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
class MaterialFunctor {
- /* ------------------------------------------------------------------------
- */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
public:
MaterialFunctor(const SolidMechanicsModel & model)
: model(model), material_index(model.getMaterialByElement()),
nb_data_per_element("nb_data_per_element", model.getID()),
spatial_dimension(model.getSpatialDimension()) {}
- /* ------------------------------------------------------------------------
- */
- /* Methods */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Methods */
+ /* ---------------------------------------------------------------------- */
/// return the material from the global element index
const Material & getMaterialFromGlobalIndex(Element global_index) {
UInt index = global_index.element;
UInt material_id = material_index(global_index.type)(index);
const Material & material = model.getMaterial(material_id);
return material;
}
/// return the type of the element from global index
ElementType
getElementTypeFromGlobalIndex( // NOLINT(readability-convert-member-functions-to-static)
Element global_index) {
return global_index.type;
}
protected:
- /* ------------------------------------------------------------------------
- */
- /* Class Members */
- /* ------------------------------------------------------------------------
- */
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
/// all material padders probably need access to solid mechanics model
const SolidMechanicsModel & model;
/// they also need an access to the map from global ids to material id and
/// local ids
const ElementTypeMapArray<UInt> & material_index;
/// the number of data per element
const ElementTypeMapArray<UInt> nb_data_per_element;
UInt spatial_dimension;
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
template <class T, class R>
class MaterialPadder : public MaterialFunctor,
public PadderGeneric<Vector<T>, R> {
public:
MaterialPadder(const SolidMechanicsModel & model)
: MaterialFunctor(model) {}
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
template <UInt spatial_dimension>
class StressPadder : public MaterialPadder<Real, Matrix<Real>> {
public:
StressPadder(const SolidMechanicsModel & model)
: MaterialPadder<Real, Matrix<Real>>(model) {
this->setPadding(3, 3);
}
inline Matrix<Real> func(const Vector<Real> & in,
Element global_element_id) override {
UInt nrows = spatial_dimension;
UInt ncols = in.size() / nrows;
UInt nb_data = in.size() / (nrows * nrows);
Matrix<Real> stress = this->pad(in, nrows, ncols, nb_data);
const Material & material =
this->getMaterialFromGlobalIndex(global_element_id);
bool plane_strain = true;
if (spatial_dimension == 2) {
plane_strain = !((bool)material.getParam("Plane_Stress"));
}
if (plane_strain) {
Real nu = material.getParam("nu");
for (UInt d = 0; d < nb_data; ++d) {
stress(2, 2 + 3 * d) =
nu * (stress(0, 0 + 3 * d) + stress(1, 1 + 3 * d));
}
}
return stress;
}
UInt getDim() override { return 9; };
UInt getNbComponent(UInt /*old_nb_comp*/) override {
return this->getDim();
};
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
template <UInt spatial_dimension>
class StrainPadder : public MaterialFunctor,
public PadderGeneric<Matrix<Real>, Matrix<Real>> {
public:
StrainPadder(const SolidMechanicsModel & model) : MaterialFunctor(model) {
this->setPadding(3, 3);
}
inline Matrix<Real> func(const Matrix<Real> & in,
Element global_element_id) override {
UInt nrows = spatial_dimension;
UInt nb_data = in.size() / (nrows * nrows);
Matrix<Real> strain = this->pad(in, nb_data);
const Material & material =
this->getMaterialFromGlobalIndex(global_element_id);
bool plane_stress = material.getParam("Plane_Stress");
if (plane_stress) {
Real nu = material.getParam("nu");
for (UInt d = 0; d < nb_data; ++d) {
strain(2, 2 + 3 * d) =
nu / (nu - 1) * (strain(0, 0 + 3 * d) + strain(1, 1 + 3 * d));
}
}
return strain;
}
UInt getDim() override { return 9; };
UInt getNbComponent(UInt /*old_nb_comp*/) override {
return this->getDim();
};
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
template <bool green_strain>
class ComputeStrain : public MaterialFunctor,
public ComputeFunctor<Vector<Real>, Matrix<Real>> {
public:
ComputeStrain(const SolidMechanicsModel & model) : MaterialFunctor(model) {}
inline Matrix<Real> func(const Vector<Real> & in,
Element /*global_element_id*/) override {
UInt nrows = spatial_dimension;
UInt ncols = in.size() / nrows;
UInt nb_data = in.size() / (nrows * nrows);
Matrix<Real> ret_all_strain(nrows, ncols);
Tensor3<Real> all_grad_u(in.storage(), nrows, nrows, nb_data);
Tensor3<Real> all_strain(ret_all_strain.storage(), nrows, nrows, nb_data);
for (UInt d = 0; d < nb_data; ++d) {
Matrix<Real> grad_u = all_grad_u(d);
Matrix<Real> strain = all_strain(d);
if (spatial_dimension == 2) {
if (green_strain) {
Material::gradUToE<2>(grad_u, strain);
} else {
Material::gradUToEpsilon<2>(grad_u, strain);
}
} else if (spatial_dimension == 3) {
if (green_strain) {
Material::gradUToE<3>(grad_u, strain);
} else {
Material::gradUToEpsilon<3>(grad_u, strain);
}
}
}
return ret_all_strain;
}
UInt getDim() override { return spatial_dimension * spatial_dimension; };
UInt getNbComponent(UInt /*old_nb_comp*/) override {
return this->getDim();
};
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
template <bool green_strain>
class ComputePrincipalStrain
: public MaterialFunctor,
public ComputeFunctor<Vector<Real>, Matrix<Real>> {
public:
ComputePrincipalStrain(const SolidMechanicsModel & model)
: MaterialFunctor(model) {}
inline Matrix<Real> func(const Vector<Real> & in,
Element /*global_element_id*/) override {
UInt nrows = spatial_dimension;
UInt nb_data = in.size() / (nrows * nrows);
Matrix<Real> ret_all_strain(nrows, nb_data);
Tensor3<Real> all_grad_u(in.storage(), nrows, nrows, nb_data);
Matrix<Real> strain(nrows, nrows);
for (UInt d = 0; d < nb_data; ++d) {
Matrix<Real> grad_u = all_grad_u(d);
if (spatial_dimension == 2) {
if (green_strain) {
Material::gradUToE<2>(grad_u, strain);
} else {
Material::gradUToEpsilon<2>(grad_u, strain);
}
} else if (spatial_dimension == 3) {
if (green_strain) {
Material::gradUToE<3>(grad_u, strain);
} else {
Material::gradUToEpsilon<3>(grad_u, strain);
}
}
Vector<Real> principal_strain(ret_all_strain(d));
strain.eig(principal_strain);
}
return ret_all_strain;
}
UInt getDim() override { return spatial_dimension; };
UInt getNbComponent(UInt /*old_nb_comp*/) override {
return this->getDim();
};
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
class ComputeVonMisesStress
: public MaterialFunctor,
public ComputeFunctor<Vector<Real>, Vector<Real>> {
public:
ComputeVonMisesStress(const SolidMechanicsModel & model)
: MaterialFunctor(model) {}
inline Vector<Real> func(const Vector<Real> & in,
Element /*global_element_id*/) override {
UInt nrows = spatial_dimension;
UInt nb_data = in.size() / (nrows * nrows);
Vector<Real> von_mises_stress(nb_data);
Matrix<Real> deviatoric_stress(3, 3);
for (UInt d = 0; d < nb_data; ++d) {
Matrix<Real> cauchy_stress(in.storage() + d * nrows * nrows, nrows,
nrows);
von_mises_stress(d) = Material::stressToVonMises(cauchy_stress);
}
return von_mises_stress;
}
UInt getDim() override { return 1; };
UInt getNbComponent(UInt /*old_nb_comp*/) override {
return this->getDim();
};
};
- /* --------------------------------------------------------------------------
- */
+ /* ------------------------------------------------------------------------ */
} // namespace dumpers
} // namespace akantu
#endif /* AKANTU_DUMPER_MATERIAL_PADDERS_HH_ */
diff --git a/src/io/dumper/dumper_nodal_field.hh b/src/io/dumper/dumper_nodal_field.hh
index 2749754a9..e8e6b4386 100644
--- a/src/io/dumper/dumper_nodal_field.hh
+++ b/src/io/dumper/dumper_nodal_field.hh
@@ -1,243 +1,238 @@
/**
* @file dumper_nodal_field.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Oct 26 2012
* @date last modification: Wed Nov 08 2017
*
* @brief Description of nodal fields
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_NODAL_FIELD_HH_
#define AKANTU_DUMPER_NODAL_FIELD_HH_
#include "dumper_field.hh"
#include <io_helper.hh>
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
-// This represents a iohelper compatible field
-template <typename T, bool filtered = false, class Container = Array<T>,
- class Filter = Array<UInt>>
-class NodalField;
+ // This represents a iohelper compatible field
+ template <typename T, bool filtered = false, class Container = Array<T>,
+ class Filter = Array<UInt>>
+ class NodalField;
-/* -------------------------------------------------------------------------- */
-template <typename T, class Container, class Filter>
-class NodalField<T, false, Container, Filter> : public dumpers::Field {
-public:
- /* ------------------------------------------------------------------------ */
- /* Typedefs */
/* ------------------------------------------------------------------------ */
-
- /// associated iterator with any nodal field (non filetered)
- class iterator : public iohelper::iterator<T, iterator, Vector<T>> {
+ template <typename T, class Container, class Filter>
+ class NodalField<T, false, Container, Filter> : public dumpers::Field {
public:
- iterator(T * vect, UInt offset, UInt n, UInt stride,
- __attribute__((unused)) const UInt * filter = nullptr)
- :
-
- internal_it(vect), offset(offset), n(n), stride(stride) {}
+ /* ---------------------------------------------------------------------- */
+ /* Typedefs */
+ /* ---------------------------------------------------------------------- */
+ using support_type = UInt;
+
+ /// associated iterator with any nodal field (non filetered)
+ class iterator : public iohelper::iterator<T, iterator, Vector<T>> {
+ public:
+ iterator(T * vect, UInt offset, UInt n, UInt stride,
+ __attribute__((unused)) const UInt * filter = nullptr)
+ : internal_it(vect), offset(offset), n(n), stride(stride) {}
+
+ bool operator!=(const iterator & it) const override {
+ return internal_it != it.internal_it;
+ }
+ iterator & operator++() override {
+ internal_it += offset;
+ return *this;
+ };
+ Vector<T> operator*() override {
+ return Vector<T>(internal_it + stride, n);
+ };
+
+ private:
+ T * internal_it;
+ UInt offset, n, stride;
+ };
- bool operator!=(const iterator & it) const override {
- return internal_it != it.internal_it;
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
+ public:
+ NodalField(const Container & field, UInt n = 0, UInt stride = 0,
+ [[gnu::unused]] const Filter * filter = nullptr)
+ : field(field), n(n), stride(stride), padding(0) {
+ AKANTU_DEBUG_ASSERT(filter == nullptr,
+ "Filter passed to unfiltered NodalField!");
+ if (n == 0) {
+ this->n = field.getNbComponent() - stride;
+ }
}
- iterator & operator++() override {
- internal_it += offset;
- return *this;
- };
- Vector<T> operator*() override {
- return Vector<T>(internal_it + stride, n);
- };
- private:
- T * internal_it;
- UInt offset, n, stride;
- };
+ /* ---------------------------------------------------------------------- */
+ /* Methods */
+ /* ---------------------------------------------------------------------- */
+ public:
+ void registerToDumper(const std::string & id,
+ iohelper::Dumper & dumper) override {
+ dumper.addNodeDataField(id, *this);
+ }
- /* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------ */
- NodalField(const Container & field, UInt n = 0, UInt stride = 0,
- [[gnu::unused]] const Filter * filter = nullptr)
- : field(field), n(n), stride(stride), padding(0) {
- AKANTU_DEBUG_ASSERT(filter == nullptr,
- "Filter passed to unfiltered NodalField!");
- if (n == 0) {
- this->n = field.getNbComponent() - stride;
+ inline iterator begin() {
+ return iterator(field.storage(), field.getNbComponent(), n, stride);
}
- }
- /* ------------------------------------------------------------------------ */
- /* Methods */
- /* ------------------------------------------------------------------------ */
- void registerToDumper(const std::string & id,
- iohelper::Dumper & dumper) override {
- dumper.addNodeDataField(id, *this);
- }
-
- inline iterator begin() {
- return iterator(field.storage(), field.getNbComponent(), n, stride);
- }
-
- inline iterator end() {
- return iterator(field.storage() + field.getNbComponent() * field.size(),
- field.getNbComponent(), n, stride);
- }
-
- bool isHomogeneous() override { return true; }
- void checkHomogeneity() override { this->homogeneous = true; }
-
- virtual UInt getDim() {
- if (this->padding) {
- return this->padding;
+ inline iterator end() {
+ return iterator(field.storage() + field.getNbComponent() * field.size(),
+ field.getNbComponent(), n, stride);
}
+
+ bool isHomogeneous() override { return true; }
+ void checkHomogeneity() override { this->homogeneous = true; }
+
+ virtual UInt getDim() {
+ if (this->padding) {
+ return this->padding;
+ }
return n;
+ }
- }
+ void setPadding(UInt padding) { this->padding = padding; }
- void setPadding(UInt padding) { this->padding = padding; }
+ UInt size() { return field.size(); }
- UInt size() { return field.size(); }
+ iohelper::DataType getDataType() { return iohelper::getDataType<T>(); }
- iohelper::DataType getDataType() { return iohelper::getDataType<T>(); }
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
+ private:
+ const Container & field;
+ UInt n, stride;
+ UInt padding;
+ };
/* ------------------------------------------------------------------------ */
- /* Class Members */
- /* ------------------------------------------------------------------------ */
+ template <typename T, class Container, class Filter>
+ class NodalField<T, true, Container, Filter> : public dumpers::Field {
+ /* ---------------------------------------------------------------------- */
+ /* Typedefs */
+ /* ---------------------------------------------------------------------- */
+ public:
+ class iterator : public iohelper::iterator<T, iterator, Vector<T>> {
-private:
- const Container & field;
- UInt n, stride;
- UInt padding;
-};
+ public:
+ iterator(T * const vect, UInt _offset, UInt _n, UInt _stride,
+ const UInt * filter)
+ :
-/* -------------------------------------------------------------------------- */
-template <typename T, class Container, class Filter>
-class NodalField<T, true, Container, Filter> : public dumpers::Field {
+ internal_it(vect), offset(_offset), n(_n), stride(_stride),
+ filter(filter) {}
- /* ------------------------------------------------------------------------ */
- /* Typedefs */
- /* ------------------------------------------------------------------------ */
+ bool operator!=(const iterator & it) const override {
+ return filter != it.filter;
+ }
-public:
- class iterator : public iohelper::iterator<T, iterator, Vector<T>> {
+ iterator & operator++() override {
+ ++filter;
+ return *this;
+ }
- public:
- iterator(T * const vect, UInt _offset, UInt _n, UInt _stride,
- const UInt * filter)
- :
+ Vector<T> operator*() override {
+ return Vector<T>(internal_it + *(filter)*offset + stride, n);
+ }
- internal_it(vect), offset(_offset), n(_n), stride(_stride),
- filter(filter) {}
+ private:
+ T * const internal_it;
+ UInt offset, n, stride;
+ const UInt * filter;
+ };
- bool operator!=(const iterator & it) const override {
- return filter != it.filter;
+ /* ---------------------------------------------------------------------- */
+ /* Constructors/Destructors */
+ /* ---------------------------------------------------------------------- */
+ public:
+ NodalField(const Container & _field, UInt _n = 0, UInt _stride = 0,
+ const Filter * filter = NULL)
+ : field(_field), n(_n), stride(_stride), filter(filter), padding(0) {
+ AKANTU_DEBUG_ASSERT(this->filter != nullptr,
+ "No filter passed to filtered NodalField!");
+
+ AKANTU_DEBUG_ASSERT(this->filter->getNbComponent() == 1,
+ "Multi-component filter given to NodalField ("
+ << this->filter->getNbComponent()
+ << " components detected, sould be 1");
+ if (n == 0) {
+ this->n = field.getNbComponent() - stride;
+ }
}
- iterator & operator++() override {
- ++filter;
- return *this;
+ /* ---------------------------------------------------------------------- */
+ /* Methods */
+ /* ---------------------------------------------------------------------- */
+ public:
+ void registerToDumper(const std::string & id,
+ iohelper::Dumper & dumper) override {
+ dumper.addNodeDataField(id, *this);
}
- Vector<T> operator*() override {
- return Vector<T>(internal_it + *(filter)*offset + stride, n);
+ inline iterator begin() {
+ return iterator(field.storage(), field.getNbComponent(), n, stride,
+ filter->storage());
}
- private:
- T * const internal_it;
- UInt offset, n, stride;
- const UInt * filter;
- };
-
- /* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------ */
-
- NodalField(const Container & _field, UInt _n = 0, UInt _stride = 0,
- const Filter * filter = NULL)
- : field(_field), n(_n), stride(_stride), filter(filter), padding(0) {
- AKANTU_DEBUG_ASSERT(this->filter != nullptr,
- "No filter passed to filtered NodalField!");
-
- AKANTU_DEBUG_ASSERT(this->filter->getNbComponent() == 1,
- "Multi-component filter given to NodalField ("
- << this->filter->getNbComponent()
- << " components detected, sould be 1");
- if (n == 0) {
- this->n = field.getNbComponent() - stride;
+ inline iterator end() {
+ return iterator(field.storage(), field.getNbComponent(), n, stride,
+ filter->storage() + filter->size());
}
- }
- /* ------------------------------------------------------------------------ */
- /* Methods */
- /* ------------------------------------------------------------------------ */
+ bool isHomogeneous() override { return true; }
+ void checkHomogeneity() override { this->homogeneous = true; }
- void registerToDumper(const std::string & id,
- iohelper::Dumper & dumper) override {
- dumper.addNodeDataField(id, *this);
- }
-
- inline iterator begin() {
- return iterator(field.storage(), field.getNbComponent(), n, stride,
- filter->storage());
- }
-
- inline iterator end() {
- return iterator(field.storage(), field.getNbComponent(), n, stride,
- filter->storage() + filter->size());
- }
-
- bool isHomogeneous() override { return true; }
- void checkHomogeneity() override { this->homogeneous = true; }
-
- virtual UInt getDim() {
- if (this->padding) {
- return this->padding;
- }
+ virtual UInt getDim() {
+ if (this->padding) {
+ return this->padding;
+ }
return n;
+ }
- }
-
- void setPadding(UInt padding) { this->padding = padding; }
-
- UInt size() { return filter->size(); }
+ void setPadding(UInt padding) { this->padding = padding; }
- iohelper::DataType getDataType() { return iohelper::getDataType<T>(); }
+ UInt size() { return filter->size(); }
- /* ------------------------------------------------------------------------ */
- /* Class Members */
- /* ------------------------------------------------------------------------ */
+ iohelper::DataType getDataType() { return iohelper::getDataType<T>(); }
-private:
- const Container & field;
- UInt n, stride;
- const Filter * filter;
+ /* ---------------------------------------------------------------------- */
+ /* Class Members */
+ /* ---------------------------------------------------------------------- */
+ private:
+ const Container & field;
+ UInt n, stride;
+ const Filter * filter;
- UInt padding;
-};
+ UInt padding;
+ };
} // namespace dumpers
} // namespace akantu
/* -------------------------------------------------------------------------- */
#endif /* AKANTU_DUMPER_NODAL_FIELD_HH_ */
diff --git a/src/io/dumper/dumper_type_traits.hh b/src/io/dumper/dumper_type_traits.hh
index cb608b329..810b78b25 100644
--- a/src/io/dumper/dumper_type_traits.hh
+++ b/src/io/dumper/dumper_type_traits.hh
@@ -1,89 +1,88 @@
/**
* @file dumper_type_traits.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date creation: Tue Sep 02 2014
* @date last modification: Wed Nov 08 2017
*
* @brief Type traits for field properties
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
#ifndef AKANTU_DUMPER_TYPE_TRAITS_HH_
#define AKANTU_DUMPER_TYPE_TRAITS_HH_
/* -------------------------------------------------------------------------- */
#include "element_type_map.hh"
#include "element_type_map_filter.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace dumpers {
/* ------------------------------------------------------------------------ */
template <class data, class ret, class field> struct TypeTraits {
-
//! the stored data (real, int, uint, ...)
using data_type = data;
//! the type returned by the operator *
using return_type = ret;
//! the field type (ElementTypeMap or ElementTypeMapFilter)
using field_type = field;
//! the type over which we iterate
using it_type = typename field_type::type;
//! the type of array (Array<T> or ArrayFilter<T>)
using array_type = typename field_type::array_type;
//! the iterator over the array
using array_iterator = typename array_type::const_vector_iterator;
};
/* ------------------------------------------------------------------------ */
// specialization for the case in which input and output types are the same
template <class T, template <class> class ret, bool filtered>
struct SingleType : public TypeTraits<T, ret<T>, ElementTypeMapArray<T>> {};
/* ------------------------------------------------------------------------ */
// same as before but for filtered data
template <class T, template <class> class ret>
struct SingleType<T, ret, true>
: public TypeTraits<T, ret<T>, ElementTypeMapArrayFilter<T>> {};
/* ------------------------------------------------------------------------ */
// specialization for the case in which input and output types are different
template <class it_type, class data_type, template <class> class ret,
bool filtered>
struct DualType : public TypeTraits<data_type, ret<data_type>,
ElementTypeMapArray<it_type>> {};
/* ------------------------------------------------------------------------ */
// same as before but for filtered data
template <class it_type, class data_type, template <class> class ret>
struct DualType<it_type, data_type, ret, true>
: public TypeTraits<data_type, ret<data_type>,
ElementTypeMapArrayFilter<it_type>> {};
/* ------------------------------------------------------------------------ */
} // namespace dumpers
} // namespace akantu
/* -------------------------------------------------------------------------- */
#endif /* AKANTU_DUMPER_TYPE_TRAITS_HH_ */
diff --git a/src/model/contact_mechanics/contact_detector_inline_impl.cc b/src/model/contact_mechanics/contact_detector_inline_impl.cc
index 96fadcb77..8004ed90d 100644
--- a/src/model/contact_mechanics/contact_detector_inline_impl.cc
+++ b/src/model/contact_mechanics/contact_detector_inline_impl.cc
@@ -1,333 +1,332 @@
/**
* @file contact_detection.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Mon Apr 29 2019
* @date last modification: Mon Apr 29 2019
*
* @brief inine implementation of the contact detector class
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "contact_detector.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__
#define __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__
namespace akantu {
/* -------------------------------------------------------------------------- */
inline bool
ContactDetector::checkValidityOfProjection(Vector<Real> & projection) {
-
- UInt nb_xi_inside = 0;
Real tolerance = 1e-3;
for (auto xi : projection) {
if ((xi < -1.0 - tolerance) or (xi > 1.0 + tolerance)) {
return false;
+ }
}
return true;
}
/* -------------------------------------------------------------------------- */
inline void ContactDetector::coordinatesOfElement(const Element & el,
Matrix<Real> & coords) {
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type);
Vector<UInt> connect = mesh.getConnectivity(el.type, _not_ghost)
.begin(nb_nodes_per_element)[el.element];
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connect[n];
for (UInt s : arange(spatial_dimension)) {
coords(s, n) = this->positions(node, s);
}
}
}
/* -------------------------------------------------------------------------- */
inline void ContactDetector::computeCellSpacing(Vector<Real> & spacing) {
for (UInt s : arange(spatial_dimension))
spacing(s) = std::sqrt(2.0) * max_dd;
}
/* -------------------------------------------------------------------------- */
inline void
ContactDetector::constructBoundingBox(BBox & bbox,
const Array<UInt> & nodes_list) {
auto to_bbox = [&](UInt node) {
Vector<Real> pos(spatial_dimension);
for (UInt s : arange(spatial_dimension)) {
pos(s) = this->positions(node, s);
}
bbox += pos;
};
std::for_each(nodes_list.begin(), nodes_list.end(), to_bbox);
auto & lower_bound = bbox.getLowerBounds();
auto & upper_bound = bbox.getUpperBounds();
for (UInt s : arange(spatial_dimension)) {
lower_bound(s) -= this->max_bb;
upper_bound(s) += this->max_bb;
}
AKANTU_DEBUG_INFO("BBox" << bbox);
}
/* -------------------------------------------------------------------------- */
inline void ContactDetector::constructGrid(SpatialGrid<UInt> & grid,
BBox & bbox,
const Array<UInt> & nodes_list) {
auto to_grid = [&](UInt node) {
Vector<Real> pos(spatial_dimension);
for (UInt s : arange(spatial_dimension)) {
pos(s) = this->positions(node, s);
}
if (bbox.contains(pos)) {
grid.insert(node, pos);
}
};
std::for_each(nodes_list.begin(), nodes_list.end(), to_grid);
}
/* -------------------------------------------------------------------------- */
inline void ContactDetector::computeMaximalDetectionDistance() {
AKANTU_DEBUG_IN();
Real elem_size;
Real max_elem_size = std::numeric_limits<Real>::min();
Real min_elem_size = std::numeric_limits<Real>::max();
auto & master_nodes = this->surface_selector->getMasterList();
for (auto & master : master_nodes) {
Array<Element> elements;
this->mesh.getAssociatedElements(master, elements);
for (auto element : elements) {
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type);
Matrix<Real> elem_coords(spatial_dimension, nb_nodes_per_element);
this->coordinatesOfElement(element, elem_coords);
elem_size = FEEngine::getElementInradius(elem_coords, element.type);
max_elem_size = std::max(max_elem_size, elem_size);
min_elem_size = std::min(min_elem_size, elem_size);
}
}
AKANTU_DEBUG_INFO("The maximum element size : " << max_elem_size);
this->min_dd = min_elem_size;
this->max_dd = max_elem_size;
this->max_bb = max_elem_size;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline Vector<UInt>
ContactDetector::constructConnectivity(UInt & slave, const Element & master) {
Vector<UInt> master_conn =
const_cast<const Mesh &>(this->mesh).getConnectivity(master);
Vector<UInt> elem_conn(master_conn.size() + 1);
elem_conn[0] = slave;
for (UInt i = 1; i < elem_conn.size(); ++i) {
elem_conn[i] = master_conn[i - 1];
}
return elem_conn;
}
/* -------------------------------------------------------------------------- */
inline void ContactDetector::computeNormalOnElement(const Element & element,
Vector<Real> & normal) {
Matrix<Real> vectors(spatial_dimension, spatial_dimension - 1);
this->vectorsAlongElement(element, vectors);
switch (this->spatial_dimension) {
case 2: {
Math::normal2(vectors.storage(), normal.storage());
break;
}
case 3: {
Math::normal3(vectors(0).storage(), vectors(1).storage(), normal.storage());
break;
}
default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); }
}
// to ensure that normal is always outwards from master surface
const auto & element_to_subelement =
mesh.getElementToSubelement(element.type)(element.element);
Vector<Real> outside(spatial_dimension);
mesh.getBarycenter(element, outside);
// check if mesh facets exists for cohesive elements contact
Vector<Real> inside(spatial_dimension);
if (mesh.isMeshFacets()) {
mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside);
} else {
mesh.getBarycenter(element_to_subelement[0], inside);
}
Vector<Real> inside_to_outside = outside - inside;
auto projection = inside_to_outside.dot(normal);
if (projection < 0) {
normal *= -1.0;
}
}
/* -------------------------------------------------------------------------- */
inline void ContactDetector::vectorsAlongElement(const Element & el,
Matrix<Real> & vectors) {
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type);
Matrix<Real> coords(spatial_dimension, nb_nodes_per_element);
this->coordinatesOfElement(el, coords);
for(auto i : arange(spatial_dimension - 1)) {
vectors(i) = Vector<Real>(coords(i + 1)) - Vector<Real>(coords(0));
}
}
/* -------------------------------------------------------------------------- */
inline Real ContactDetector::computeGap(Vector<Real> & slave,
Vector<Real> & master) {
Vector<Real> slave_to_master(spatial_dimension);
slave_to_master = master - slave;
Real gap = slave_to_master.norm();
return gap;
}
/* -------------------------------------------------------------------------- */
inline void
ContactDetector::filterBoundaryElements(Array<Element> & elements,
Array<Element> & boundary_elements) {
for (auto elem : elements) {
const auto & element_to_subelement =
mesh.getElementToSubelement(elem.type)(elem.element);
// for regular boundary elements
if (element_to_subelement.size() == 1 and
element_to_subelement[0].kind() == _ek_regular) {
boundary_elements.push_back(elem);
continue;
}
// for cohesive boundary elements
UInt nb_subelements_regular = 0;
for (auto subelem : element_to_subelement) {
if (subelem == ElementNull)
continue;
if (subelem.kind() == _ek_regular)
++nb_subelements_regular;
}
auto nb_subelements = element_to_subelement.size();
if (nb_subelements_regular < nb_subelements)
boundary_elements.push_back(elem);
}
}
/* -------------------------------------------------------------------------- */
inline bool
ContactDetector::isValidSelfContact(const UInt & slave_node, const Real & gap,
const Vector<Real> & normal) {
UInt master_node;
// finding the master node corresponding to slave node
for (auto & pair : contact_pairs) {
if (pair.first == slave_node) {
master_node = pair.second;
break;
}
}
Array<Element> slave_elements;
this->mesh.getAssociatedElements(slave_node, slave_elements);
// Check 1 : master node is not connected to elements connected to
// slave node
Vector<Real> slave_normal(spatial_dimension);
for (auto & element : slave_elements) {
if (element.kind() != _ek_regular)
continue;
Vector<UInt> connectivity =
const_cast<const Mesh &>(this->mesh).getConnectivity(element);
// finding the normal at slave node by averaging of normals
Vector<Real> normal(spatial_dimension);
GeometryUtils::normal(mesh, positions, element, normal);
slave_normal = slave_normal + normal;
auto node_iter =
std::find(connectivity.begin(), connectivity.end(), master_node);
if (node_iter != connectivity.end())
return false;
}
// Check 2 : if gap is twice the size of smallest element
if (std::abs(gap) > 2.0 * min_dd)
return false;
// Check 3 : check the directions of normal at slave node and at
// master element, should be in opposite directions
auto norm = slave_normal.norm();
if (norm != 0)
slave_normal /= norm;
auto product = slave_normal.dot(normal);
if (product >= 0)
return false;
return true;
}
} // namespace akantu
#endif /* __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__ */
diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc
index 333b51d99..54b27e4b1 100644
--- a/src/model/contact_mechanics/contact_mechanics_model.cc
+++ b/src/model/contact_mechanics/contact_mechanics_model.cc
@@ -1,754 +1,760 @@
/**
* @file coontact_mechanics_model.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Tue May 08 2012
* @date last modification: Wed Feb 21 2018
*
* @brief Contact mechanics model
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "contact_mechanics_model.hh"
#include "boundary_condition_functor.hh"
#include "dumpable_inline_impl.hh"
+#include "group_manager_inline_impl.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
-#include "group_manager_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
/* -------------------------------------------------------------------------- */
#include <algorithm>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
ContactMechanicsModel::ContactMechanicsModel(
Mesh & mesh, UInt dim, const ID & id,
std::shared_ptr<DOFManager> dof_manager, const ModelType model_type)
: Model(mesh, model_type, dof_manager, dim, id) {
AKANTU_DEBUG_IN();
this->registerFEEngineObject<MyFEEngineType>("ContactMechanicsModel", mesh,
Model::spatial_dimension);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("contact_mechanics", id, true);
this->mesh.addDumpMeshToDumper("contact_mechanics", mesh,
Model::spatial_dimension, _not_ghost,
_ek_regular);
#endif
this->registerDataAccessor(*this);
this->detector =
std::make_unique<ContactDetector>(this->mesh, id + ":contact_detector");
registerFEEngineObject<MyFEEngineType>("ContactFacetsFEEngine", mesh,
- Model::spatial_dimension - 1);
-
+ Model::spatial_dimension - 1);
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
ContactMechanicsModel::~ContactMechanicsModel() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::initFullImpl(const ModelOptions & options) {
Model::initFullImpl(options);
// initalize the resolutions
if (this->parser.getLastParsedFile() != "") {
this->instantiateResolutions();
this->initResolutions();
}
this->initBC(*this, *displacement, *displacement_increment, *external_force);
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::instantiateResolutions() {
ParserSection model_section;
bool is_empty;
std::tie(model_section, is_empty) = this->getParserSection();
if (not is_empty) {
auto model_resolutions =
model_section.getSubSections(ParserType::_contact_resolution);
for (const auto & section : model_resolutions) {
this->registerNewResolution(section);
}
}
auto sub_sections =
this->parser.getSubSections(ParserType::_contact_resolution);
for (const auto & section : sub_sections) {
this->registerNewResolution(section);
}
if (resolutions.empty())
AKANTU_EXCEPTION("No contact resolutions where instantiated for the model"
<< getID());
are_resolutions_instantiated = true;
}
/* -------------------------------------------------------------------------- */
Resolution &
ContactMechanicsModel::registerNewResolution(const ParserSection & section) {
std::string res_name;
std::string res_type = section.getName();
std::string opt_param = section.getOption();
try {
std::string tmp = section.getParameter("name");
res_name = tmp; /** this can seem weird, but there is an ambiguous operator
* overload that i couldn't solve. @todo remove the
* weirdness of this code
*/
} catch (debug::Exception &) {
AKANTU_ERROR("A contact resolution of type \'"
<< res_type
<< "\' in the input file has been defined without a name!");
}
Resolution & res = this->registerNewResolution(res_name, res_type, opt_param);
res.parseSection(section);
return res;
}
/* -------------------------------------------------------------------------- */
Resolution & ContactMechanicsModel::registerNewResolution(
const ID & res_name, const ID & res_type, const ID & opt_param) {
AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) ==
resolutions_names_to_id.end(),
"A resolution with this name '"
<< res_name << "' has already been registered. "
<< "Please use unique names for resolutions");
UInt res_count = resolutions.size();
resolutions_names_to_id[res_name] = res_count;
std::stringstream sstr_res;
sstr_res << this->id << ":" << res_count << ":" << res_type;
ID res_id = sstr_res.str();
std::unique_ptr<Resolution> resolution =
ResolutionFactory::getInstance().allocate(res_type, spatial_dimension,
opt_param, *this, res_id);
resolutions.push_back(std::move(resolution));
return *(resolutions.back());
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::initResolutions() {
AKANTU_DEBUG_ASSERT(resolutions.size() != 0,
"No resolutions to initialize !");
if (!are_resolutions_instantiated)
instantiateResolutions();
// \TODO check if each resolution needs a initResolution() method
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::initModel() {
AKANTU_DEBUG_IN();
-
+
getFEEngine("ContactMechanicsModel").initShapeFunctions(_not_ghost);
getFEEngine("ContactMechanicsModel").initShapeFunctions(_ghost);
getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_not_ghost);
getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_ghost);
- AKANTU_DEBUG_OUT();
+ AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) {
return dynamic_cast<FEEngine &>(
getFEEngineClassBoundary<MyFEEngineType>(name));
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::initSolver(
TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType) {
// for alloc type of solvers
this->allocNodalField(this->displacement, spatial_dimension, "displacement");
this->allocNodalField(this->displacement_increment, spatial_dimension,
"displacement_increment");
this->allocNodalField(this->internal_force, spatial_dimension,
"internal_force");
this->allocNodalField(this->external_force, spatial_dimension,
"external_force");
- this->allocNodalField(this->normal_force, spatial_dimension,
- "normal_force");
+ this->allocNodalField(this->normal_force, spatial_dimension, "normal_force");
this->allocNodalField(this->tangential_force, spatial_dimension,
- "tangential_force");
-
+ "tangential_force");
+
this->allocNodalField(this->gaps, 1, "gaps");
this->allocNodalField(this->nodal_area, 1, "areas");
this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs");
this->allocNodalField(this->contact_state, 1, "contact_state");
- this->allocNodalField(this->previous_master_elements, 1, "previous_master_elements");
-
+ this->allocNodalField(this->previous_master_elements, 1,
+ "previous_master_elements");
+
this->allocNodalField(this->normals, spatial_dimension, "normals");
auto surface_dimension = spatial_dimension - 1;
- this->allocNodalField(this->tangents, surface_dimension*spatial_dimension,
- "tangents");
- this->allocNodalField(this->projections, surface_dimension,
- "projections");
+ this->allocNodalField(this->tangents, surface_dimension * spatial_dimension,
+ "tangents");
+ this->allocNodalField(this->projections, surface_dimension, "projections");
this->allocNodalField(this->previous_projections, surface_dimension,
- "previous_projections");
- this->allocNodalField(this->previous_tangents, surface_dimension*spatial_dimension,
- "previous_tangents");
+ "previous_projections");
+ this->allocNodalField(this->previous_tangents,
+ surface_dimension * spatial_dimension,
+ "previous_tangents");
this->allocNodalField(this->tangential_tractions, surface_dimension,
- "tangential_tractions");
+ "tangential_tractions");
this->allocNodalField(this->previous_tangential_tractions, surface_dimension,
- "previous_tangential_tractions");
-
+ "previous_tangential_tractions");
+
// todo register multipliers as dofs for lagrange multipliers
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_lumped_mass: {
return std::make_tuple("explicit_lumped",
TimeStepSolverType::_dynamic_lumped);
}
case _explicit_consistent_mass: {
return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
}
case _static: {
return std::make_tuple("static", TimeStepSolverType::_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
}
default:
return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
}
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_dynamic_lumped: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type =
NonLinearSolverType::_newton_raphson_contact;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
break;
}
return options;
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleResidual() {
AKANTU_DEBUG_IN();
/* ------------------------------------------------------------------------ */
// computes the internal forces
this->assembleInternalForces();
/* ------------------------------------------------------------------------ */
this->getDOFManager().assembleToResidual("displacement",
*this->internal_force, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the contact forces");
UInt nb_nodes = mesh.getNbNodes();
this->internal_force->clear();
this->normal_force->clear();
this->tangential_force->clear();
internal_force->resize(nb_nodes, 0.);
normal_force->resize(nb_nodes, 0.);
tangential_force->resize(nb_nodes, 0.);
// assemble the forces due to contact
auto assemble = [&](auto && ghost_type) {
for (auto & resolution : resolutions) {
resolution->assembleInternalForces(ghost_type);
}
};
AKANTU_DEBUG_INFO("Assemble residual for local elements");
assemble(_not_ghost);
// assemble the stresses due to ghost elements
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
- //assemble(_ghost);
+ // assemble(_ghost);
// TODO : uncomment when developing code for parallelization,
// currently it addes the force twice for not ghost elements
// hence source of error
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::search() {
- // save the previous state
+ // save the previous state
this->savePreviousState();
-
+
contact_elements.clear();
- // this needs to be resized if cohesive elements are added
+ // this needs to be resized if cohesive elements are added
UInt nb_nodes = mesh.getNbNodes();
auto resize_arrays = [&](auto & internal_array) {
internal_array->resize(nb_nodes);
internal_array->zero();
};
resize_arrays(gaps);
resize_arrays(normals);
resize_arrays(tangents);
resize_arrays(projections);
resize_arrays(tangential_tractions);
resize_arrays(contact_state);
resize_arrays(nodal_area);
resize_arrays(external_force);
-
- this->detector->search(contact_elements, *gaps,
- *normals, *tangents, *projections);
+
+ this->detector->search(contact_elements, *gaps, *normals, *tangents,
+ *projections);
// intepenetration value must be positive for contact mechanics
// model to work by default the gap value from detector is negative
- std::for_each((*gaps).begin(), (*gaps).end(),
- [](Real & gap){ gap *= -1.; });
-
+ std::for_each((*gaps).begin(), (*gaps).end(), [](Real & gap) { gap *= -1.; });
+
if (contact_elements.size() != 0) {
- this->computeNodalAreas();
+ this->computeNodalAreas();
}
-
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::savePreviousState() {
AKANTU_DEBUG_IN();
- // saving previous natural projections
+ // saving previous natural projections
(*previous_projections).copy(*projections);
// saving previous tangents
(*previous_tangents).copy(*tangents);
// saving previous tangential traction
(*previous_tangential_tractions).copy(*tangential_tractions);
previous_master_elements->clear();
previous_master_elements->resize(projections->size());
previous_master_elements->set(ElementNull);
for (auto & element : contact_elements) {
(*previous_master_elements)[element.slave] = element.master;
}
AKANTU_DEBUG_OUT();
}
-
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::computeNodalAreas(GhostType ghost_type) {
UInt nb_nodes = mesh.getNbNodes();
nodal_area->resize(nb_nodes);
nodal_area->zero();
external_force->resize(nb_nodes);
external_force->zero();
-
- auto & fem_boundary = getFEEngineClassBoundary<MyFEEngineType>("ContactMechanicsModel");
+
+ auto & fem_boundary =
+ getFEEngineClassBoundary<MyFEEngineType>("ContactMechanicsModel");
fem_boundary.initShapeFunctions(getPositions(), _not_ghost);
fem_boundary.initShapeFunctions(getPositions(), _ghost);
-
+
fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_ghost);
-
-
+
IntegrationPoint quad_point;
quad_point.ghost_type = ghost_type;
-
+
auto & group = mesh.getElementGroup("contact_surface");
UInt nb_degree_of_freedom = external_force->getNbComponent();
for (auto && type : group.elementTypes(spatial_dimension - 1, ghost_type)) {
const auto & element_ids = group.getElements(type, ghost_type);
-
- UInt nb_quad_points =
- fem_boundary.getNbIntegrationPoints(type, ghost_type);
+
+ UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type, ghost_type);
UInt nb_elements = element_ids.size();
-
- UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
+
+ UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
Array<Real> dual_before_integ(nb_elements * nb_quad_points,
- nb_degree_of_freedom, 0.);
+ nb_degree_of_freedom, 0.);
Array<Real> quad_coords(nb_elements * nb_quad_points, spatial_dimension);
const auto & normals_on_quad =
- fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type);
+ fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type);
auto normals_begin = normals_on_quad.begin(spatial_dimension);
decltype(normals_begin) normals_iter;
auto quad_coords_iter = quad_coords.begin(spatial_dimension);
auto dual_iter = dual_before_integ.begin(nb_degree_of_freedom);
quad_point.type = type;
Element subelement;
subelement.type = type;
subelement.ghost_type = ghost_type;
for (auto el : element_ids) {
subelement.element = el;
const auto & element_to_subelement =
- mesh.getElementToSubelement(type)(el);
+ mesh.getElementToSubelement(type)(el);
Vector<Real> outside(spatial_dimension);
mesh.getBarycenter(subelement, outside);
Vector<Real> inside(spatial_dimension);
if (mesh.isMeshFacets()) {
- mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside);
- }
- else {
- mesh.getBarycenter(element_to_subelement[0], inside);
+ mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside);
+ } else {
+ mesh.getBarycenter(element_to_subelement[0], inside);
}
Vector<Real> inside_to_outside(spatial_dimension);
inside_to_outside = outside - inside;
-
+
normals_iter = normals_begin + el * nb_quad_points;
-
+
quad_point.element = el;
for (auto q : arange(nb_quad_points)) {
- quad_point.num_point = q;
- auto ddot = inside_to_outside.dot(*normals_iter);
- Vector<Real> normal(*normals_iter);
- if (ddot < 0)
- normal *= -1.0;
-
- (*dual_iter).mul<false>(Matrix<Real>::eye(spatial_dimension, 1),
- normal);
- ++dual_iter;
- ++quad_coords_iter;
- ++normals_iter;
+ quad_point.num_point = q;
+ auto ddot = inside_to_outside.dot(*normals_iter);
+ Vector<Real> normal(*normals_iter);
+ if (ddot < 0)
+ normal *= -1.0;
+
+ (*dual_iter)
+ .mul<false>(Matrix<Real>::eye(spatial_dimension, 1), normal);
+ ++dual_iter;
+ ++quad_coords_iter;
+ ++normals_iter;
}
}
Array<Real> dual_by_shapes(nb_elements * nb_quad_points,
- nb_degree_of_freedom * nb_nodes_per_element);
+ nb_degree_of_freedom * nb_nodes_per_element);
- fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type,
- ghost_type, element_ids);
+ fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type, ghost_type,
+ element_ids);
Array<Real> dual_by_shapes_integ(nb_elements, nb_degree_of_freedom *
- nb_nodes_per_element);
+ nb_nodes_per_element);
fem_boundary.integrate(dual_by_shapes, dual_by_shapes_integ,
- nb_degree_of_freedom * nb_nodes_per_element, type,
- ghost_type, element_ids);
+ nb_degree_of_freedom * nb_nodes_per_element, type,
+ ghost_type, element_ids);
this->getDOFManager().assembleElementalArrayLocalArray(
- dual_by_shapes_integ, *external_force, type, ghost_type, 1., element_ids);
+ dual_by_shapes_integ, *external_force, type, ghost_type, 1.,
+ element_ids);
}
for (auto && tuple :
- zip(*nodal_area,
- make_view(*external_force, spatial_dimension))) {
-
+ zip(*nodal_area, make_view(*external_force, spatial_dimension))) {
+
auto & area = std::get<0>(tuple);
Vector<Real> force(std::get<1>(tuple));
area = force.norm();
}
-
+
this->external_force->clear();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space(indent, AKANTU_INDENT);
stream << space << "Contact Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << Model::spatial_dimension
<< std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + resolutions [" << std::endl;
for (auto & resolution : resolutions) {
resolution->printself(stream, indent + 1);
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) {
if (matrix_id == "K")
return _symmetric;
return _mt_not_defined;
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
}
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix");
if (!this->getDOFManager().hasMatrix("K")) {
this->getDOFManager().getNewMatrix("K", getMatrixType("K"));
}
for (auto & resolution : resolutions) {
resolution->assembleStiffnessMatrix(_not_ghost);
}
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {
AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::beforeSolveStep() {
- for (auto & resolution : resolutions)
+ for (auto & resolution : resolutions) {
resolution->beforeSolveStep();
+ }
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::afterSolveStep(bool converged) {
- for (auto & resolution : resolutions)
+ for (auto & resolution : resolutions) {
resolution->afterSolveStep(converged);
+ }
}
-
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
ContactMechanicsModel::createNodalFieldBool(const std::string &,
- const std::string &, bool) {
-
+ const std::string &, bool) {
return nullptr;
}
-
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
ContactMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
-
std::map<std::string, Array<Real> *> real_nodal_fields;
- real_nodal_fields["contact_force"] = this->internal_force.get();
- real_nodal_fields["normal_force"] = this->normal_force.get();
+ real_nodal_fields["contact_force"] = this->internal_force.get();
+ real_nodal_fields["normal_force"] = this->normal_force.get();
real_nodal_fields["tangential_force"] = this->tangential_force.get();
- real_nodal_fields["blocked_dofs"] = this->blocked_dofs.get();
- real_nodal_fields["normals"] = this->normals.get();
- real_nodal_fields["tangents"] = this->tangents.get();
- real_nodal_fields["gaps"] = this->gaps.get();
- real_nodal_fields["areas"] = this->nodal_area.get();
- real_nodal_fields["contact_state"] = this->contact_state.get();
+ real_nodal_fields["blocked_dofs"] = this->blocked_dofs.get();
+ real_nodal_fields["normals"] = this->normals.get();
+ real_nodal_fields["tangents"] = this->tangents.get();
+ real_nodal_fields["gaps"] = this->gaps.get();
+ real_nodal_fields["areas"] = this->nodal_area.get();
real_nodal_fields["tangential_traction"] = this->tangential_tractions.get();
-
+
std::shared_ptr<dumpers::Field> field;
- if (padding_flag)
- field = this->mesh.createNodalField(real_nodal_fields[field_name],
- group_name, 3);
- else
+ if (padding_flag) {
field = this->mesh.createNodalField(real_nodal_fields[field_name],
- group_name);
+ group_name, 3);
+ } else {
+ field =
+ this->mesh.createNodalField(real_nodal_fields[field_name], group_name);
+ }
+ return field;
+}
+
+/* -------------------------------------------------------------------------- */
+std::shared_ptr<dumpers::Field>
+ContactMechanicsModel::createNodalFieldUInt(const std::string & field_name,
+ const std::string & group_name,
+ bool /*padding_flag*/) {
+ std::shared_ptr<dumpers::Field> field;
+ if (field_name == "contact_state") {
+ auto && func =
+ std::make_unique<dumpers::ComputeUIntFromEnum<ContactState>>();
+ field = mesh.createNodalField(this->contact_state.get(), group_name);
+ field =
+ dumpers::FieldComputeProxy::createFieldCompute(field, std::move(func));
+ }
return field;
}
#else
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
ContactMechanicsModel::createNodalFieldBool(const std::string &,
- const std::string &, bool) {
-
+ const std::string &, bool) {
return nullptr;
}
-
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
-ContactMechanicsModel::createNodalFieldReal(const std::string & ,
- const std::string & , bool) {
+ContactMechanicsModel::createNodalFieldReal(const std::string &,
+ const std::string &, bool) {
return nullptr;
}
-
#endif
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(const std::string & dumper_name) {
mesh.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(const std::string & dumper_name, UInt step) {
mesh.dump(dumper_name, step);
}
/* ------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(const std::string & dumper_name, Real time,
UInt step) {
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump() { mesh.dump(); }
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(UInt step) { mesh.dump(step); }
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::dump(Real time, UInt step) {
mesh.dump(time, step);
}
/* -------------------------------------------------------------------------- */
UInt ContactMechanicsModel::getNbData(
const Array<Element> & elements, const SynchronizationTag & /*tag*/) const {
AKANTU_DEBUG_IN();
UInt size = 0;
UInt nb_nodes_per_element = 0;
for (const Element & el : elements) {
nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/,
const Array<Element> & /*elements*/,
const SynchronizationTag & /*tag*/) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/,
const Array<Element> & /*elements*/,
const SynchronizationTag & /*tag*/) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt ContactMechanicsModel::getNbData(
const Array<UInt> & dofs, const SynchronizationTag & /*tag*/) const {
AKANTU_DEBUG_IN();
UInt size = 0;
AKANTU_DEBUG_OUT();
return size * dofs.size();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/,
const Array<UInt> & /*dofs*/,
const SynchronizationTag & /*tag*/) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/,
const Array<UInt> & /*dofs*/,
const SynchronizationTag & /*tag*/) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
} // namespace akantu
diff --git a/src/model/contact_mechanics/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh
index b51a8cdb4..d0efb2b5d 100644
--- a/src/model/contact_mechanics/contact_mechanics_model.hh
+++ b/src/model/contact_mechanics/contact_mechanics_model.hh
@@ -1,379 +1,384 @@
/**
* @file contact_mechanics_model.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Model of Contact Mechanics
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition.hh"
#include "contact_detector.hh"
#include "data_accessor.hh"
#include "fe_engine.hh"
#include "model.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__
#define __AKANTU_CONTACT_MECHANICS_MODEL_HH__
namespace akantu {
class Resolution;
template <ElementKind kind, class IntegrationOrderFunctor>
class IntegratorGauss;
template <ElementKind kind> class ShapeLagrange;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
class ContactMechanicsModel : public Model,
public DataAccessor<Element>,
public DataAccessor<UInt>,
public BoundaryCondition<ContactMechanicsModel> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
public:
ContactMechanicsModel(
Mesh & mesh, UInt spatial_dimension = _all_dimensions,
- const ID & id = "contact_mechanics_model", std::shared_ptr<DOFManager> dof_manager = nullptr,
+ const ID & id = "contact_mechanics_model",
+ std::shared_ptr<DOFManager> dof_manager = nullptr,
const ModelType model_type = ModelType::_contact_mechanics_model);
~ContactMechanicsModel() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// initialize completely the model
void initFullImpl(const ModelOptions & options) override;
/// allocate all vectors
void initSolver(TimeStepSolverType, NonLinearSolverType) override;
/// initialize all internal arrays for resolutions
void initResolutions();
/// initialize the modelType
void initModel() override;
/// call back for the solver, computes the force residual
void assembleResidual() override;
/// get the type of matrix needed
MatrixType getMatrixType(const ID & matrix_id) override;
/// callback for the solver, this assembles different matrices
void assembleMatrix(const ID & matrix_id) override;
/// callback for the solver, this assembles the stiffness matrix
void assembleLumpedMatrix(const ID & matrix_id) override;
/// get some default values for derived classes
std::tuple<ID, TimeStepSolverType>
getDefaultSolverID(const AnalysisMethod & method) override;
ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const;
/// callback for the solver, this is called at beginning of solve
void beforeSolveStep() override;
/// callback for the solver, this is called at end of solve
void afterSolveStep(bool converted = true) override;
/// function to print the containt of the class
void printself(std::ostream & stream, int indent = 0) const override;
/* ------------------------------------------------------------------------ */
/* Contact Detection */
/* ------------------------------------------------------------------------ */
public:
void search();
void computeNodalAreas(GhostType ghost_type = _not_ghost);
/* ------------------------------------------------------------------------ */
/* Contact Resolution */
/* ------------------------------------------------------------------------ */
public:
/// register an empty contact resolution of a given type
Resolution & registerNewResolution(const ID & res_name, const ID & res_type,
const ID & opt_param);
protected:
/// register a resolution in the dynamic database
Resolution & registerNewResolution(const ParserSection & res_section);
/// read the resolution files to instantiate all the resolutions
void instantiateResolutions();
/// save the parameters from previous state
void savePreviousState();
/* ------------------------------------------------------------------------ */
/* Solver Interface */
/* ------------------------------------------------------------------------ */
public:
/// assembles the contact stiffness matrix
virtual void assembleStiffnessMatrix();
/// assembles the contant internal forces
virtual void assembleInternalForces();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
FEEngine & getFEEngineBoundary(const ID & name = "") override;
/* ------------------------------------------------------------------------ */
/* Dumpable interface */
/* ------------------------------------------------------------------------ */
public:
std::shared_ptr<dumpers::Field>
createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
+ std::shared_ptr<dumpers::Field>
+ createNodalFieldUInt(const std::string & field_name,
+ const std::string & group_name,
+ bool padding_flag) override;
+
std::shared_ptr<dumpers::Field>
createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
void dump() override;
virtual void dump(UInt step);
virtual void dump(Real time, UInt step);
virtual void dump(const std::string & dumper_name);
virtual void dump(const std::string & dumper_name, UInt step);
virtual void dump(const std::string & dumper_name, Real time, UInt step);
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) override;
UInt getNbData(const Array<UInt> & dofs,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
const SynchronizationTag & tag) override;
protected:
/// contact detection class
friend class ContactDetector;
/// contact resolution class
friend class Resolution;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt);
/// get the ContactMechanicsModel::displacement vector
AKANTU_GET_MACRO(Displacement, *displacement, Array<Real> &);
/// get the ContactMechanicsModel::increment vector \warn only consistent
/// if ContactMechanicsModel::setIncrementFlagOn has been called before
AKANTU_GET_MACRO(Increment, *displacement_increment, Array<Real> &);
/// get the ContactMechanics::internal_force vector (internal forces)
AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &);
/// get the ContactMechanicsModel::external_force vector (external forces)
AKANTU_GET_MACRO(ExternalForce, *external_force, Array<Real> &);
/// get the ContactMechanics::normal_force vector (normal forces)
AKANTU_GET_MACRO(NormalForce, *normal_force, Array<Real> &);
/// get the ContactMechanics::tangential_force vector (friction forces)
AKANTU_GET_MACRO(TangentialForce, *tangential_force, Array<Real> &);
-
+
/// get the ContactMechanics::traction vector (friction traction)
AKANTU_GET_MACRO(TangentialTractions, *tangential_tractions, Array<Real> &);
/// get the ContactMechanics::previous_tangential_tractions vector
AKANTU_GET_MACRO(PreviousTangentialTractions, *previous_tangential_tractions,
- Array<Real> &);
+ Array<Real> &);
/// get the ContactMechanicsModel::force vector (external forces)
Array<Real> & getForce() {
AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, "
"use getExternalForce instead");
return *external_force;
}
/// get the ContactMechanics::blocked_dofs vector
AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<Real> &);
/// get the ContactMechanics::gaps (contact gaps)
AKANTU_GET_MACRO(Gaps, *gaps, Array<Real> &);
/// get the ContactMechanics::normals (normals on slave nodes)
AKANTU_GET_MACRO(Normals, *normals, Array<Real> &);
/// get the ContactMechanics::tangents (tangents on slave nodes)
AKANTU_GET_MACRO(Tangents, *tangents, Array<Real> &);
/// get the ContactMechanics::previous_tangents (tangents on slave nodes)
AKANTU_GET_MACRO(PreviousTangents, *previous_tangents, Array<Real> &);
-
+
/// get the ContactMechanics::areas (nodal areas)
AKANTU_GET_MACRO(NodalArea, *nodal_area, Array<Real> &);
/// get the ContactMechanics::previous_projections (previous_projections)
AKANTU_GET_MACRO(PreviousProjections, *previous_projections, Array<Real> &);
/// get the ContactMechanics::projections (projections)
AKANTU_GET_MACRO(Projections, *projections, Array<Real> &);
-
+
/// get the ContactMechanics::contact_state vector (no_contact/stick/slip
/// state)
- AKANTU_GET_MACRO(ContactState, *contact_state, Array<Real> &);
+ AKANTU_GET_MACRO(ContactState, *contact_state, Array<ContactState> &);
/// get the ContactMechanics::previous_master_elements
- AKANTU_GET_MACRO(PreviousMasterElements, *previous_master_elements, Array<Element> &);
-
+ AKANTU_GET_MACRO(PreviousMasterElements, *previous_master_elements,
+ Array<Element> &);
+
/// get contact detector
AKANTU_GET_MACRO_NOT_CONST(ContactDetector, *detector, ContactDetector &);
/// get the contact elements
inline Array<ContactElement> & getContactElements() {
return contact_elements;
}
/// get the current positions of the nodes
- inline Array<Real> & getPositions() {
- return detector->getPositions();
- }
+ inline Array<Real> & getPositions() { return detector->getPositions(); }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// tells if the resolutions are instantiated
bool are_resolutions_instantiated;
/// displacements array
std::unique_ptr<Array<Real>> displacement;
/// increment of displacement
std::unique_ptr<Array<Real>> displacement_increment;
/// contact forces array
std::unique_ptr<Array<Real>> internal_force;
/// external forces array
std::unique_ptr<Array<Real>> external_force;
/// normal force array
std::unique_ptr<Array<Real>> normal_force;
/// friction force array
std::unique_ptr<Array<Real>> tangential_force;
/// friction traction array
std::unique_ptr<Array<Real>> tangential_tractions;
-
+
/// previous friction traction array
std::unique_ptr<Array<Real>> previous_tangential_tractions;
-
+
/// boundary vector
std::unique_ptr<Array<Real>> blocked_dofs;
/// array to store gap between slave and master
std::unique_ptr<Array<Real>> gaps;
-
+
/// array to store normals from master to slave
std::unique_ptr<Array<Real>> normals;
/// array to store tangents on the master element
std::unique_ptr<Array<Real>> tangents;
/// array to store previous tangents on the master element
std::unique_ptr<Array<Real>> previous_tangents;
-
+
/// array to store nodal areas
std::unique_ptr<Array<Real>> nodal_area;
/// array to store stick/slip state :
- std::unique_ptr<Array<Real>> contact_state;
+ std::unique_ptr<Array<ContactState>> contact_state;
/// array to store previous projections in covariant basis
std::unique_ptr<Array<Real>> previous_projections;
-
+
// array to store projections in covariant basis
std::unique_ptr<Array<Real>> projections;
-
+
/// contact detection
std::unique_ptr<ContactDetector> detector;
/// list of contact resolutions
std::vector<std::unique_ptr<Resolution>> resolutions;
/// mapping between resolution name and resolution internal id
std::map<std::string, UInt> resolutions_names_to_id;
/// array to store contact elements
Array<ContactElement> contact_elements;
/// array to store previous master elements
std::unique_ptr<Array<Element>> previous_master_elements;
};
} // namespace akantu
/* ------------------------------------------------------------------------ */
/* inline functions */
/* ------------------------------------------------------------------------ */
#include "parser.hh"
#include "resolution.hh"
/* ------------------------------------------------------------------------ */
#endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */
diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc
index 2c349c683..0c952233a 100644
--- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc
+++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc
@@ -1,839 +1,839 @@
/**
* @file resolution_penalty.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Mon Jan 7 2019
* @date last modification: Mon Jan 7 2019
*
* @brief Specialization of the resolution class for the penalty method
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "resolution_penalty.hh"
#include "element_class_helper.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model,
const ID & id)
: Resolution(model, id) {
AKANTU_DEBUG_IN();
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::initialize() {
this->registerParam("epsilon_n", epsilon_n, Real(0.),
_pat_parsable | _pat_modifiable,
"Normal penalty parameter");
this->registerParam("epsilon_t", epsilon_t, Real(0.),
_pat_parsable | _pat_modifiable,
"Tangential penalty parameter");
}
/* -------------------------------------------------------------------------- */
Real ResolutionPenalty::computeNormalTraction(Real & gap) {
return epsilon_n * macaulay(gap);
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeNormalForce(const ContactElement & element,
Vector<Real> & force) {
force.zero();
auto & gaps = model.getGaps();
auto & projections = model.getProjections();
auto & normals = model.getNormals();
auto surface_dimension = spatial_dimension - 1;
Real gap(gaps.begin()[element.slave]);
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & nodal_area = const_cast<Array<Real> &>(model.getNodalArea());
// compute normal traction
Real p_n = computeNormalTraction(gap);
p_n *= nodal_area[element.slave];
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> shape_matric(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
ResolutionUtils::computeShapeFunctionMatric(element, projection,
shape_matric);
force.mul<true>(shape_matric, normal, p_n);
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTangentialForce(const ContactElement & element,
Vector<Real> & force) {
if (mu == 0)
return;
force.zero();
UInt surface_dimension = spatial_dimension - 1;
// compute covariant basis
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
// check for no-contact to contact condition
// need a better way to check if new node added is not presnt in the
// previous master elemets
auto & previous_master_elements = model.getPreviousMasterElements();
if (element.slave >= previous_master_elements.size())
return;
auto & previous_element = previous_master_elements[element.slave];
if (previous_element.type == _not_defined)
return;
// compute tangential traction using return map algorithm
auto & tangential_tractions = model.getTangentialTractions();
Vector<Real> tangential_traction(
tangential_tractions.begin(surface_dimension)[element.slave]);
this->computeTangentialTraction(element, covariant_basis,
tangential_traction);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> shape_matric(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
ResolutionUtils::computeShapeFunctionMatric(element, projection,
shape_matric);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
auto & nodal_area = const_cast<Array<Real> &>(model.getNodalArea());
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent_alpha = std::get<1>(values1);
for (auto && values2 : enumerate(tangential_traction)) {
auto & beta = std::get<0>(values2);
auto & traction_beta = std::get<1>(values2);
Vector<Real> tmp(force.size());
tmp.mul<true>(shape_matric, tangent_alpha, traction_beta);
tmp *=
contravariant_metric_tensor(alpha, beta) * nodal_area[element.slave];
force += tmp;
}
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTangentialTraction(
const ContactElement & element, const Matrix<Real> & covariant_basis,
Vector<Real> & traction_tangential) {
UInt surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
// Return map algorithm is employed
// compute trial traction
Vector<Real> traction_trial(surface_dimension);
this->computeTrialTangentialTraction(element, covariant_basis,
traction_trial);
// compute norm of trial traction
Real traction_trial_norm = 0;
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto i : arange(surface_dimension)) {
for (auto j : arange(surface_dimension)) {
traction_trial_norm += traction_trial[i] * traction_trial[j] *
contravariant_metric_tensor(i, j);
}
}
traction_trial_norm = sqrt(traction_trial_norm);
// check stick or slip condition
auto & contact_state = model.getContactState();
auto & state = contact_state.begin()[element.slave];
Real p_n = computeNormalTraction(gap);
- bool stick = (traction_trial_norm <= mu * p_n) ? true : false;
+ bool stick = (traction_trial_norm <= mu * p_n);
if (stick) {
state = ContactState::_stick;
computeStickTangentialTraction(element, traction_trial,
traction_tangential);
} else {
state = ContactState::_slip;
computeSlipTangentialTraction(element, covariant_basis, traction_trial,
traction_tangential);
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTrialTangentialTraction(
const ContactElement & element, const Matrix<Real> & covariant_basis,
Vector<Real> & traction) {
UInt surface_dimension = spatial_dimension - 1;
auto & projections = model.getProjections();
Vector<Real> current_projection(
projections.begin(surface_dimension)[element.slave]);
auto & previous_projections = model.getPreviousProjections();
Vector<Real> previous_projection(
previous_projections.begin(surface_dimension)[element.slave]);
// method from Laursen et. al.
/*auto covariant_metric_tensor =
GeometryUtils::covariantMetricTensor(covariant_basis); auto
increment_projection = current_projection - previous_projection;
traction.mul<false>(covariant_metric_tensor, increment_projection, epsilon_t);
auto & previous_tangential_tractions = model.getPreviousTangentialTractions();
Vector<Real>
previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]);
traction = previous_traction + traction;*/
// method from Schweizerhof
auto covariant_metric_tensor =
GeometryUtils::covariantMetricTensor(covariant_basis);
auto & previous_tangential_tractions = model.getPreviousTangentialTractions();
Vector<Real> previous_traction(
previous_tangential_tractions.begin(surface_dimension)[element.slave]);
auto & previous_tangents = model.getPreviousTangents();
Matrix<Real> previous_covariant_basis(previous_tangents.begin(
surface_dimension, spatial_dimension)[element.slave]);
auto previous_contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(previous_covariant_basis);
auto current_tangent = covariant_basis.transpose();
auto previous_tangent = previous_covariant_basis.transpose();
for (auto alpha : arange(surface_dimension)) {
Vector<Real> tangent_alpha(current_tangent(alpha));
for (auto gamma : arange(surface_dimension)) {
for (auto beta : arange(surface_dimension)) {
Vector<Real> tangent_beta(previous_tangent(beta));
auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha);
traction[alpha] += previous_traction[gamma] *
previous_contravariant_metric_tensor(gamma, beta) *
t_alpha_t_beta;
}
}
}
auto & previous_master_elements = model.getPreviousMasterElements();
auto & previous_element = previous_master_elements[element.slave];
Vector<Real> previous_real_projection(spatial_dimension);
GeometryUtils::realProjection(
model.getMesh(), model.getContactDetector().getPositions(),
previous_element, previous_projection, previous_real_projection);
Vector<Real> current_real_projection(spatial_dimension);
GeometryUtils::realProjection(
model.getMesh(), model.getContactDetector().getPositions(),
element.master, current_projection, current_real_projection);
auto increment_real = current_real_projection - previous_real_projection;
Vector<Real> increment_xi(surface_dimension);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
// increment in natural coordinate
for (auto beta : arange(surface_dimension)) {
for (auto gamma : arange(surface_dimension)) {
auto temp = increment_real.dot(current_tangent(gamma));
temp *= contravariant_metric_tensor(beta, gamma);
increment_xi[beta] += temp;
}
}
Vector<Real> temp(surface_dimension);
temp.mul<false>(covariant_metric_tensor, increment_xi, epsilon_t);
traction -= temp;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeStickTangentialTraction(
const ContactElement & /*element*/, Vector<Real> & traction_trial,
Vector<Real> & traction_tangential) {
traction_tangential = traction_trial;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeSlipTangentialTraction(
const ContactElement & element, const Matrix<Real> & covariant_basis,
Vector<Real> & traction_trial, Vector<Real> & traction_tangential) {
UInt surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
// compute norm of trial traction
Real traction_trial_norm = 0;
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto alpha : arange(surface_dimension)) {
for (auto beta : arange(surface_dimension)) {
traction_trial_norm += traction_trial[alpha] * traction_trial[beta] *
contravariant_metric_tensor(alpha, beta);
}
}
traction_trial_norm = sqrt(traction_trial_norm);
auto slip_direction = traction_trial;
slip_direction /= traction_trial_norm;
Real p_n = computeNormalTraction(gap);
traction_tangential = slip_direction;
traction_tangential *= mu * p_n;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeNormalModuli(const ContactElement & element,
Matrix<Real> & stiffness) {
auto surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
Real gap(gaps.begin()[element.slave]);
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
// method from Schweizerhof and A. Konyukhov, K. Schweizerhof
// DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3
// construct A matrix
const ElementType & type = element.master.type;
auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> A(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
A(j, i * spatial_dimension + j) = 1;
continue;
}
A(j, i * spatial_dimension + j) = -shapes[i - 1];
}
}
// construct the main part of normal matrix
Matrix<Real> k_main(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> n_outer_n(spatial_dimension, spatial_dimension);
Matrix<Real> mat_n(normal.storage(), normal.size(), 1.);
n_outer_n.mul<false, true>(mat_n, mat_n);
Matrix<Real> tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(n_outer_n, A);
k_main.mul<true, false>(A, tmp);
k_main *= epsilon_n * heaviside(gap) * nodal_area;
// construct the rotational part of the normal matrix
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
// computing shape derivatives
auto && shape_derivatives =
ElementClassHelper<_ek_regular>::getDNDS(projection, type);
// consists of 2 rotational parts
Matrix<Real> k_rot1(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_rot2(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
auto construct_Aj = [&](auto && dnds) {
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
Aj(j, i * spatial_dimension + j) = 0;
continue;
}
Aj(j, i * spatial_dimension + j) = dnds(i - 1);
}
}
};
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent = std::get<1>(values1);
Matrix<Real> n_outer_t(spatial_dimension, spatial_dimension);
Matrix<Real> mat_t(tangent.storage(), tangent.size(), 1.);
n_outer_t.mul<false, true>(mat_n, mat_t);
Matrix<Real> t_outer_n(spatial_dimension, spatial_dimension);
t_outer_n.mul<false, true>(mat_t, mat_n);
for (auto && values2 : enumerate(shape_derivatives.transpose())) {
auto & beta = std::get<0>(values2);
auto & dnds = std::get<1>(values2);
// construct Aj from shape function wrt to jth natural
// coordinate
construct_Aj(dnds);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(n_outer_t, A);
tmp1.mul<true, false>(Aj, tmp);
tmp1 *= contravariant_metric_tensor(alpha, beta);
k_rot1 += tmp1;
tmp.mul<false, false>(t_outer_n, Aj);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, beta);
k_rot2 += tmp1;
}
}
k_rot1 *= -epsilon_n * heaviside(gap) * gap * nodal_area;
k_rot2 *= -epsilon_n * heaviside(gap) * gap * nodal_area;
stiffness += k_main + k_rot1 + k_rot2;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeTangentialModuli(const ContactElement & element,
Matrix<Real> & stiffness) {
-
- if (mu == 0)
+ if (mu == 0) {
return;
+ }
stiffness.zero();
auto & contact_state = model.getContactState();
- UInt state = contact_state.begin()[element.slave];
+ auto state = contact_state.begin()[element.slave];
switch (state) {
case ContactState::_stick: {
computeStickModuli(element, stiffness);
break;
}
case ContactState::_slip: {
computeSlipModuli(element, stiffness);
break;
}
default:
break;
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeStickModuli(const ContactElement & element,
Matrix<Real> & stiffness) {
auto surface_dimension = spatial_dimension - 1;
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
// method from Schweizerhof and A. Konyukhov, K. Schweizerhof
// DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3
// construct A matrix
const ElementType & type = element.master.type;
auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> A(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
A(j, i * spatial_dimension + j) = 1;
continue;
}
A(j, i * spatial_dimension + j) = -shapes[i - 1];
}
}
// computing shape derivatives
auto && shape_derivatives =
ElementClassHelper<_ek_regular>::getDNDS(projection, type);
Matrix<Real> Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
auto construct_Aj = [&](auto && dnds) {
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
Aj(j, i * spatial_dimension + j) = 0;
continue;
}
Aj(j, i * spatial_dimension + j) = dnds(i - 1);
}
}
};
// tangents should have been calculated in normal modulii
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
// construct 1st part of the stick modulii
Matrix<Real> k_main(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent_alpha = std::get<1>(values1);
Matrix<Real> t_outer_t(spatial_dimension, spatial_dimension);
Matrix<Real> mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.);
for (auto && values2 : enumerate(covariant_basis.transpose())) {
auto & beta = std::get<0>(values2);
auto & tangent_beta = std::get<1>(values2);
Matrix<Real> mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.);
t_outer_t.mul<false, true>(mat_t_alpha, mat_t_beta);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(t_outer_t, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, beta);
k_main += tmp1;
}
}
k_main *= -epsilon_t;
// construct 2nd part of the stick modulii
auto & tangential_tractions = model.getTangentialTractions();
Vector<Real> tangential_traction(
tangential_tractions.begin(surface_dimension)[element.slave]);
Matrix<Real> k_second(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto alpha : arange(surface_dimension)) {
Matrix<Real> k_sum(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto && values1 : enumerate(shape_derivatives.transpose())) {
auto & beta = std::get<0>(values1);
auto & dnds = std::get<1>(values1);
// construct Aj from shape function wrt to jth natural
// coordinate
construct_Aj(dnds);
for (auto && values2 : enumerate(covariant_basis.transpose())) {
auto & gamma = std::get<0>(values2);
auto & tangent_gamma = std::get<1>(values2);
Matrix<Real> t_outer_t(spatial_dimension, spatial_dimension);
Matrix<Real> mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(),
1.);
for (auto && values3 : enumerate(covariant_basis.transpose())) {
auto & theta = std::get<0>(values3);
auto & tangent_theta = std::get<1>(values3);
Matrix<Real> mat_t_theta(tangent_theta.storage(),
tangent_theta.size(), 1.);
t_outer_t.mul<false, true>(mat_t_gamma, mat_t_theta);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(t_outer_t, Aj);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, theta) *
contravariant_metric_tensor(beta, gamma);
Matrix<Real> tmp2(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp3(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp2.mul<false, false>(t_outer_t, A);
tmp3.mul<true, false>(Aj, tmp2);
tmp3 *= contravariant_metric_tensor(alpha, gamma) *
contravariant_metric_tensor(beta, theta);
k_sum += tmp1 + tmp3;
}
}
}
k_second += tangential_traction[alpha] * k_sum;
}
stiffness += k_main * nodal_area - k_second * nodal_area;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::computeSlipModuli(const ContactElement & element,
Matrix<Real> & stiffness) {
auto surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
Real gap(gaps.begin()[element.slave]);
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
// compute normal traction
Real p_n = computeNormalTraction(gap);
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
// restructure normal as a matrix for an outer product
Matrix<Real> mat_n(normal.storage(), normal.size(), 1.);
// method from Schweizerhof and A. Konyukhov, K. Schweizerhof
// DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3
// construct A matrix
const ElementType & type = element.master.type;
auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> A(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
A(j, i * spatial_dimension + j) = 1;
continue;
}
A(j, i * spatial_dimension + j) = -shapes[i - 1];
}
}
// computing shape derivatives
auto && shape_derivatives =
ElementClassHelper<_ek_regular>::getDNDS(projection, type);
Matrix<Real> Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
auto construct_Aj = [&](auto && dnds) {
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
Aj(j, i * spatial_dimension + j) = 0;
continue;
}
Aj(j, i * spatial_dimension + j) = dnds(i - 1);
}
}
};
// tangents should have been calculated in normal modulii
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
auto & tangential_tractions = model.getTangentialTractions();
Vector<Real> tangential_traction(
tangential_tractions.begin(surface_dimension)[element.slave]);
// compute norm of trial traction
Real traction_norm = 0;
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto i : arange(surface_dimension)) {
for (auto j : arange(surface_dimension)) {
traction_norm += tangential_traction[i] * tangential_traction[j] *
contravariant_metric_tensor(i, j);
}
}
traction_norm = sqrt(traction_norm);
// construct four parts of stick modulii (eq 107,107a-c)
Matrix<Real> k_first(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_second(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_third(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_fourth(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent_alpha = std::get<1>(values1);
Matrix<Real> mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.);
Matrix<Real> t_outer_n(spatial_dimension, spatial_dimension);
Matrix<Real> t_outer_t(spatial_dimension, spatial_dimension);
for (auto && values2 :
zip(arange(surface_dimension), covariant_basis.transpose(),
shape_derivatives.transpose())) {
auto & beta = std::get<0>(values2);
auto & tangent_beta = std::get<1>(values2);
auto & dnds = std::get<2>(values2);
// construct Aj from shape function wrt to jth natural
// coordinate
construct_Aj(dnds);
// eq 107
Matrix<Real> mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.);
t_outer_n.mul<false, true>(mat_t_beta, mat_n);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(t_outer_n, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= epsilon_n * mu * tangential_traction[alpha] *
contravariant_metric_tensor(alpha, beta);
tmp1 /= traction_norm;
k_first += tmp1 * nodal_area;
// eq 107a
t_outer_t.mul<false, true>(mat_t_alpha, mat_t_beta);
tmp.mul<false, false>(t_outer_t, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta);
tmp1 /= traction_norm;
k_second += tmp1 * nodal_area;
for (auto && values3 : enumerate(covariant_basis.transpose())) {
auto & gamma = std::get<0>(values3);
auto & tangent_gamma = std::get<1>(values3);
Matrix<Real> mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(),
1.);
for (auto && values4 : enumerate(covariant_basis.transpose())) {
auto & theta = std::get<0>(values4);
auto & tangent_theta = std::get<1>(values4);
Matrix<Real> mat_t_theta(tangent_theta.storage(),
tangent_theta.size(), 1.);
t_outer_t.mul<false, true>(mat_t_gamma, mat_t_theta);
// eq 107b
tmp.mul<false, false>(t_outer_t, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= epsilon_t * mu * p_n * tangential_traction[alpha] *
tangential_traction[beta];
tmp1 *= contravariant_metric_tensor(alpha, gamma) *
contravariant_metric_tensor(beta, theta);
tmp1 /= pow(traction_norm, 3);
k_third += tmp1 * nodal_area;
// eq 107c
tmp.mul<false, false>(t_outer_t, Aj);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, theta) *
contravariant_metric_tensor(beta, gamma);
tmp1 *= mu * p_n * tangential_traction[alpha];
tmp1 /= traction_norm;
Matrix<Real> tmp2(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp3(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp2.mul<false, false>(t_outer_t, A);
tmp3.mul<true, false>(Aj, tmp2);
tmp3 *= contravariant_metric_tensor(alpha, gamma) *
contravariant_metric_tensor(beta, theta);
tmp3 *= mu * p_n * tangential_traction[alpha];
tmp3 /= traction_norm;
k_fourth += (tmp1 + tmp3) * nodal_area;
}
}
}
}
stiffness += k_third + k_fourth - k_first - k_second;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::beforeSolveStep() {}
/* -------------------------------------------------------------------------- */
void ResolutionPenalty::afterSolveStep(__attribute__((unused)) bool converged) {}
INSTANTIATE_RESOLUTION(penalty_linear, ResolutionPenalty);
} // namespace akantu
diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.hh b/src/model/contact_mechanics/resolutions/resolution_penalty.hh
index 4cae174f4..61f7ba197 100644
--- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh
+++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh
@@ -1,126 +1,122 @@
/**
* @file contact_resolution_penalty.hh
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Mon Jan 29 2018
*
* @brief Linear Penalty Resolution for Contact Mechanics Model
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "resolution.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_RESOLUTION_PENALTY_HH__
#define __AKANTU_RESOLUTION_PENALTY_HH__
namespace akantu {
class ResolutionPenalty : public Resolution {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ResolutionPenalty(ContactMechanicsModel & model, const ID & id = "");
~ResolutionPenalty() override = default;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// initialize the resolution
void initialize();
-
/* ------------------------------------------------------------------------ */
/* Methods for stiffness computation */
/* ------------------------------------------------------------------------ */
protected:
-
/// local computaion of stiffness matrix due to stick state
void computeStickModuli(const ContactElement &, Matrix<Real> &);
- /// local computation of stiffness matrix due to slip state
+ /// local computation of stiffness matrix due to slip state
void computeSlipModuli(const ContactElement &, Matrix<Real> &);
public:
/// local computation of tangent moduli due to normal traction
void computeNormalModuli(const ContactElement &, Matrix<Real> &) override;
-
+
/// local computation of tangent moduli due to tangential traction
void computeTangentialModuli(const ContactElement &, Matrix<Real> &) override;
-
+
/* ------------------------------------------------------------------------ */
/* Methods for force computation */
/* ------------------------------------------------------------------------ */
public:
/// local computation of normal force due to normal contact
void computeNormalForce(const ContactElement &, Vector<Real> &) override;
-
- /// local computation of tangential force due to frictional traction
+
+ /// local computation of tangential force due to frictional traction
void computeTangentialForce(const ContactElement &, Vector<Real> &) override;
protected:
/// local computation of normal traction due to penetration
Real computeNormalTraction(Real &);
-
+
/// local computation of trial tangential traction due to friction
- void computeTrialTangentialTraction(const ContactElement &, const Matrix<Real> &,
- Vector<Real> &);
+ void computeTrialTangentialTraction(const ContactElement &,
+ const Matrix<Real> &, Vector<Real> &);
- /// local computation of tangential traction due to stick
+ /// local computation of tangential traction due to stick
void computeStickTangentialTraction(const ContactElement &, Vector<Real> &,
- Vector<Real> &);
+ Vector<Real> &);
/// local computation of tangential traction due to slip
- void computeSlipTangentialTraction(const ContactElement &, const Matrix<Real> &,
- Vector<Real> &, Vector<Real> &);
+ void computeSlipTangentialTraction(const ContactElement &,
+ const Matrix<Real> &, Vector<Real> &,
+ Vector<Real> &);
/// local computation of tangential traction due to friction
void computeTangentialTraction(const ContactElement &, const Matrix<Real> &,
- Vector<Real> &);
+ Vector<Real> &);
public:
-
void beforeSolveStep() override;
void afterSolveStep(bool converged = true) override;
-
+
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// penalty parameter for normal traction
Real epsilon_n;
/// penalty parameter for tangential traction
Real epsilon_t;
};
-
-
-} // akantu
+} // namespace akantu
#endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */
diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc
index f0fe8e7e4..615d71fab 100644
--- a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc
+++ b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc
@@ -1,832 +1,832 @@
/**
* @file resolution_penalty_quadratic.cc
*
* @author Mohit Pundir <mohit.pundir@epfl.ch>
*
* @date creation: Sun Aug 2 2020
* @date last modification: Sun Aug 2 2020
*
* @brief Specialization of the resolution class for the quadratic penalty
* method
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "resolution_penalty_quadratic.hh"
#include "element_class_helper.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
ResolutionPenaltyQuadratic::ResolutionPenaltyQuadratic(
ContactMechanicsModel & model, const ID & id)
: Parent(model, id) {
AKANTU_DEBUG_IN();
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::initialize() {}
/* -------------------------------------------------------------------------- */
Real ResolutionPenaltyQuadratic::computeNormalTraction(Real & gap) {
return epsilon_n * (macaulay(gap) * macaulay(gap) + macaulay(gap));
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeNormalForce(
const ContactElement & element, Vector<Real> & force) {
force.zero();
auto & gaps = model.getGaps();
auto & projections = model.getProjections();
auto & normals = model.getNormals();
auto surface_dimension = spatial_dimension - 1;
Real gap(gaps.begin()[element.slave]);
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & nodal_area = const_cast<Array<Real> &>(model.getNodalArea());
// compute normal traction
Real p_n = computeNormalTraction(gap);
p_n *= nodal_area[element.slave];
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> shape_matric(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
ResolutionUtils::computeShapeFunctionMatric(element, projection,
shape_matric);
force.mul<true>(shape_matric, normal, p_n);
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeTangentialForce(
const ContactElement & element, Vector<Real> & force) {
if (mu == 0)
return;
force.zero();
UInt surface_dimension = spatial_dimension - 1;
// compute tangents
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
// check for no-contact to contact condition
// need a better way to check if new node added is not presnt in the
// previous master elemets
auto & previous_master_elements = model.getPreviousMasterElements();
if (element.slave >= previous_master_elements.size())
return;
auto & previous_element = previous_master_elements[element.slave];
if (previous_element.type == _not_defined)
return;
// compute tangential traction using return map algorithm
auto & tangential_tractions = model.getTangentialTractions();
Vector<Real> tangential_traction(
tangential_tractions.begin(surface_dimension)[element.slave]);
this->computeTangentialTraction(element, covariant_basis,
tangential_traction);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> shape_matric(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
ResolutionUtils::computeShapeFunctionMatric(element, projection,
shape_matric);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
auto & nodal_area = const_cast<Array<Real> &>(model.getNodalArea());
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent_alpha = std::get<1>(values1);
for (auto && values2 : enumerate(tangential_traction)) {
auto & beta = std::get<0>(values2);
auto & traction_beta = std::get<1>(values2);
Vector<Real> tmp(force.size());
tmp.mul<true>(shape_matric, tangent_alpha, traction_beta);
tmp *=
contravariant_metric_tensor(alpha, beta) * nodal_area[element.slave];
force += tmp;
}
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeTangentialTraction(
const ContactElement & element, const Matrix<Real> & covariant_basis,
Vector<Real> & traction_tangential) {
UInt surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
// Return map algorithm is employed
// compute trial traction
Vector<Real> traction_trial(surface_dimension);
this->computeTrialTangentialTraction(element, covariant_basis,
traction_trial);
// compute norm of trial traction
Real traction_trial_norm = 0;
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto i : arange(surface_dimension)) {
for (auto j : arange(surface_dimension)) {
traction_trial_norm += traction_trial[i] * traction_trial[j] *
contravariant_metric_tensor(i, j);
}
}
traction_trial_norm = sqrt(traction_trial_norm);
// check stick or slip condition
auto & contact_state = model.getContactState();
auto & state = contact_state.begin()[element.slave];
Real p_n = computeNormalTraction(gap);
bool stick = (traction_trial_norm <= mu * p_n) ? true : false;
if (stick) {
state = ContactState::_stick;
computeStickTangentialTraction(element, traction_trial,
traction_tangential);
} else {
state = ContactState::_slip;
computeSlipTangentialTraction(element, covariant_basis, traction_trial,
traction_tangential);
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeTrialTangentialTraction(
const ContactElement & element, const Matrix<Real> & covariant_basis,
Vector<Real> & traction) {
UInt surface_dimension = spatial_dimension - 1;
auto & projections = model.getProjections();
Vector<Real> current_projection(
projections.begin(surface_dimension)[element.slave]);
auto & previous_projections = model.getPreviousProjections();
Vector<Real> previous_projection(
previous_projections.begin(surface_dimension)[element.slave]);
// method from Laursen et. al.
/*auto covariant_metric_tensor =
GeometryUtils::covariantMetricTensor(covariant_basis); auto
increment_projection = current_projection - previous_projection;
traction.mul<false>(covariant_metric_tensor, increment_projection, epsilon_t);
auto & previous_tangential_tractions = model.getPreviousTangentialTractions();
Vector<Real>
previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]);
traction = previous_traction + traction;*/
// method from Schweizerhof
auto covariant_metric_tensor =
GeometryUtils::covariantMetricTensor(covariant_basis);
auto & previous_tangential_tractions = model.getPreviousTangentialTractions();
Vector<Real> previous_traction(
previous_tangential_tractions.begin(surface_dimension)[element.slave]);
auto & previous_tangents = model.getPreviousTangents();
Matrix<Real> previous_covariant_basis(previous_tangents.begin(
surface_dimension, spatial_dimension)[element.slave]);
auto previous_contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(previous_covariant_basis);
auto current_tangent = covariant_basis.transpose();
auto previous_tangent = previous_covariant_basis.transpose();
for (auto alpha : arange(surface_dimension)) {
Vector<Real> tangent_alpha(current_tangent(alpha));
for (auto gamma : arange(surface_dimension)) {
for (auto beta : arange(surface_dimension)) {
Vector<Real> tangent_beta(previous_tangent(beta));
auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha);
traction[alpha] += previous_traction[gamma] *
previous_contravariant_metric_tensor(gamma, beta) *
t_alpha_t_beta;
}
}
}
auto & previous_master_elements = model.getPreviousMasterElements();
auto & previous_element = previous_master_elements[element.slave];
Vector<Real> previous_real_projection(spatial_dimension);
GeometryUtils::realProjection(
model.getMesh(), model.getContactDetector().getPositions(),
previous_element, previous_projection, previous_real_projection);
Vector<Real> current_real_projection(spatial_dimension);
GeometryUtils::realProjection(
model.getMesh(), model.getContactDetector().getPositions(),
element.master, current_projection, current_real_projection);
auto increment_real = current_real_projection - previous_real_projection;
Vector<Real> increment_xi(surface_dimension);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
// increment in natural coordinate
for (auto beta : arange(surface_dimension)) {
for (auto gamma : arange(surface_dimension)) {
auto temp = increment_real.dot(current_tangent(gamma));
temp *= contravariant_metric_tensor(beta, gamma);
increment_xi[beta] += temp;
}
}
Vector<Real> temp(surface_dimension);
temp.mul<false>(covariant_metric_tensor, increment_xi, epsilon_t);
traction -= temp;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeStickTangentialTraction(
const ContactElement & /*element*/, Vector<Real> & traction_trial,
Vector<Real> & traction_tangential) {
traction_tangential = traction_trial;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeSlipTangentialTraction(
const ContactElement & element, const Matrix<Real> & covariant_basis,
Vector<Real> & traction_trial, Vector<Real> & traction_tangential) {
UInt surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
auto & gap = gaps.begin()[element.slave];
// compute norm of trial traction
Real traction_trial_norm = 0;
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto i : arange(surface_dimension)) {
for (auto j : arange(surface_dimension)) {
traction_trial_norm += traction_trial[i] * traction_trial[j] *
contravariant_metric_tensor(i, j);
}
}
traction_trial_norm = sqrt(traction_trial_norm);
auto slip_direction = traction_trial;
slip_direction /= traction_trial_norm;
Real p_n = computeNormalTraction(gap);
traction_tangential = slip_direction;
traction_tangential *= mu * p_n;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeNormalModuli(
const ContactElement & element, Matrix<Real> & stiffness) {
auto surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
Real gap(gaps.begin()[element.slave]);
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
// method from Schweizerhof and A. Konyukhov, K. Schweizerhof
// DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3
// construct A matrix
const ElementType & type = element.master.type;
auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> A(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
A(j, i * spatial_dimension + j) = 1;
continue;
}
A(j, i * spatial_dimension + j) = -shapes[i - 1];
}
}
// construct the main part of normal matrix
Matrix<Real> k_main(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> n_outer_n(spatial_dimension, spatial_dimension);
Matrix<Real> mat_n(normal.storage(), normal.size(), 1.);
n_outer_n.mul<false, true>(mat_n, mat_n);
Matrix<Real> tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(n_outer_n, A);
k_main.mul<true, false>(A, tmp);
k_main *= epsilon_n * heaviside(gap) * (2 * gap + 1) * nodal_area;
// construct the rotational part of the normal matrix
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
// computing shape derivatives
auto && shape_derivatives =
ElementClassHelper<_ek_regular>::getDNDS(projection, type);
// consists of 2 rotational parts
Matrix<Real> k_rot1(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_rot2(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
auto construct_Aj = [&](auto && dnds) {
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
Aj(j, i * spatial_dimension + j) = 0;
continue;
}
Aj(j, i * spatial_dimension + j) = dnds(i - 1);
}
}
};
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent = std::get<1>(values1);
Matrix<Real> n_outer_t(spatial_dimension, spatial_dimension);
Matrix<Real> mat_t(tangent.storage(), tangent.size(), 1.);
n_outer_t.mul<false, true>(mat_n, mat_t);
Matrix<Real> t_outer_n(spatial_dimension, spatial_dimension);
t_outer_n.mul<false, true>(mat_t, mat_n);
for (auto && values2 : enumerate(shape_derivatives.transpose())) {
auto & beta = std::get<0>(values2);
auto & dnds = std::get<1>(values2);
// construct Aj from shape function wrt to jth natural
// coordinate
construct_Aj(dnds);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(n_outer_t, A);
tmp1.mul<true, false>(Aj, tmp);
tmp1 *= contravariant_metric_tensor(alpha, beta);
k_rot1 += tmp1;
tmp.mul<false, false>(t_outer_n, Aj);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, beta);
k_rot2 += tmp1;
}
}
k_rot1 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area;
k_rot2 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area;
stiffness += k_main + k_rot1 + k_rot2;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeTangentialModuli(
const ContactElement & element, Matrix<Real> & stiffness) {
-
- if (mu == 0)
+ if (mu == 0) {
return;
+ }
stiffness.zero();
auto & contact_state = model.getContactState();
- UInt state = contact_state.begin()[element.slave];
+ auto state = contact_state.begin()[element.slave];
switch (state) {
case ContactState::_stick: {
computeStickModuli(element, stiffness);
break;
}
case ContactState::_slip: {
computeSlipModuli(element, stiffness);
break;
}
default:
break;
}
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeStickModuli(
const ContactElement & element, Matrix<Real> & stiffness) {
auto surface_dimension = spatial_dimension - 1;
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
// method from Schweizerhof and A. Konyukhov, K. Schweizerhof
// DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3
// construct A matrix
const ElementType & type = element.master.type;
auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> A(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
A(j, i * spatial_dimension + j) = 1;
continue;
}
A(j, i * spatial_dimension + j) = -shapes[i - 1];
}
}
// computing shape derivatives
auto && shape_derivatives =
ElementClassHelper<_ek_regular>::getDNDS(projection, type);
Matrix<Real> Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
auto construct_Aj = [&](auto && dnds) {
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
Aj(j, i * spatial_dimension + j) = 0;
continue;
}
Aj(j, i * spatial_dimension + j) = dnds(i - 1);
}
}
};
// tangents should have been calculated in normal modulii
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
// construct 1st part of the stick modulii
Matrix<Real> k_main(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent_alpha = std::get<1>(values1);
Matrix<Real> t_outer_t(spatial_dimension, spatial_dimension);
Matrix<Real> mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.);
for (auto && values2 : enumerate(covariant_basis.transpose())) {
auto & beta = std::get<0>(values2);
auto & tangent_beta = std::get<1>(values2);
Matrix<Real> mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.);
t_outer_t.mul<false, true>(mat_t_alpha, mat_t_beta);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(t_outer_t, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, beta);
k_main += tmp1;
}
}
k_main *= -epsilon_t;
// construct 2nd part of the stick modulii
auto & tangential_tractions = model.getTangentialTractions();
Vector<Real> tangential_traction(
tangential_tractions.begin(surface_dimension)[element.slave]);
Matrix<Real> k_second(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto alpha : arange(surface_dimension)) {
Matrix<Real> k_sum(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto && values1 : enumerate(shape_derivatives.transpose())) {
auto & beta = std::get<0>(values1);
auto & dnds = std::get<1>(values1);
// construct Aj from shape function wrt to jth natural
// coordinate
construct_Aj(dnds);
for (auto && values2 : enumerate(covariant_basis.transpose())) {
auto & gamma = std::get<0>(values2);
auto & tangent_gamma = std::get<1>(values2);
Matrix<Real> t_outer_t(spatial_dimension, spatial_dimension);
Matrix<Real> mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(),
1.);
for (auto && values3 : enumerate(covariant_basis.transpose())) {
auto & theta = std::get<0>(values3);
auto & tangent_theta = std::get<1>(values3);
Matrix<Real> mat_t_theta(tangent_theta.storage(),
tangent_theta.size(), 1.);
t_outer_t.mul<false, true>(mat_t_gamma, mat_t_theta);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(t_outer_t, Aj);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, theta) *
contravariant_metric_tensor(beta, gamma);
Matrix<Real> tmp2(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp3(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp2.mul<false, false>(t_outer_t, A);
tmp3.mul<true, false>(Aj, tmp2);
tmp3 *= contravariant_metric_tensor(alpha, gamma) *
contravariant_metric_tensor(beta, theta);
k_sum += tmp1 + tmp3;
}
}
}
k_second += tangential_traction[alpha] * k_sum;
}
stiffness += k_main * nodal_area - k_second * nodal_area;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::computeSlipModuli(
const ContactElement & element, Matrix<Real> & stiffness) {
auto surface_dimension = spatial_dimension - 1;
auto & gaps = model.getGaps();
Real gap(gaps.begin()[element.slave]);
auto & nodal_areas = model.getNodalArea();
auto & nodal_area = nodal_areas.begin()[element.slave];
// compute normal traction
Real p_n = computeNormalTraction(gap);
auto & projections = model.getProjections();
Vector<Real> projection(projections.begin(surface_dimension)[element.slave]);
auto & normals = model.getNormals();
Vector<Real> normal(normals.begin(spatial_dimension)[element.slave]);
// restructure normal as a matrix for an outer product
Matrix<Real> mat_n(normal.storage(), normal.size(), 1.);
// method from Schweizerhof and A. Konyukhov, K. Schweizerhof
// DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3
// construct A matrix
const ElementType & type = element.master.type;
auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type);
UInt nb_nodes_per_contact = element.getNbNodes();
Matrix<Real> A(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
A(j, i * spatial_dimension + j) = 1;
continue;
}
A(j, i * spatial_dimension + j) = -shapes[i - 1];
}
}
// computing shape derivatives
auto && shape_derivatives =
ElementClassHelper<_ek_regular>::getDNDS(projection, type);
Matrix<Real> Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact);
auto construct_Aj = [&](auto && dnds) {
for (auto i : arange(nb_nodes_per_contact)) {
for (auto j : arange(spatial_dimension)) {
if (i == 0) {
Aj(j, i * spatial_dimension + j) = 0;
continue;
}
Aj(j, i * spatial_dimension + j) = dnds(i - 1);
}
}
};
// tangents should have been calculated in normal modulii
auto & tangents = model.getTangents();
Matrix<Real> covariant_basis(
tangents.begin(surface_dimension, spatial_dimension)[element.slave]);
auto & tangential_tractions = model.getTangentialTractions();
Vector<Real> tangential_traction(
tangential_tractions.begin(surface_dimension)[element.slave]);
// compute norm of trial traction
Real traction_norm = 0;
auto contravariant_metric_tensor =
GeometryUtils::contravariantMetricTensor(covariant_basis);
for (auto i : arange(surface_dimension)) {
for (auto j : arange(surface_dimension)) {
traction_norm += tangential_traction[i] * tangential_traction[j] *
contravariant_metric_tensor(i, j);
}
}
traction_norm = sqrt(traction_norm);
// construct four parts of stick modulii (eq 107,107a-c)
Matrix<Real> k_first(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_second(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_third(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
Matrix<Real> k_fourth(nb_nodes_per_contact * spatial_dimension,
nb_nodes_per_contact * spatial_dimension);
for (auto && values1 : enumerate(covariant_basis.transpose())) {
auto & alpha = std::get<0>(values1);
auto & tangent_alpha = std::get<1>(values1);
Matrix<Real> mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.);
Matrix<Real> t_outer_n(spatial_dimension, spatial_dimension);
Matrix<Real> t_outer_t(spatial_dimension, spatial_dimension);
for (auto && values2 :
zip(arange(surface_dimension), covariant_basis.transpose(),
shape_derivatives.transpose())) {
auto & beta = std::get<0>(values2);
auto & tangent_beta = std::get<1>(values2);
auto & dnds = std::get<2>(values2);
// construct Aj from shape function wrt to jth natural
// coordinate
construct_Aj(dnds);
// eq 107
Matrix<Real> mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.);
t_outer_n.mul<false, true>(mat_t_beta, mat_n);
Matrix<Real> tmp(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp1(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp.mul<false, false>(t_outer_n, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= epsilon_n * mu * tangential_traction[alpha] *
contravariant_metric_tensor(alpha, beta);
tmp1 /= traction_norm;
k_first += tmp1 * nodal_area;
// eq 107a
t_outer_t.mul<false, true>(mat_t_alpha, mat_t_beta);
tmp.mul<false, false>(t_outer_t, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta);
tmp1 /= traction_norm;
k_second += tmp1 * nodal_area;
for (auto && values3 : enumerate(covariant_basis.transpose())) {
auto & gamma = std::get<0>(values3);
auto & tangent_gamma = std::get<1>(values3);
Matrix<Real> mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(),
1.);
for (auto && values4 : enumerate(covariant_basis.transpose())) {
auto & theta = std::get<0>(values4);
auto & tangent_theta = std::get<1>(values4);
Matrix<Real> mat_t_theta(tangent_theta.storage(),
tangent_theta.size(), 1.);
t_outer_t.mul<false, true>(mat_t_gamma, mat_t_theta);
// eq 107b
tmp.mul<false, false>(t_outer_t, A);
tmp1.mul<true, false>(A, tmp);
tmp1 *= epsilon_t * mu * p_n * tangential_traction[alpha] *
tangential_traction[beta];
tmp1 *= contravariant_metric_tensor(alpha, gamma) *
contravariant_metric_tensor(beta, theta);
tmp1 /= pow(traction_norm, 3);
k_third += tmp1 * nodal_area;
// eq 107c
tmp.mul<false, false>(t_outer_t, Aj);
tmp1.mul<true, false>(A, tmp);
tmp1 *= contravariant_metric_tensor(alpha, theta) *
contravariant_metric_tensor(beta, gamma);
tmp1 *= mu * p_n * tangential_traction[alpha];
tmp1 /= traction_norm;
Matrix<Real> tmp2(spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
Matrix<Real> tmp3(nb_nodes_per_contact * spatial_dimension,
spatial_dimension * nb_nodes_per_contact);
tmp2.mul<false, false>(t_outer_t, A);
tmp3.mul<true, false>(Aj, tmp2);
tmp3 *= contravariant_metric_tensor(alpha, gamma) *
contravariant_metric_tensor(beta, theta);
tmp3 *= mu * p_n * tangential_traction[alpha];
tmp3 /= traction_norm;
k_fourth += (tmp1 + tmp3) * nodal_area;
}
}
}
}
stiffness += k_third + k_fourth - k_first - k_second;
}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::beforeSolveStep() {}
/* -------------------------------------------------------------------------- */
void ResolutionPenaltyQuadratic::afterSolveStep(
__attribute__((unused)) bool converged) {}
INSTANTIATE_RESOLUTION(penalty_quadratic, ResolutionPenaltyQuadratic);
} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_io.cc b/src/model/solid_mechanics/solid_mechanics_model_io.cc
index b6fc9c8da..f6d6d5d50 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_io.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_io.cc
@@ -1,338 +1,337 @@
/**
* @file solid_mechanics_model_io.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Sun Jul 09 2017
* @date last modification: Sun Dec 03 2017
*
* @brief Dumpable part of the SolidMechnicsModel
*
*
* Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
#include "group_manager_inline_impl.hh"
#include "dumpable_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_element_partition.hh"
#include "dumper_elemental_field.hh"
#include "dumper_field.hh"
#include "dumper_homogenizing_field.hh"
#include "dumper_internal_material_field.hh"
#include "dumper_iohelper.hh"
#include "dumper_material_padders.hh"
#include "dumper_paraview.hh"
#endif
namespace akantu {
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::isInternal(const std::string & field_name,
ElementKind element_kind) {
/// check if at least one material contains field_id as an internal
for (auto & material : materials) {
bool is_internal = material->isInternal<Real>(field_name, element_kind);
if (is_internal) {
return true;
}
}
return false;
}
/* -------------------------------------------------------------------------- */
ElementTypeMap<UInt>
SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name,
ElementKind element_kind) {
-
if (!(this->isInternal(field_name, element_kind))) {
AKANTU_EXCEPTION("unknown internal " << field_name);
}
for (auto & material : materials) {
if (material->isInternal<Real>(field_name, element_kind)) {
return material->getInternalDataPerElem<Real>(field_name, element_kind);
}
}
return ElementTypeMap<UInt>();
}
/* -------------------------------------------------------------------------- */
ElementTypeMapArray<Real> &
SolidMechanicsModel::flattenInternal(const std::string & field_name,
ElementKind kind,
const GhostType ghost_type) {
auto key = std::make_pair(field_name, kind);
ElementTypeMapArray<Real> * internal_flat;
auto it = this->registered_internals.find(key);
if (it == this->registered_internals.end()) {
auto internal = std::make_unique<ElementTypeMapArray<Real>>(
field_name, this->id);
internal_flat = internal.get();
this->registered_internals[key] = std::move(internal);
} else {
internal_flat = it->second.get();
}
for (auto type :
mesh.elementTypes(Model::spatial_dimension, ghost_type, kind)) {
if (internal_flat->exists(type, ghost_type)) {
auto & internal = (*internal_flat)(type, ghost_type);
internal.resize(0);
}
}
for (auto & material : materials) {
if (material->isInternal<Real>(field_name, kind)) {
material->flattenInternal(field_name, *internal_flat, ghost_type, kind);
}
}
return *internal_flat;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::flattenAllRegisteredInternals(
ElementKind kind) {
ElementKind _kind;
ID _id;
for (auto & internal : this->registered_internals) {
std::tie(_id, _kind) = internal.first;
if (kind == _kind) {
this->flattenInternal(_id, kind);
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onDump() {
this->flattenAllRegisteredInternals(_ek_regular);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
std::shared_ptr<dumpers::Field> SolidMechanicsModel::createElementalField(
const std::string & field_name, const std::string & group_name,
bool padding_flag, UInt spatial_dimension,
ElementKind kind) {
std::shared_ptr<dumpers::Field> field;
if (field_name == "partitions") {
field = mesh.createElementalField<UInt, dumpers::ElementPartitionField>(
mesh.getConnectivities(), group_name, spatial_dimension, kind);
} else if (field_name == "material_index") {
field = mesh.createElementalField<UInt, Vector, dumpers::ElementalField>(
material_index, group_name, spatial_dimension, kind);
} else {
// this copy of field_name is used to compute derivated data such as
// strain and von mises stress that are based on grad_u and stress
std::string field_name_copy(field_name);
if (field_name == "strain" || field_name == "Green strain" ||
field_name == "principal strain" ||
field_name == "principal Green strain") {
field_name_copy = "grad_u";
} else if (field_name == "Von Mises stress") {
field_name_copy = "stress";
}
bool is_internal = this->isInternal(field_name_copy, kind);
if (is_internal) {
auto nb_data_per_elem =
this->getInternalDataPerElem(field_name_copy, kind);
auto & internal_flat = this->flattenInternal(field_name_copy, kind);
field = mesh.createElementalField<Real, dumpers::InternalMaterialField>(
internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem);
std::unique_ptr<dumpers::ComputeFunctorInterface> func;
if (field_name == "strain") {
func = std::make_unique<dumpers::ComputeStrain<false>>(*this);
} else if (field_name == "Von Mises stress") {
func = std::make_unique<dumpers::ComputeVonMisesStress>(*this);
} else if (field_name == "Green strain") {
func = std::make_unique<dumpers::ComputeStrain<true>>(*this);
} else if (field_name == "principal strain") {
func = std::make_unique<dumpers::ComputePrincipalStrain<false>>(*this);
} else if (field_name == "principal Green strain") {
func = std::make_unique<dumpers::ComputePrincipalStrain<true>>(*this);
}
if (func) {
field = dumpers::FieldComputeProxy::createFieldCompute(field,
std::move(func));
}
// treat the paddings
if (padding_flag) {
if (field_name == "stress") {
if (spatial_dimension == 2) {
auto foo = std::make_unique<dumpers::StressPadder<2>>(*this);
field = dumpers::FieldComputeProxy::createFieldCompute(
field, std::move(foo));
}
} else if (field_name == "strain" || field_name == "Green strain") {
if (spatial_dimension == 2) {
auto foo = std::make_unique<dumpers::StrainPadder<2>>(*this);
field = dumpers::FieldComputeProxy::createFieldCompute(
field, std::move(foo));
}
}
}
// homogenize the field
auto foo = dumpers::HomogenizerProxy::createHomogenizer(*field);
field =
dumpers::FieldComputeProxy::createFieldCompute(field, std::move(foo));
}
}
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string, Array<Real> *> real_nodal_fields;
real_nodal_fields["displacement"] = this->displacement.get();
real_nodal_fields["mass"] = this->mass.get();
real_nodal_fields["velocity"] = this->velocity.get();
real_nodal_fields["acceleration"] = this->acceleration.get();
real_nodal_fields["external_force"] = this->external_force.get();
real_nodal_fields["internal_force"] = this->internal_force.get();
real_nodal_fields["increment"] = this->displacement_increment.get();
if (field_name == "force") {
AKANTU_EXCEPTION("The 'force' field has been renamed in 'external_force'");
} else if (field_name == "residual") {
AKANTU_EXCEPTION(
"The 'residual' field has been replaced by 'internal_force'");
}
std::shared_ptr<dumpers::Field> field;
if (padding_flag) {
field = this->mesh.createNodalField(real_nodal_fields[field_name],
group_name, 3);
} else {
field =
this->mesh.createNodalField(real_nodal_fields[field_name], group_name);
}
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field> SolidMechanicsModel::createNodalFieldBool(
const std::string & field_name, const std::string & group_name,
__attribute__((unused)) bool padding_flag) {
std::map<std::string, Array<bool> *> uint_nodal_fields;
uint_nodal_fields["blocked_dofs"] = blocked_dofs.get();
std::shared_ptr<dumpers::Field> field;
field = mesh.createNodalField(uint_nodal_fields[field_name], group_name);
return field;
}
/* -------------------------------------------------------------------------- */
#else
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
SolidMechanicsModel::createElementalField(const std::string &,
const std::string &, bool,
const UInt &, ElementKind) {
return nullptr;
}
/* --------------------------------------------------------------------------
*/
std::shaed_ptr<dumpers::Field>
SolidMechanicsModel::createNodalFieldReal(const std::string &,
const std::string &, bool) {
return nullptr;
}
/* --------------------------------------------------------------------------
*/
std::shared_ptr<dumpers::Field>
SolidMechanicsModel::createNodalFieldBool(const std::string &,
const std::string &, bool) {
return nullptr;
}
#endif
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, Real time,
UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump() {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(Real time, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(time, step);
}
} // namespace akantu