diff --git a/packages/core.cmake b/packages/core.cmake
index 52ad62250..8cab12b7a 100644
--- a/packages/core.cmake
+++ b/packages/core.cmake
@@ -1,442 +1,443 @@
#===============================================================================
# @file core.cmake
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date creation: Mon Nov 21 2011
# @date last modification: Fri Apr 09 2021
#
# @brief package description for core
#
#
# @section LICENSE
#
# Copyright (©) 2010-2021 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/>.
#
#===============================================================================
package_declare(core NOT_OPTIONAL
DESCRIPTION "core package for Akantu"
FEATURES_PUBLIC cxx_strong_enums cxx_defaulted_functions
cxx_deleted_functions cxx_auto_type cxx_decltype_auto
FEATURES_PRIVATE cxx_lambdas cxx_nullptr cxx_range_for
cxx_delegating_constructors
DEPENDS INTERFACE akantu_iterators Boost
)
if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU" OR CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
package_set_compile_flags(core "-Wall -Wextra -pedantic")
else()
package_set_compile_flags(core "-Wall")
endif()
package_declare_sources(core
common/aka_array.cc
common/aka_array.hh
common/aka_array_filter.hh
common/aka_array_tmpl.hh
common/aka_array_printer.hh
common/aka_bbox.hh
common/aka_blas_lapack.hh
common/aka_circular_array.hh
common/aka_circular_array_inline_impl.hh
common/aka_common.cc
common/aka_common.hh
common/aka_common_inline_impl.hh
common/aka_csr.hh
common/aka_element_classes_info_inline_impl.hh
common/aka_enum_macros.hh
common/aka_error.cc
common/aka_error.hh
common/aka_event_handler_manager.hh
common/aka_extern.cc
common/aka_factory.hh
common/aka_fwd.hh
common/aka_grid_dynamic.hh
common/aka_math.cc
common/aka_math.hh
common/aka_math_tmpl.hh
common/aka_named_argument.hh
common/aka_random_generator.hh
common/aka_safe_enum.hh
common/aka_types.hh
common/aka_voigthelper.hh
common/aka_voigthelper_tmpl.hh
common/aka_voigthelper.cc
common/aka_warning.hh
common/aka_warning_restore.hh
fe_engine/element_class.hh
fe_engine/element_class_helper.hh
fe_engine/element_class_tmpl.hh
fe_engine/element_classes/element_class_hexahedron_8_inline_impl.hh
fe_engine/element_classes/element_class_hexahedron_20_inline_impl.hh
fe_engine/element_classes/element_class_pentahedron_6_inline_impl.hh
fe_engine/element_classes/element_class_pentahedron_15_inline_impl.hh
fe_engine/element_classes/element_class_point_1_inline_impl.hh
fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh
fe_engine/element_classes/element_class_quadrangle_8_inline_impl.hh
fe_engine/element_classes/element_class_segment_2_inline_impl.hh
fe_engine/element_classes/element_class_segment_3_inline_impl.hh
fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.hh
fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.hh
fe_engine/element_classes/element_class_triangle_3_inline_impl.hh
fe_engine/element_classes/element_class_triangle_6_inline_impl.hh
fe_engine/element_type_conversion.hh
fe_engine/fe_engine.cc
fe_engine/fe_engine.hh
fe_engine/fe_engine_inline_impl.hh
fe_engine/fe_engine_template.hh
fe_engine/fe_engine_template_tmpl_field.hh
fe_engine/fe_engine_template_tmpl.hh
fe_engine/geometrical_element_property.hh
fe_engine/geometrical_element_property.cc
fe_engine/gauss_integration.cc
fe_engine/gauss_integration_tmpl.hh
fe_engine/integrator.hh
fe_engine/integrator_gauss.hh
fe_engine/integrator_gauss_inline_impl.hh
fe_engine/interpolation_element_tmpl.hh
fe_engine/integration_point.hh
fe_engine/shape_functions.hh
fe_engine/shape_functions.cc
fe_engine/shape_functions_inline_impl.hh
fe_engine/shape_lagrange_base.cc
fe_engine/shape_lagrange_base.hh
fe_engine/shape_lagrange_base_inline_impl.hh
fe_engine/shape_lagrange.hh
fe_engine/shape_lagrange_inline_impl.hh
fe_engine/element.hh
io/dumper/dumpable.hh
io/dumper/dumpable.cc
io/dumper/dumpable_dummy.hh
io/dumper/dumpable_inline_impl.hh
io/dumper/dumper_field.hh
io/dumper/dumper_material_padders.hh
io/dumper/dumper_filtered_connectivity.hh
io/dumper/dumper_element_partition.hh
io/mesh_io.cc
io/mesh_io.hh
io/mesh_io/mesh_io_diana.cc
io/mesh_io/mesh_io_diana.hh
io/mesh_io/mesh_io_msh.cc
io/mesh_io/mesh_io_msh.hh
io/parser/algebraic_parser.hh
io/parser/input_file_parser.hh
io/parser/parsable.cc
io/parser/parsable.hh
io/parser/parser.cc
io/parser/parser_real.cc
io/parser/parser_random.cc
io/parser/parser_types.cc
io/parser/parser_input_files.cc
io/parser/parser.hh
io/parser/parser_tmpl.hh
io/parser/parser_grammar_tmpl.hh
io/parser/cppargparse/cppargparse.hh
io/parser/cppargparse/cppargparse.cc
io/parser/cppargparse/cppargparse_tmpl.hh
io/parser/parameter_registry.cc
io/parser/parameter_registry.hh
io/parser/parameter_registry_tmpl.hh
mesh/element_group.cc
mesh/element_group.hh
mesh/element_group_inline_impl.hh
mesh/element_type_map.cc
mesh/element_type_map.hh
mesh/element_type_map_tmpl.hh
mesh/element_type_map_filter.hh
mesh/group_manager.cc
mesh/group_manager.hh
mesh/group_manager_inline_impl.hh
mesh/mesh.cc
mesh/mesh.hh
mesh/mesh_periodic.cc
mesh/mesh_accessor.hh
mesh/mesh_events.hh
mesh/mesh_filter.hh
mesh/mesh_global_data_updater.hh
mesh/mesh_data.cc
mesh/mesh_data.hh
mesh/mesh_data_tmpl.hh
mesh/mesh_inline_impl.hh
mesh/node_group.cc
mesh/node_group.hh
mesh/node_group_inline_impl.hh
mesh/mesh_iterators.hh
mesh_utils/mesh_partition.cc
mesh_utils/mesh_partition.hh
mesh_utils/mesh_partition/mesh_partition_mesh_data.cc
mesh_utils/mesh_partition/mesh_partition_mesh_data.hh
mesh_utils/mesh_partition/mesh_partition_scotch.hh
mesh_utils/mesh_utils_pbc.cc
mesh_utils/mesh_utils.cc
mesh_utils/mesh_utils.hh
mesh_utils/mesh_utils_distribution.cc
mesh_utils/mesh_utils_distribution.hh
mesh_utils/mesh_utils.hh
mesh_utils/mesh_utils_inline_impl.hh
mesh_utils/global_ids_updater.hh
mesh_utils/global_ids_updater.cc
mesh_utils/global_ids_updater_inline_impl.hh
model/common/boundary_condition/boundary_condition.hh
model/common/boundary_condition/boundary_condition_functor.hh
model/common/boundary_condition/boundary_condition_functor_inline_impl.hh
model/common/boundary_condition/boundary_condition_tmpl.hh
model/common/constitutive_laws/constitutive_law_non_local_interface.hh
model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
model/common/constitutive_laws/constitutive_law.hh
+ model/common/constitutive_laws/constitutive_law_selector.hh
model/common/constitutive_laws/constitutive_law_tmpl.hh
model/common/constitutive_laws/constitutive_laws_handler.hh
model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh
model/common/constitutive_laws/internal_field.hh
model/common/constitutive_laws/internal_field_tmpl.hh
model/common/constitutive_laws/random_internal_field.hh
model/common/constitutive_laws/random_internal_field_tmpl.hh
model/common/non_local_toolbox/neighborhood_base.hh
model/common/non_local_toolbox/neighborhood_base.cc
model/common/non_local_toolbox/neighborhood_base_inline_impl.hh
model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh
model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc
model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.hh
model/common/non_local_toolbox/non_local_manager.hh
model/common/non_local_toolbox/non_local_manager.cc
model/common/non_local_toolbox/non_local_manager_inline_impl.hh
model/common/non_local_toolbox/non_local_manager_callback.hh
model/common/non_local_toolbox/non_local_neighborhood_base.hh
model/common/non_local_toolbox/non_local_neighborhood_base.cc
model/common/non_local_toolbox/non_local_neighborhood.hh
model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh
model/common/non_local_toolbox/non_local_neighborhood_inline_impl.hh
model/common/non_local_toolbox/base_weight_function.hh
model/common/non_local_toolbox/base_weight_function_inline_impl.hh
model/common/model_solver.cc
model/common/model_solver.hh
model/common/solver_callback.hh
model/common/solver_callback.cc
model/common/dof_manager/dof_manager.cc
model/common/dof_manager/dof_manager.hh
model/common/dof_manager/dof_manager_default.cc
model/common/dof_manager/dof_manager_default.hh
model/common/dof_manager/dof_manager_default_inline_impl.hh
model/common/dof_manager/dof_manager_inline_impl.hh
model/common/non_linear_solver/non_linear_solver.cc
model/common/non_linear_solver/non_linear_solver.hh
model/common/non_linear_solver/non_linear_solver_default.hh
model/common/non_linear_solver/non_linear_solver_lumped.cc
model/common/non_linear_solver/non_linear_solver_lumped.hh
model/common/time_step_solvers/time_step_solver.hh
model/common/time_step_solvers/time_step_solver.cc
model/common/time_step_solvers/time_step_solver_default.cc
model/common/time_step_solvers/time_step_solver_default.hh
model/common/time_step_solvers/time_step_solver_default_explicit.hh
model/common/integration_scheme/generalized_trapezoidal.cc
model/common/integration_scheme/generalized_trapezoidal.hh
model/common/integration_scheme/integration_scheme.cc
model/common/integration_scheme/integration_scheme.hh
model/common/integration_scheme/integration_scheme_1st_order.cc
model/common/integration_scheme/integration_scheme_1st_order.hh
model/common/integration_scheme/integration_scheme_2nd_order.cc
model/common/integration_scheme/integration_scheme_2nd_order.hh
model/common/integration_scheme/newmark-beta.cc
model/common/integration_scheme/newmark-beta.hh
model/common/integration_scheme/pseudo_time.cc
model/common/integration_scheme/pseudo_time.hh
model/model.cc
model/model.hh
model/model_inline_impl.hh
model/model_options.hh
solver/solver_vector.hh
solver/solver_vector_default.hh
solver/solver_vector_default_tmpl.hh
solver/solver_vector_distributed.cc
solver/solver_vector_distributed.hh
solver/sparse_matrix.cc
solver/sparse_matrix.hh
solver/sparse_matrix_aij.cc
solver/sparse_matrix_aij.hh
solver/sparse_matrix_aij_inline_impl.hh
solver/sparse_matrix_inline_impl.hh
solver/sparse_solver.cc
solver/sparse_solver.hh
solver/sparse_solver_inline_impl.hh
solver/terms_to_assemble.hh
synchronizer/communication_buffer_inline_impl.hh
synchronizer/communication_descriptor.hh
synchronizer/communication_descriptor_tmpl.hh
synchronizer/communication_request.hh
synchronizer/communication_tag.hh
synchronizer/communications.hh
synchronizer/communications_tmpl.hh
synchronizer/communicator.cc
synchronizer/communicator.hh
synchronizer/communicator_dummy_inline_impl.hh
synchronizer/communicator_event_handler.hh
synchronizer/communicator_inline_impl.hh
synchronizer/data_accessor.cc
synchronizer/data_accessor.hh
synchronizer/dof_synchronizer.cc
synchronizer/dof_synchronizer.hh
synchronizer/dof_synchronizer_inline_impl.hh
synchronizer/element_info_per_processor.cc
synchronizer/element_info_per_processor.hh
synchronizer/element_info_per_processor_tmpl.hh
synchronizer/element_synchronizer.cc
synchronizer/element_synchronizer.hh
synchronizer/facet_synchronizer.cc
synchronizer/facet_synchronizer.hh
synchronizer/facet_synchronizer_inline_impl.hh
synchronizer/grid_synchronizer.cc
synchronizer/grid_synchronizer.hh
synchronizer/grid_synchronizer_tmpl.hh
synchronizer/master_element_info_per_processor.cc
synchronizer/node_info_per_processor.cc
synchronizer/node_info_per_processor.hh
synchronizer/node_synchronizer.cc
synchronizer/node_synchronizer.hh
synchronizer/node_synchronizer_inline_impl.hh
synchronizer/periodic_node_synchronizer.cc
synchronizer/periodic_node_synchronizer.hh
synchronizer/slave_element_info_per_processor.cc
synchronizer/synchronizer.cc
synchronizer/synchronizer.hh
synchronizer/synchronizer_impl.hh
synchronizer/synchronizer_impl_tmpl.hh
synchronizer/synchronizer_registry.cc
synchronizer/synchronizer_registry.hh
synchronizer/synchronizer_tmpl.hh
synchronizer/communication_buffer.hh
)
set(AKANTU_SPIRIT_SOURCES
io/mesh_io/mesh_io_abaqus.cc
io/parser/parser_real.cc
io/parser/parser_random.cc
io/parser/parser_types.cc
io/parser/parser_input_files.cc
PARENT_SCOPE
)
package_declare_elements(core
ELEMENT_TYPES
_point_1
_segment_2
_segment_3
_triangle_3
_triangle_6
_quadrangle_4
_quadrangle_8
_tetrahedron_4
_tetrahedron_10
_pentahedron_6
_pentahedron_15
_hexahedron_8
_hexahedron_20
KIND regular
GEOMETRICAL_TYPES
_gt_point
_gt_segment_2
_gt_segment_3
_gt_triangle_3
_gt_triangle_6
_gt_quadrangle_4
_gt_quadrangle_8
_gt_tetrahedron_4
_gt_tetrahedron_10
_gt_hexahedron_8
_gt_hexahedron_20
_gt_pentahedron_6
_gt_pentahedron_15
INTERPOLATION_TYPES
_itp_lagrange_point_1
_itp_lagrange_segment_2
_itp_lagrange_segment_3
_itp_lagrange_triangle_3
_itp_lagrange_triangle_6
_itp_lagrange_quadrangle_4
_itp_serendip_quadrangle_8
_itp_lagrange_tetrahedron_4
_itp_lagrange_tetrahedron_10
_itp_lagrange_hexahedron_8
_itp_serendip_hexahedron_20
_itp_lagrange_pentahedron_6
_itp_lagrange_pentahedron_15
GEOMETRICAL_SHAPES
_gst_point
_gst_triangle
_gst_square
_gst_prism
GAUSS_INTEGRATION_TYPES
_git_point
_git_segment
_git_triangle
_git_tetrahedron
_git_pentahedron
INTERPOLATION_KIND _itk_lagrangian
FE_ENGINE_LISTS
gradient_on_integration_points
interpolate_on_integration_points
interpolate
compute_normals_on_integration_points
inverse_map
contains
compute_shapes
compute_shapes_derivatives
get_N
compute_dnds
compute_d2nds2
compute_jmat
get_shapes_derivatives
lagrange_base
assemble_fields
)
find_program(READLINK_COMMAND readlink)
find_program(ADDR2LINE_COMMAND addr2line)
find_program(PATCH_COMMAND patch)
mark_as_advanced(READLINK_COMMAND)
mark_as_advanced(ADDR2LINE_COMMAND)
package_declare_extra_files_to_package(core
SOURCES
common/aka_element_classes_info.hh.in
common/aka_config.hh.in
)
if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang" AND (NOT CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.9))
package_set_compile_flags(core CXX "-Wno-undefined-var-template")
endif()
diff --git a/src/model/common/constitutive_laws/constitutive_law.hh b/src/model/common/constitutive_laws/constitutive_law.hh
index eb931a22a..fd97e75a1 100644
--- a/src/model/common/constitutive_laws/constitutive_law.hh
+++ b/src/model/common/constitutive_laws/constitutive_law.hh
@@ -1,238 +1,262 @@
/**
* @file constitutive_law.hh
*
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Mother class for all constitutive laws
*
*
* 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 "data_accessor.hh"
#include "mesh_events.hh"
#include "parsable.hh"
#include "parser.hh"
/* -------------------------------------------------------------------------- */
#include "internal_field.hh"
#include "random_internal_field.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_CONSTITUTIVE_LAW_HH__
#define __AKANTU_CONSTITUTIVE_LAW_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
template <class ConstitutiveLawsHandler_>
class ConstitutiveLaw : public DataAccessor<Element>,
public MeshEventHandler,
public Parsable {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using ConstitutiveLawsHandler = ConstitutiveLawsHandler_;
ConstitutiveLaw(const ConstitutiveLaw & law) = delete;
ConstitutiveLaw & operator=(const ConstitutiveLaw & law) = delete;
/// Initialize constitutive law with defaults
ConstitutiveLaw(ConstitutiveLawsHandler & handler, const ID & id = "",
ElementKind element_kind = _ek_regular);
/// Destructor
~ConstitutiveLaw() override;
protected:
void initialize();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
- template <typename T, int fps = 0>
- void registerInternal(InternalField<T> & /*vect*/) {
- AKANTU_TO_IMPLEMENT();
- }
+ template <typename T>
+ inline void registerInternal(std::shared_ptr<InternalField<T>> & vect);
- template <typename T, int fps = 0>
- void unregisterInternal(InternalField<T> & /*vect*/) {
- AKANTU_TO_IMPLEMENT();
- }
+ inline void unregisterInternal(const ID & id);
/// initialize the constitutive law computed parameter
virtual void initConstitutiveLaw();
/// save the internals in the previous_state if needed
virtual void savePreviousState();
/// restore the internals from previous_state if needed
virtual void restorePreviousState();
/// add an element to the local mesh filter
inline UInt addElement(ElementType type, UInt element, GhostType ghost_type);
inline UInt addElement(const Element & element);
/// add many elements at once
void addElements(const Array<Element> & elements_to_add);
/// remove many element at once
void removeElements(const Array<Element> & elements_to_remove);
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
protected:
/// resize the internals arrrays
virtual void resizeInternals();
/// function called to update the internal parameters when the
/// modifiable parameters are modified
virtual void updateInternalParameters() {}
/// converts global element to local element
inline Element convertToLocalElement(const Element & global_element) const;
/// converts local element to global element
inline Element convertToGlobalElement(const Element & local_element) const;
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
template <typename T>
inline void packElementDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID()) const;
template <typename T>
inline void unpackElementDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID());
/* ------------------------------------------------------------------------ */
/* MeshEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual void
onNodesAdded(const Array<UInt> & /*unused*/,
const NewNodesEvent & /*unused*/) override{};
virtual void
onNodesRemoved(const Array<UInt> & /*unused*/,
const Array<UInt> & /*unused*/,
const RemovedNodesEvent & /*unused*/) override{};
virtual void
onElementsChanged(const Array<Element> & /*unused*/,
const Array<Element> & /*unused*/,
const ElementTypeMapArray<UInt> & /*unused*/,
const ChangedElementsEvent & /*unused*/) override{};
virtual void
onElementsAdded(const Array<Element> & /*unused*/,
const NewElementsEvent & /*unused*/);
virtual void
onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event);
public:
- template <typename T, int fps = 0>
- const InternalField<T> & getInternal(const ID & /*id*/) const {
- AKANTU_TO_IMPLEMENT();
- return NULL;
+ template <typename T>
+ const InternalField<T> & getInternal(const ID & id) const {
+ auto it = internal_vectors.find(getID() + ":" + id);
+ if (it != internal_vectors.end() and
+ aka::is_of_type<InternalField<T>>(*it->second)) {
+ return aka::as_type<InternalField<T>>(*it->second);
+ }
+
+ AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
+ << ") does not contain an internal "
+ << id << " ("
+ << (getID() + ":" + id) << ")");
}
- template <typename T, int fps = 0>
- InternalField<T> & getInternal(const ID & /*id*/) {
- AKANTU_TO_IMPLEMENT();
- return NULL;
+ template <typename T> InternalField<T> & getInternal(const ID & id) {
+ auto it = internal_vectors.find(getID() + ":" + id);
+ if (it != internal_vectors.end() and
+ aka::is_of_type<InternalField<T>>(*it->second)) {
+ return aka::as_type<InternalField<T>>(*it->second);
+ }
+
+ AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
+ << ") does not contain an internal "
+ << id << " ("
+ << (getID() + ":" + id) << ")");
}
- template <typename T, int fps = 0>
- inline bool isInternal(const ID & /*id*/,
- const ElementKind & /*element_kind*/) const {
- return false;
+ template <typename T>
+ inline bool isInternal(const ID & id,
+ const ElementKind & element_kind) const {
+ auto it = internal_vectors.find(this->getID() + ":" + id);
+
+ return (it != internal_vectors.end() and
+ aka::is_of_type<InternalField<T>>(*it->second) and
+ aka::as_type<InternalField<T>>(*it->second).getElementKind() !=
+ element_kind);
}
+ template <typename T>
+ const Array<T> & getArray(const ID & id, ElementType type,
+ GhostType ghost_type = _not_ghost) const;
+ template <typename T>
+ Array<T> & getArray(const ID & id, ElementType type,
+ GhostType ghost_type = _not_ghost);
+
template <typename T> inline void setParam(const ID & param, T value);
inline const Parameter & getParam(const ID & param) const;
template <typename T>
void flattenInternal(const std::string & field_id,
ElementTypeMapArray<T> & internal_flat,
const GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) const;
- /* ------------------------------------------------------------------------ */
- /* Accessors */
- /* ------------------------------------------------------------------------ */
+ /* ------------------------------------------------------------------------
+ */
+ /* Accessors */
+ /* ------------------------------------------------------------------------
+ */
public:
AKANTU_GET_MACRO(Name, name, const std::string &);
AKANTU_GET_MACRO(ID, id, const ID &);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt);
template <typename T>
ElementTypeMap<UInt> getInternalDataPerElem(const ID & id,
ElementKind element_kind) const;
protected:
bool isInit() const { return is_init; }
- /* ------------------------------------------------------------------------ */
- /* Class Members */
- /* ------------------------------------------------------------------------ */
+ /* ------------------------------------------------------------------------
+ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------
+ */
private:
/// boolean to know if the constitutive law has been initialized
bool is_init{false};
- std::map<ID, InternalField<Real> *> internal_vectors_real;
- std::map<ID, InternalField<UInt> *> internal_vectors_uint;
- std::map<ID, InternalField<bool> *> internal_vectors_bool;
+ std::map<ID, std::shared_ptr<InternalFieldBase>> internal_vectors;
protected:
ID id;
/// constitutive law name
std::string name;
/// list of element handled by the constitutive law
ElementTypeMapArray<UInt> element_filter;
// Constitutive law handler for which this is a constitutive law
ConstitutiveLawsHandler & handler;
};
} // namespace akantu
#include "constitutive_law_tmpl.hh"
#endif
diff --git a/src/model/common/constitutive_laws/constitutive_law_tmpl.hh b/src/model/common/constitutive_laws/constitutive_law_tmpl.hh
index baf0fb5ed..d39f106b6 100644
--- a/src/model/common/constitutive_laws/constitutive_law_tmpl.hh
+++ b/src/model/common/constitutive_laws/constitutive_law_tmpl.hh
@@ -1,666 +1,467 @@
/**
* @file constitutive_law_tmpl.hh *
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Implementation of the templated part of the constitutive law class
*
*
* 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 "constitutive_law.hh"
#include "constitutive_laws_handler.hh"
+#include "fe_engine.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_CONSTITUTIVE_LAW_TMPL_HH
#define AKANTU_CONSTITUTIVE_LAW_TMPL_HH
namespace akantu {
template <class ConstitutiveLawsHandler_>
ConstitutiveLaw<ConstitutiveLawsHandler_>::ConstitutiveLaw(
ConstitutiveLawsHandler_ & handler, const ID & id, ElementKind element_kind)
: Parsable(ParserType::_constitutive_law, id), handler(handler), id(id),
element_filter("element_filter", id) {
/// for each connectivity types allocate the element filer array of the
/// constitutive law
element_filter.initialize(handler.getMesh(),
_spatial_dimension = handler.getSpatialDimension(),
_element_kind = element_kind);
this->initialize();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::initialize() {
registerParam("name", name, std::string(), _pat_parsable | _pat_readable);
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::initConstitutiveLaw() {
this->resizeInternals();
this->updateInternalParmaters();
is_init = true;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::savePreviousState() {
- for (auto pair : internal_vectors_real) {
+ for (auto pair : internal_vectors) {
if (pair.second->hasHistory()) {
pair.second->saveCurrentValues();
}
}
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::restorePreviousState() {
- for (auto pair : internal_vectors_real) {
+ for (auto pair : internal_vectors) {
if (pair.second->hasHistory()) {
pair.second->restorePreviousValues();
}
}
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::resizeInternals() {
- for (auto it = internal_vectors_real.begin();
- it != internal_vectors_real.end(); ++it) {
- it->second->resize();
- }
-
- for (auto it = internal_vectors_uint.begin();
- it != internal_vectors_uint.end(); ++it) {
- it->second->resize();
- }
-
- for (auto it = internal_vectors_bool.begin();
- it != internal_vectors_bool.end(); ++it) {
- it->second->resize();
+ for (auto && pair : internal_vectors) {
+ pair.second->resize();
}
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::addElements(
const Array<Element> & elements_to_add) {
AKANTU_DEBUG_IN();
UInt law_id = handler.getConstitutiveLawIndex(name);
for (const auto & element : elements_to_add) {
auto index = this->addElement(element);
handler.constitutive_law_index(element) = law_id;
handler.constitutive_law_local_numbering(element) = index;
}
this->resizeInternals();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::removeElements(
const Array<Element> & elements_to_remove) {
AKANTU_DEBUG_IN();
auto el_begin = elements_to_remove.begin();
auto el_end = elements_to_remove.end();
if (elements_to_remove.empty()) {
return;
}
auto & mesh = handler.getMesh();
ElementTypeMapArray<UInt> constitutive_law_local_new_numbering(
"remove constitutive law filter elem", id);
constitutive_law_local_new_numbering.initialize(
mesh, _element_filter = &element_filter, _element_kind = _ek_not_defined,
_with_nb_element = true);
ElementTypeMapArray<UInt> element_filter_tmp("element_filter_tmp", id);
element_filter_tmp.initialize(mesh, _element_filter = &element_filter,
_element_kind = _ek_not_defined);
ElementTypeMap<UInt> new_ids, element_ids;
for_each_element(
mesh,
[&](auto && el) {
if (not new_ids(el.type, el.ghost_type)) {
element_ids(el.type, el.ghost_type) = 0;
}
auto & element_id = element_ids(el.type, el.ghost_type);
auto l_el = Element{el.type, element_id, el.ghost_type};
if (std::find(el_begin, el_end, el) != el_end) {
constitutive_law_local_new_numbering(l_el) = UInt(-1);
return;
}
element_filter_tmp(el.type, el.ghost_type).push_back(el.element);
if (not new_ids(el.type, el.ghost_type)) {
new_ids(el.type, el.ghost_type) = 0;
}
auto & new_id = new_ids(el.type, el.ghost_type);
constitutive_law_local_new_numbering(l_el) = new_id;
handler.constitutive_law_local_numbering(el) = new_id;
++new_id;
++element_id;
},
_element_filter = &element_filter, _element_kind = _ek_not_defined);
for (auto ghost_type : ghost_types) {
for (const auto & type : element_filter.elementTypes(
_ghost_type = ghost_type, _element_kind = _ek_not_defined)) {
element_filter(type, ghost_type)
.copy(element_filter_tmp(type, ghost_type));
}
}
- for (auto it = internal_vectors_real.begin();
- it != internal_vectors_real.end(); ++it) {
- it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
- }
-
- for (auto it = internal_vectors_uint.begin();
- it != internal_vectors_uint.end(); ++it) {
- it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
- }
-
- for (auto it = internal_vectors_bool.begin();
- it != internal_vectors_bool.end(); ++it) {
- it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
+ for (auto && pair : internal_vectors) {
+ pair.second->removeIntegrationPoints(constitutive_law_local_new_numbering);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::onElementsAdded(
const Array<Element> & /*unused*/, const NewElementsEvent & /*unused*/) {
this->resizeInternals();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::onElementsRemoved(
const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & /*event*/) {
UInt my_num = handler.getInternalIndexFromID(getID());
ElementTypeMapArray<UInt> constitutive_law_local_new_numbering(
"remove constitutive law filter elem", getID());
auto el_begin = element_list.begin();
auto el_end = element_list.end();
for (auto && gt : ghost_types) {
for (auto && type :
new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) {
if (not element_filter.exists(type, gt) ||
element_filter(type, gt).empty()) {
continue;
}
auto & elem_filter = element_filter(type, gt);
auto & law_indexes = this->handler.constitutive_law_index(type, gt);
auto & law_loc_num =
this->handler.constitutive_law_local_numbering(type, gt);
auto nb_element = this->handler.getMesh().getNbElement(type, gt);
// all constitutive laws will resized to the same size...
law_indexes.resize(nb_element);
law_loc_num.resize(nb_element);
if (!constitutive_law_local_new_numbering.exists(type, gt)) {
constitutive_law_local_new_numbering.alloc(elem_filter.size(), 1, type,
gt);
}
auto & law_renumbering = constitutive_law_local_new_numbering(type, gt);
const auto & renumbering = new_numbering(type, gt);
Array<UInt> elem_filter_tmp;
UInt ni = 0;
Element el{type, 0, gt};
for (UInt i = 0; i < elem_filter.size(); ++i) {
el.element = elem_filter(i);
if (std::find(el_begin, el_end, el) == el_end) {
UInt new_el = renumbering(el.element);
AKANTU_DEBUG_ASSERT(new_el != UInt(-1),
"A not removed element as been badly renumbered");
elem_filter_tmp.push_back(new_el);
law_renumbering(i) = ni;
law_indexes(new_el) = my_num;
law_loc_num(new_el) = ni;
++ni;
} else {
law_renumbering(i) = UInt(-1);
}
}
elem_filter.resize(elem_filter_tmp.size());
elem_filter.copy(elem_filter_tmp);
}
}
- for (auto it = internal_vectors_real.begin();
- it != internal_vectors_real.end(); ++it) {
- it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
- }
-
- for (auto it = internal_vectors_uint.begin();
- it != internal_vectors_uint.end(); ++it) {
- it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
- }
-
- for (auto it = internal_vectors_bool.begin();
- it != internal_vectors_bool.end(); ++it) {
- it->second->removeIntegrationPoints(constitutive_law_local_new_numbering);
+ for (auto pair : internal_vectors) {
+ pair.second->removeIntegrationPoints(constitutive_law_local_new_numbering);
}
}
-
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
inline void ConstitutiveLaw<ConstitutiveLawsHandler_>::packElementDataHelper(
const ElementTypeMapArray<T> & data_to_pack, CommunicationBuffer & buffer,
const Array<Element> & elements, const ID & fem_id) const {
DataAccessor::packElementalDataHelper<T>(data_to_pack, buffer, elements, true,
handler.getFEEngine(fem_id));
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
inline void ConstitutiveLaw<ConstitutiveLawsHandler_>::unpackElementDataHelper(
ElementTypeMapArray<T> & data_to_unpack, CommunicationBuffer & buffer,
const Array<Element> & elements, const ID & fem_id) {
DataAccessor::unpackElementalDataHelper<T>(data_to_unpack, buffer, elements,
true, handler.getFEEngine(fem_id));
-}
-
+}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
-inline Element
-ConstitutiveLaw<ConstitutiveLawsHandler_>::convertToLocalElement(const Element & global_element) const {
+inline Element ConstitutiveLaw<ConstitutiveLawsHandler_>::convertToLocalElement(
+ const Element & global_element) const {
UInt ge = global_element.element;
#ifndef AKANTU_NDEBUG
UInt model_law_index = handler.getConstitutiveLawByElement(
global_element.type, global_element.ghost_type)(ge);
UInt law_index = handler.getConstitutiveLawIndex(this->name);
AKANTU_DEBUG_ASSERT(model_law_index == law_index,
"Conversion of a global element in a local element for "
"the wrong constitutive law "
<< this->name << std::endl);
#endif
UInt le = handler.getConstitutiveLawLocalNumbering(
global_element.type, global_element.ghost_type)(ge);
Element tmp_quad{global_element.type, le, global_element.ghost_type};
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline Element
-ConstitutiveLaw<ConstitutiveLawsHandler_>::convertToGlobalElement(const Element & local_element) const {
+ConstitutiveLaw<ConstitutiveLawsHandler_>::convertToGlobalElement(
+ const Element & local_element) const {
UInt le = local_element.element;
UInt ge =
this->element_filter(local_element.type, local_element.ghost_type)(le);
Element tmp_quad{local_element.type, ge, local_element.ghost_type};
return tmp_quad;
}
-
-
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline UInt ConstitutiveLaw<ConstitutiveLawsHandler_>::addElement(
ElementType type, UInt element, GhostType ghost_type) {
Array<UInt> & el_filter = this->element_filter(type, ghost_type);
el_filter.push_back(element);
return el_filter.size() - 1;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline UInt
ConstitutiveLaw<ConstitutiveLawsHandler_>::addElement(const Element & element) {
return this->addElement(element.type, element.element, element.ghost_type);
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
inline const Parameter &
ConstitutiveLaw<ConstitutiveLawsHandler_>::getParam(const ID & param) const {
try {
return get(param);
} catch (...) {
AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law"
<< getID());
}
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
template <typename T>
inline void
ConstitutiveLaw<ConstitutiveLawsHandler_>::setParam(const ID & param, T value) {
try {
set<T>(param, value);
} catch (...) {
AKANTU_EXCEPTION("No parameter " << param << " in the constitutive law "
<< getID());
}
updateInternalParameters();
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline void
-ConstitutiveLaw<ConstitutiveLawsHandler_>::registerInternal<Real, fps>(
- InternalField<Real> & vect) {
- internal_vectors_real[vect.getID()] = &vect;
-}
-
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline void
-ConstitutiveLaw<ConstitutiveLawsHandler_>::registerInternal<UInt, fps>(
- InternalField<UInt> & vect) {
- internal_vectors_uint[vect.getID()] = &vect;
-}
-
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline void
-ConstitutiveLaw<ConstitutiveLawsHandler_>::registerInternal<bool, fps>(
- InternalField<bool> & vect) {
- internal_vectors_bool[vect.getID()] = &vect;
+template <typename T>
+inline void ConstitutiveLaw<ConstitutiveLawsHandler_>::registerInternal(
+ std::shared_ptr<InternalField<T>> & vect) {
+ internal_vectors[vect.getID()] = vect;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
-template <int fps>
inline void
-ConstitutiveLaw<ConstitutiveLawsHandler_>::unregisterInternal<Real, fps>(
- InternalField<Real> & vect) {
- internal_vectors_real.erase(vect.getID());
-}
-
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline void
-ConstitutiveLaw<ConstitutiveLawsHandler_>::unregisterInternal<UInt, fps>(
- InternalField<UInt> & vect) {
- internal_vectors_uint.erase(vect.getID());
-}
-
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline void
-ConstitutiveLaw<ConstitutiveLawsHandler_>::unregisterInternal<bool, fps>(
- InternalField<bool> & vect) {
- internal_vectors_bool.erase(vect.getID());
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-const InternalField<Real> &
-ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal<Real, fps>(
- const ID & int_id) const {
- auto it = internal_vectors_real.find(getID() + ":" + int_id);
- if (it == internal_vectors_real.end()) {
- AKANTU_SILENT_EXCEPTION("The constitutive law " << name << "(" << getID()
- << ") does not contain an internal "
- << int_id << " ("
- << (getID() + ":" + int_id) << ")");
- }
- return *it->second;
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <>
-InternalField<Real, fps> &
-ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal<Real, fps>(
- const ID & int_id) {
- auto it = internal_vectors_real.find(getID() + ":" + int_id);
- if (it == internal_vectors_real.end()) {
- AKANTU_SILENT_EXCEPTION("The constitutive law " << name << "(" << getID()
- << ") does not contain an internal "
- << int_id << " ("
- << (getID() + ":" + int_id) << ")");
- }
- return *it->second;
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-const InternalField<UInt> &
-ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal<UInt, fps>(
- const ID & int_id) const {
- auto it = internal_vectors_uint.find(getID() + ":" + int_id);
- if (it == internal_vectors_uint.end()) {
- AKANTU_SILENT_EXCEPTION("The constitutive law " << name << "(" << getID()
- << ") does not contain an internal "
- << int_id << " ("
- << (getID() + ":" + int_id) << ")");
- }
- return *it->second;
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-InternalField<UInt> &
-ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal<UInt, fps>(
- const ID & int_id) {
- auto it = internal_vectors_uint.find(getID() + ":" + int_id);
- if (it == internal_vectors_uint.end()) {
- AKANTU_SILENT_EXCEPTION("The constitutive law " << name << "(" << getID()
- << ") does not contain an internal "
- << int_id << " ("
- << (getID() + ":" + int_id) << ")");
- }
- return *it->second;
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-const InternalField<bool> &
-ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal<bool, fps>(
- const ID & int_id) const {
- auto it = internal_vectors_bool.find(getID() + ":" + int_id);
- if (it == internal_vectors_bool.end()) {
- AKANTU_SILENT_EXCEPTION("The constitutive law " << name << "(" << getID()
- << ") does not contain an internal "
- << int_id << " ("
- << (getID() + ":" + int_id) << ")");
- }
- return *it->second;
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-InternalField<bool> &
-ConstituitveLaw<ConstitutiveLawsHandler_>::getInternal<bool, fps>(
- const ID & int_id) {
- auto it = internal_vectors_bool.find(getID() + ":" + int_id);
- if (it == internal_vectors_bool.end()) {
- AKANTU_SILENT_EXCEPTION("The constitutive law " << name << "(" << getID()
- << ") does not contain an internal "
- << int_id << " ("
- << (getID() + ":" + int_id) << ")");
- }
- return *it->second;
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline bool ConstitutiveLaw<ConstitutiveLawsHandler_>::isInternal<Real, fps>(
- const ID & id, ElementKind element_kind) const {
- auto internal_array = internal_vectors_real.find(this->getID() + ":" + id);
-
- return not(internal_array == internal_vectors_real.end() ||
- internal_array->second->getElementKind() != element_kind);
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline bool ConstitutiveLaw<ConstitutiveLawsHandler_>::isInternal<UInt, fps>(
- const ID & id, ElementKind element_kind) const {
- auto internal_array = internal_vectors_uint.find(this->getID() + ":" + id);
-
- return not(internal_array == internal_vectors_uint.end() ||
- internal_array->second->getElementKind() != element_kind);
-}
-
-/* -------------------------------------------------------------------------- */
-template <class ConstitutiveLawsHandler_>
-template <int fps>
-inline bool ConstitutiveLaw<ConstitutiveLawsHandler_>::isInternal<bool, fps>(
- const ID & id, ElementKind element_kind) const {
- auto internal_array = internal_vectors_bool.find(this->getID() + ":" + id);
-
- return not(internal_array == internal_vectors_bool.end() ||
- internal_array->second->getElementKind() != element_kind);
+ConstitutiveLaw<ConstitutiveLawsHandler_>::unregisterInternal(const ID & id) {
+ internal_vectors.erase(id);
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
+template <typename T>
inline ElementTypeMap<UInt>
ConstitutiveLaw<ConstitutiveLawsHandler_>::getInternalDataPerElem(
const ID & field_id, ElementKind element_kind) const {
- if (!this->template isInternal<T>(field_id, element_kind)) {
+ if (not this->template isInternal<T>(field_id, element_kind)) {
AKANTU_EXCEPTION("Cannot find internal field "
<< id << " in the constitutive law " << this->name);
}
- const InternalField<T> & internal_field =
- this->template getInternal<T>(field_id);
+ const auto & internal_field = this->template getInternal<T>(field_id);
const FEEngine & fe_engine = internal_field.getFEEngine();
UInt nb_data_per_quad = internal_field.getNbComponent();
ElementTypeMap<UInt> res;
for (auto ghost_type : ghost_types) {
for (auto & type : internal_field.elementTypes(ghost_type)) {
UInt nb_quadrature_points =
fe_engine.getNbIntegrationPoints(type, ghost_type);
res(type, ghost_type) = nb_data_per_quad * nb_quadrature_points;
}
}
return res;
}
/* -------------------------------------------------------------------------- */
template <class ConstitutiveLawsHandler_>
+template <typename T>
void ConstitutiveLaw<ConstitutiveLawsHandler_>::flattenInternal(
const std::string & field_id, ElementTypeMapArray<T> & internal_flat,
const GhostType ghost_type, ElementKind element_kind) const {
if (!this->template isInternal<T>(field_id, element_kind)) {
AKANTU_EXCEPTION("Cannot find internal field "
<< id << " in the constitutive law " << this->name);
}
const InternalField<T> & internal_field =
this->template getInternal<T>(field_id);
const FEEngine & fe_engine = internal_field.getFEEngine();
const Mesh & mesh = fe_engine.getMesh();
for (auto && type : internal_field.filterTypes(ghost_type)) {
const Array<Real> & src_vect = internal_field(type, ghost_type);
const Array<UInt> & filter = internal_field.getFilter(type, ghost_type);
// total number of elements in the corresponding mesh
UInt nb_element_dst = mesh.getNbElement(type, ghost_type);
// number of element in the internal field
UInt nb_element_src = filter.size();
// number of quadrature points per elem
UInt nb_quad_per_elem = fe_engine.getNbIntegrationPoints(type);
// number of data per quadrature point
UInt nb_data_per_quad = internal_field.getNbComponent();
if (!internal_flat.exists(type, ghost_type)) {
internal_flat.alloc(nb_element_dst * nb_quad_per_elem, nb_data_per_quad,
type, ghost_type);
}
if (nb_element_src == 0) {
continue;
}
// number of data per element
UInt nb_data = nb_quad_per_elem * nb_data_per_quad;
Array<Real> & dst_vect = internal_flat(type, ghost_type);
dst_vect.resize(nb_element_dst * nb_quad_per_elem);
auto it_dst = make_view(dst_vect, nb_data).begin();
for (auto && data : zip(filter, make_view(src_vect, nb_data))) {
it_dst[std::get<0>(data)] = std::get<1>(data);
}
}
}
-
} // namespace akantu
#endif // AKANTU_CONSTITUTIVE_LAW_TMPL_HH
diff --git a/src/model/common/constitutive_laws/internal_field.hh b/src/model/common/constitutive_laws/internal_field.hh
index 8632a87a0..a8f5c1a89 100644
--- a/src/model/common/constitutive_laws/internal_field.hh
+++ b/src/model/common/constitutive_laws/internal_field.hh
@@ -1,252 +1,277 @@
/**
* @file internal_field.hh
*
* @author Lucas Frerot <lucas.frerot@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Mar 26 2021
*
* @brief Constitutive law internal properties
*
*
* @section LICENSE
*
* Copyright (©) 2010-2021 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 "element_type_map.hh"
/* -------------------------------------------------------------------------- */
+#include <memory>
+/* -------------------------------------------------------------------------- */
#ifndef AKANTU_INTERNAL_FIELD_HH_
#define AKANTU_INTERNAL_FIELD_HH_
namespace akantu {
class ConstitutiveLaw;
class FEEngine;
+class InternalFieldBase
+ : public std::enable_shared_from_this<InternalFieldBase> {
+public:
+ /// activate the history of this field
+ virtual void initializeHistory() = 0;
+
+ /// resize the arrays and set the new element to 0
+ virtual void resize() = 0;
+
+ /// save the current values in the history
+ virtual void saveCurrentValues() = 0;
+
+ /// restore the previous values from the history
+ virtual void restorePreviousValues() = 0;
+
+ /// remove the quadrature points corresponding to suppressed elements
+ virtual void
+ removeIntegrationPoints(const ElementTypeMapArray<UInt> & new_numbering) = 0;
+
+ virtual bool hasHistory() const = 0;
+};
+
/**
* class for the internal fields of constitutive law
* to store values for each quadrature
*/
template <class ConstitutiveLaw_, typename T>
-class InternalFieldTmpl : public ElementTypeMapArray<T> {
+class InternalFieldTmpl : public InternalFieldBase,
+ public ElementTypeMapArray<T> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using ConstitutiveLaw = ConstitutiveLaw_;
InternalFieldTmpl(const ID & id, ConstitutiveLaw & constitutive_law);
~InternalFieldTmpl() override;
/// This constructor is only here to let cohesive elements compile
InternalFieldTmpl(const ID & id, ConstitutiveLaw & constitutive_law,
FEEngine & fem,
const ElementTypeMapArray<UInt> & element_filter);
/// More general constructor
InternalFieldTmpl(const ID & id, ConstitutiveLaw & constitutive_law, UInt dim,
FEEngine & fem,
const ElementTypeMapArray<UInt> & element_filter);
InternalFieldTmpl(const ID & id,
const InternalFieldTmpl<ConstitutiveLaw_, T> & other);
InternalFieldTmpl operator=(const InternalFieldTmpl &) = delete;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// function to reset the FEEngine for the internal fieldx
virtual void setFEEngine(FEEngine & fe_engine);
/// function to reset the element kind for the internal
virtual void setElementKind(ElementKind element_kind);
/// initialize the field to a given number of component
virtual void initialize(UInt nb_component);
/// activate the history of this field
- virtual void initializeHistory();
+ void initializeHistory() override;
/// resize the arrays and set the new element to 0
- virtual void resize();
+ void resize() override;
/// set the field to a given value v
virtual void setDefaultValue(const T & v);
/// reset all the fields to the default value
virtual void reset();
/// save the current values in the history
- virtual void saveCurrentValues();
+ void saveCurrentValues() override;
/// restore the previous values from the history
- virtual void restorePreviousValues();
+ void restorePreviousValues() override;
/// remove the quadrature points corresponding to suppressed elements
- virtual void
- removeIntegrationPoints(const ElementTypeMapArray<UInt> & new_numbering);
+ void removeIntegrationPoints(
+ const ElementTypeMapArray<UInt> & new_numbering) override;
/// print the content
void printself(std::ostream & stream, int /*indent*/ = 0) const override;
/// get the default value
inline operator T() const;
virtual FEEngine & getFEEngine() { return *fem; }
virtual const FEEngine & getFEEngine() const { return *fem; }
protected:
/// initialize the arrays in the ElementTypeMapArray<T>
void internalInitialize(UInt nb_component);
/// set the values for new internals
virtual void setArrayValues(T * begin, T * end);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get filter types for range loop
decltype(auto) elementTypes(GhostType ghost_type = _not_ghost) const {
return ElementTypeMapArray<T>::elementTypes(
_spatial_dimension = this->spatial_dimension,
_element_kind = this->element_kind, _ghost_type = ghost_type);
}
/// get filter types for range loop
decltype(auto) filterTypes(GhostType ghost_type = _not_ghost) const {
return this->element_filter.elementTypes(
_spatial_dimension = this->spatial_dimension,
_element_kind = this->element_kind, _ghost_type = ghost_type);
}
/// get the array for a given type of the element_filter
const Array<UInt> & getFilter(ElementType type,
GhostType ghost_type = _not_ghost) const {
return this->element_filter(type, ghost_type);
}
/// get the Array corresponding to the type en ghost_type specified
virtual Array<T> & operator()(ElementType type,
GhostType ghost_type = _not_ghost) {
return ElementTypeMapArray<T>::operator()(type, ghost_type);
}
virtual const Array<T> & operator()(ElementType type,
GhostType ghost_type = _not_ghost) const {
return ElementTypeMapArray<T>::operator()(type, ghost_type);
}
virtual Array<T> & previous(ElementType type,
GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_ASSERT(previous_values != nullptr,
"The history of the internal "
<< this->getID() << " has not been activated");
return this->previous_values->operator()(type, ghost_type);
}
virtual const Array<T> & previous(ElementType type,
GhostType ghost_type = _not_ghost) const {
AKANTU_DEBUG_ASSERT(previous_values != nullptr,
"The history of the internal "
<< this->getID() << " has not been activated");
return this->previous_values->operator()(type, ghost_type);
}
virtual InternalFieldTmpl<ConstitutiveLaw_, T> & previous() {
AKANTU_DEBUG_ASSERT(previous_values != nullptr,
"The history of the internal "
<< this->getID() << " has not been activated");
return *(this->previous_values);
}
virtual const InternalFieldTmpl<ConstitutiveLaw_, T> & previous() const {
AKANTU_DEBUG_ASSERT(previous_values != nullptr,
"The history of the internal "
<< this->getID() << " has not been activated");
return *(this->previous_values);
}
/// check if the history is used or not
- bool hasHistory() const { return (previous_values != nullptr); }
+ bool hasHistory() const override { return (previous_values != nullptr); }
/// get the kind treated by the internal
ElementKind getElementKind() const { return element_kind; }
/// return the number of components
UInt getNbComponent() const { return nb_component; }
/// return the spatial dimension corresponding to the internal element type
/// loop filter
UInt getSpatialDimension() const { return this->spatial_dimension; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// the constitutive_law for which this is an internal parameter
ConstitutiveLaw & constitutive_law;
/// the fem containing the mesh and the element informations
FEEngine * fem{nullptr};
/// Element filter if needed
const ElementTypeMapArray<UInt> & element_filter;
/// default value
T default_value{};
/// spatial dimension of the element to consider
UInt spatial_dimension{0};
/// ElementKind of the element to consider
ElementKind element_kind{_ek_regular};
/// Number of component of the internal field
UInt nb_component{0};
/// Is the field initialized
bool is_init{false};
/// previous values
std::unique_ptr<InternalFieldTmpl<ConstitutiveLaw_, T>> previous_values;
};
/// standard output stream operator
template <class ConstitutiveLaw_, typename T>
inline std::ostream &
operator<<(std::ostream & stream,
const InternalFieldTmpl<ConstitutiveLaw_, T> & _this) {
_this.printself(stream);
return stream;
}
template <typename T>
using InternalField = InternalFieldTmpl<ConstitutiveLaw_, T>;
} // namespace akantu
#endif /* AKANTU_INTERNAL_FIELD_HH_ */
diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index 49256cfdf..c50fabcdc 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,595 +1,594 @@
/**
* @file material.hh
*
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Lucas Frerot <lucas.frerot@epfl.ch>
* @author Enrico Milanese <enrico.milanese@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Apr 09 2021
*
* @brief Mother class for all materials
*
*
* @section LICENSE
*
* Copyright (©) 2010-2021 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_factory.hh"
#include "aka_voigthelper.hh"
#include "constitutive_law.hh"
#include "data_accessor.hh"
#include "integration_point.hh"
#include "parsable.hh"
#include "parser.hh"
/* -------------------------------------------------------------------------- */
#include "internal_field.hh"
#include "random_internal_field.hh"
/* -------------------------------------------------------------------------- */
#include "mesh_events.hh"
#include "solid_mechanics_model_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_HH_
#define AKANTU_MATERIAL_HH_
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class SolidMechanicsModel;
class Material;
} // namespace akantu
namespace akantu {
using MaterialFactory =
Factory<Material, ID, UInt, const ID &, SolidMechanicsModel &, const ID &>;
/**
* Interface of all materials
* Prerequisites for a new material
* - inherit from this class
* - implement the following methods:
* \code
* virtual Real getStableTimeStep(Real h, const Element & element =
* ElementNull);
*
* virtual void computeStress(ElementType el_type,
* GhostType ghost_type = _not_ghost);
*
* virtual void computeTangentStiffness(ElementType el_type,
* Array<Real> & tangent_matrix,
* GhostType ghost_type = _not_ghost);
* \endcode
*
*/
-class Material
- : public ConstitutiveLaw<SolidMechanicsModel> public DataAccessor<Element>,
- public Parsable,
- protected SolidMechanicsModelEventHandler {
+class Material : public ConstitutiveLaw<SolidMechanicsModel>,
+ public MeshEventHandler,
+ protected SolidMechanicsModelEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Material(const Material & mat) = delete;
Material & operator=(const Material & mat) = delete;
/// Initialize material with defaults
Material(SolidMechanicsModel & model, const ID & id = "");
/// Initialize material with custom mesh & fe_engine
Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh,
FEEngine & fe_engine, const ID & id = "");
/// Destructor
~Material() override;
protected:
void initialize();
/* ------------------------------------------------------------------------ */
/* Function that materials can/should reimplement */
/* ------------------------------------------------------------------------ */
protected:
/// constitutive law
virtual void computeStress(ElementType /* el_type */,
GhostType /* ghost_type */ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the tangent stiffness matrix
virtual void computeTangentModuli(ElementType /*el_type*/,
Array<Real> & /*tangent_matrix*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the potential energy
virtual void computePotentialEnergy(ElementType el_type);
/// compute the potential energy for an element
virtual void
computePotentialEnergyByElement(ElementType /*type*/, UInt /*index*/,
Vector<Real> & /*epot_on_quad_points*/) {
AKANTU_TO_IMPLEMENT();
}
virtual void updateEnergies(ElementType /*el_type*/) {}
virtual void updateEnergiesAfterDamage(ElementType /*el_type*/) {}
/// set the material to steady state (to be implemented for materials that
/// need it)
virtual void setToSteadyState(ElementType /*el_type*/,
GhostType /*ghost_type*/ = _not_ghost) {}
/// function called to update the internal parameters when the
/// modifiable parameters are modified
virtual void updateInternalParameters();
public:
/// extrapolate internal values
virtual void extrapolateInternal(const ID & id, const Element & element,
const Matrix<Real> & points,
Matrix<Real> & extrapolated);
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed(const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed(const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
/// get a material celerity to compute the stable time step (default: is the
/// push wave speed)
virtual Real getCelerity(const Element & element) const {
return getPushWaveSpeed(element);
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material computed parameter
virtual void initMaterial();
/// compute the residual for this material
// virtual void updateResidual(GhostType ghost_type = _not_ghost);
/// assemble the residual for this material
virtual void assembleInternalForces(GhostType ghost_type);
/// compute the stresses for this material
virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
// virtual void
// computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost);
virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost);
/// set material to steady state
void setToSteadyState(GhostType ghost_type = _not_ghost);
/// compute the stiffness matrix
virtual void assembleStiffnessMatrix(GhostType ghost_type);
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points
*/
void interpolateStress(ElementTypeMapArray<Real> & result,
GhostType ghost_type = _not_ghost);
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points and store the
* results per facet
*/
void interpolateStressOnFacets(ElementTypeMapArray<Real> & result,
ElementTypeMapArray<Real> & by_elem_result,
GhostType ghost_type = _not_ghost);
/**
* function to initialize the elemental field interpolation
* function by inverting the quadrature points' coordinates
*/
void initElementalFieldInterpolation(
const ElementTypeMapArray<Real> & interpolation_points_coordinates);
/* ------------------------------------------------------------------------ */
/* Common part */
/* ------------------------------------------------------------------------ */
protected:
/* ------------------------------------------------------------------------ */
static inline UInt getTangentStiffnessVoigtSize(UInt dim);
/// compute the potential energy by element
void computePotentialEnergyByElements();
/* ------------------------------------------------------------------------ */
/* Finite deformation functions */
/* This functions area implementing what is described in the paper of Bathe */
/* et al, in IJNME, Finite Element Formulations For Large Deformation */
/* Dynamic Analysis, Vol 9, 353-386, 1975 */
/* ------------------------------------------------------------------------ */
protected:
/// assemble the residual
template <UInt dim> void assembleInternalForces(GhostType ghost_type);
template <UInt dim>
void computeAllStressesFromTangentModuli(ElementType type,
GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrix(ElementType type, GhostType ghost_type);
/// assembling in finite deformation
template <UInt dim>
void assembleStiffnessMatrixNL(ElementType type, GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrixL2(ElementType type, GhostType ghost_type);
/* ------------------------------------------------------------------------ */
/* Conversion functions */
/* ------------------------------------------------------------------------ */
public:
/// Size of the Stress matrix for the case of finite deformation see: Bathe et
/// al, IJNME, Vol 9, 353-386, 1975
static inline UInt getCauchyStressMatrixSize(UInt dim);
/// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386,
/// 1975
template <UInt dim>
static inline void setCauchyStressMatrix(const Matrix<Real> & S_t,
Matrix<Real> & sigma);
/// write the stress tensor in the Voigt notation.
template <UInt dim>
static inline decltype(auto) stressToVoigt(const Matrix<Real> & stress) {
return VoigtHelper<dim>::matrixToVoigt(stress);
}
/// write the strain tensor in the Voigt notation.
template <UInt dim>
static inline decltype(auto) strainToVoigt(const Matrix<Real> & strain) {
return VoigtHelper<dim>::matrixToVoigtWithFactors(strain);
}
/// write a voigt vector to stress
template <UInt dim>
static inline void voigtToStress(const Vector<Real> & voigt,
Matrix<Real> & stress) {
VoigtHelper<dim>::voigtToMatrix(voigt, stress);
}
/// Computation of Cauchy stress tensor in the case of finite deformation from
/// the 2nd Piola-Kirchhoff for a given element type
template <UInt dim>
void StoCauchy(ElementType el_type, GhostType ghost_type = _not_ghost);
/// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation
/// gradient
template <UInt dim>
inline void StoCauchy(const Matrix<Real> & F, const Matrix<Real> & S,
Matrix<Real> & sigma, const Real & C33 = 1.0) const;
template <UInt dim>
static inline void gradUToF(const Matrix<Real> & grad_u, Matrix<Real> & F);
template <UInt dim>
static inline decltype(auto) gradUToF(const Matrix<Real> & grad_u);
static inline void rightCauchy(const Matrix<Real> & F, Matrix<Real> & C);
static inline void leftCauchy(const Matrix<Real> & F, Matrix<Real> & B);
template <UInt dim>
static inline void gradUToEpsilon(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
template <UInt dim>
static inline decltype(auto) gradUToEpsilon(const Matrix<Real> & grad_u);
template <UInt dim>
static inline void gradUToE(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
template <UInt dim>
static inline decltype(auto) gradUToE(const Matrix<Real> & grad_u);
static inline Real stressToVonMises(const Matrix<Real> & stress);
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) override;
/* ------------------------------------------------------------------------ */
/* SolidMechanicsModelEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual void beforeSolveStep();
virtual void afterSolveStep(bool converged = true);
void onDamageIteration() override;
void onDamageUpdate() override;
void onDump() override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &)
AKANTU_GET_MACRO(Rho, rho, Real);
AKANTU_SET_MACRO(Rho, rho, Real);
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// return the potential energy for the subset of elements contained by the
/// material
Real getPotentialEnergy();
/// return the potential energy for the provided element
Real getPotentialEnergy(ElementType & type, UInt index);
/// return the energy (identified by id) for the subset of elements contained
/// by the material
virtual Real getEnergy(const std::string & type);
/// return the energy (identified by id) for the provided element
virtual Real getEnergy(const std::string & energy_id, ElementType type,
UInt index);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy,
Real);
AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(ElementFilter, element_filter,
const ElementTypeMapArray<UInt> &);
AKANTU_GET_MACRO(FEEngine, fem, FEEngine &);
bool isNonLocal() const { return is_non_local; }
bool isFiniteDeformation() const { return finite_deformation; }
bool isInelasticDeformation() const { return inelastic_deformation; }
template <typename T>
void inflateInternal(const std::string & field_id,
const ElementTypeMapArray<T> & field,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined);
/// apply a constant eigengrad_u everywhere in the material
virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
GhostType /*ghost_type*/ = _not_ghost);
bool hasMatrixChanged(const ID & id) {
if (id == "K") {
return hasStiffnessMatrixChanged() or finite_deformation;
}
return true;
}
MatrixType getMatrixType(const ID & id) {
if (id == "K") {
return getTangentType();
}
if (id == "M") {
return _symmetric;
}
return _mt_not_defined;
}
/// specify if the matrix need to be recomputed for this material
virtual bool hasStiffnessMatrixChanged() { return true; }
/// specify the type of matrix, if not overloaded the material is not valid
/// for static or implicit computations
virtual MatrixType getTangentType() { return _mt_not_defined; }
/// static method to reteive the material factory
static MaterialFactory & getFactory();
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// Link to the fem object in the model
FEEngine & fem;
/// Finite deformation
bool finite_deformation{false};
/// Finite deformation
bool inelastic_deformation{false};
/// The model to witch the material belong
SolidMechanicsModel & model;
/// density : rho
Real rho{0.};
/// spatial dimension
UInt spatial_dimension;
/// stresses arrays ordered by element types
InternalField<Real> stress;
/// eigengrad_u arrays ordered by element types
InternalField<Real> eigengradu;
/// grad_u arrays ordered by element types
InternalField<Real> gradu;
/// Green Lagrange strain (Finite deformation)
InternalField<Real> green_strain;
/// Second Piola-Kirchhoff stress tensor arrays ordered by element types
/// (Finite deformation)
InternalField<Real> piola_kirchhoff_2;
/// potential energy by element
InternalField<Real> potential_energy;
/// tell if using in non local mode or not
bool is_non_local{false};
/// tell if the material need the previous stress state
bool use_previous_stress{false};
/// tell if the material need the previous strain state
bool use_previous_gradu{false};
/// elemental field interpolation coordinates
InternalField<Real> interpolation_inverse_coordinates;
/// elemental field interpolation points
InternalField<Real> interpolation_points_matrices;
private:
/// eigen_grad_u for the parser
Matrix<Real> eigen_grad_u;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const Material & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "material_inline_impl.hh"
#include "internal_field_tmpl.hh"
#include "random_internal_field_tmpl.hh"
/* -------------------------------------------------------------------------- */
/* Auto loop */
/* -------------------------------------------------------------------------- */
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeStress. This macro in addition to write the loop
/// provides two tensors (matrices) sigma and grad_u
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \
auto && grad_u_view = \
make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto stress_view = \
make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
if (this->isFiniteDeformation()) { \
stress_view = make_view(this->piola_kirchhoff_2(el_type, ghost_type), \
this->spatial_dimension, this->spatial_dimension); \
} \
\
for (auto && data : zip(grad_u_view, stress_view)) { \
[[gnu::unused]] Matrix<Real> & grad_u = std::get<0>(data); \
[[gnu::unused]] Matrix<Real> & sigma = std::get<1>(data)
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END }
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeTangentModuli. This macro in addition to write the
/// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix
/// where the elemental tangent moduli should be stored in Voigt Notation
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \
auto && grad_u_view = \
make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto && stress_view = \
make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto tangent_size = \
this->getTangentStiffnessVoigtSize(this->spatial_dimension); \
\
auto && tangent_view = make_view(tangent_mat, tangent_size, tangent_size); \
\
for (auto && data : zip(grad_u_view, stress_view, tangent_view)) { \
[[gnu::unused]] Matrix<Real> & grad_u = std::get<0>(data); \
[[gnu::unused]] Matrix<Real> & sigma = std::get<1>(data); \
Matrix<Real> & tangent = std::get<2>(data);
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END }
/* -------------------------------------------------------------------------- */
#define INSTANTIATE_MATERIAL_ONLY(mat_name) \
template class mat_name<1>; /* NOLINT */ \
template class mat_name<2>; /* NOLINT */ \
template class mat_name<3> /* NOLINT */
#define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name) \
[](UInt dim, const ID &, SolidMechanicsModel & model, \
const ID & id) /* NOLINT */ \
-> std::unique_ptr< \
Material> { /* NOLINT */ \
switch (dim) { \
case 1: \
return std::make_unique<mat_name<1>>(/* NOLINT */ \
model, id); \
case 2: \
return std::make_unique<mat_name<2>>(/* NOLINT */ \
model, id); \
case 3: \
return std::make_unique<mat_name<3>>(/* NOLINT */ \
model, id); \
default: \
AKANTU_EXCEPTION( \
"The dimension " \
<< dim \
<< "is not a valid dimension for the material " \
<< #id); \
} \
}
#define INSTANTIATE_MATERIAL(id, mat_name) \
INSTANTIATE_MATERIAL_ONLY(mat_name); \
static bool material_is_alocated_##id = \
MaterialFactory::getInstance().registerAllocator( \
#id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name))
#endif /* AKANTU_MATERIAL_HH_ */