diff --git a/packages/core.cmake b/packages/core.cmake index 10660f10e..c8d23637e 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,433 +1,440 @@ #=============================================================================== # @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/model_io.cc - #io/model_io.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_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/packages/solid_mechanics.cmake b/packages/solid_mechanics.cmake index c7b8776ac..4b86b201f 100644 --- a/packages/solid_mechanics.cmake +++ b/packages/solid_mechanics.cmake @@ -1,118 +1,110 @@ #=============================================================================== # @file solid_mechanics.cmake # # @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> # @author Nicolas Richart <nicolas.richart@epfl.ch> # # @date creation: Mon Dec 04 2017 # @date last modification: Fri Mar 26 2021 # # @brief package description for core # # # @section LICENSE # # Copyright (©) 2016-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(solid_mechanics DEFAULT ON DESCRIPTION "Solid mechanics model" DEPENDS core lapack ) package_declare_sources(solid_mechanics model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.hh model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh - model/solid_mechanics/materials/internal_field.hh - model/solid_mechanics/materials/internal_field_tmpl.hh - model/solid_mechanics/materials/random_internal_field.hh - model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh - model/solid_mechanics/solid_mechanics_model_inline_impl.hh model/solid_mechanics/solid_mechanics_model_io.cc model/solid_mechanics/solid_mechanics_model_mass.cc - model/solid_mechanics/solid_mechanics_model_material.cc - model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.hh model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.cc model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.hh model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_phasefield.cc model/solid_mechanics/materials/material_damage/material_phasefield.hh model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.hh model/solid_mechanics/materials/material_damage/material_von_mises_mazars.cc model/solid_mechanics/materials/material_damage/material_von_mises_mazars.hh model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh model/solid_mechanics/materials/material_non_local.hh - model/solid_mechanics/materials/material_non_local_tmpl.hh model/solid_mechanics/materials/material_non_local_includes.hh ) package_declare_material_infos(solid_mechanics LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_extra_files_to_package(solid_mechanics SOURCES model/solid_mechanics/material_list.hh.in ) diff --git a/src/io/parser/parser.hh b/src/io/parser/parser.hh index fa1966f8f..3302f921f 100644 --- a/src/io/parser/parser.hh +++ b/src/io/parser/parser.hh @@ -1,543 +1,544 @@ /** * @file parser.hh * * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Wed Nov 13 2013 * @date last modification: Fri Apr 02 2021 * * @brief File parser interface * * * @section LICENSE * * Copyright (©) 2014-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 "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ #include <map> /* -------------------------------------------------------------------------- */ #ifndef AKANTU_PARSER_HH_ #define AKANTU_PARSER_HH_ namespace akantu { #if !defined(DOXYGEN) // clang-format off #define AKANTU_SECTION_TYPES \ (cohesive_inserter) \ (contact) \ (embedded_interface) \ (friction) \ (global) \ (heat) \ (integration_scheme) \ (material) \ - (phasefield) \ + (phasefield) \ + (constitutive_law) \ (mesh) \ (model) \ (model_solver) \ (neighborhood) \ (neighborhoods) \ (non_linear_solver) \ (non_local) \ (rules) \ (solver) \ (time_step_solver) \ (user) \ (weight_function) \ (contact_detector) \ (contact_resolution) \ (not_defined) // clang-format on /// Defines the possible section types AKANTU_CLASS_ENUM_DECLARE(ParserType, AKANTU_SECTION_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) #else enum class ParserType { cohesive_inserter, contact, embedded_interface, friction, global, heat, integration_scheme, material, phasefield, mesh, model, model_solver, neighborhood, neighborhoods, non_linear_solver, non_local, rules, solver, time_step_solver, user, weight_function, not_defined }; #endif /// Defines the possible search contexts/scopes (for parameter search) enum ParserParameterSearchCxt { _ppsc_current_scope = 0x1, _ppsc_parent_scope = 0x2, _ppsc_current_and_parent_scope = 0x3 }; /* ------------------------------------------------------------------------ */ /* Parameters Class */ /* ------------------------------------------------------------------------ */ class ParserSection; /// @brief The ParserParameter objects represent the end of tree branches as /// they /// are the different informations contained in the input file. class ParserParameter { public: ParserParameter() : name(std::string()), value(std::string()), dbg_filename(std::string()) { } ParserParameter(const std::string & name, const std::string & value, const ParserSection & parent_section) : parent_section(&parent_section), name(name), value(value), dbg_filename(std::string()) {} ParserParameter(const ParserParameter & param) = default; virtual ~ParserParameter() = default; /// Get parameter name const std::string & getName() const { return name; } /// Get parameter value const std::string & getValue() const { return value; } /// Set info for debug output void setDebugInfo(const std::string & filename, UInt line, UInt column) { dbg_filename = filename; dbg_line = line; dbg_column = column; } template <typename T> inline operator T() const; // template <typename T> inline operator Vector<T>() const; // template <typename T> inline operator Matrix<T>() const; /// Print parameter info in stream void printself(std::ostream & stream, __attribute__((unused)) unsigned int indent = 0) const { stream << name << ": " << value << " (" << dbg_filename << ":" << dbg_line << ":" << dbg_column << ")"; } private: void setParent(const ParserSection & sect) { parent_section = § } friend class ParserSection; private: /// Pointer to the parent section const ParserSection * parent_section{nullptr}; /// Name of the parameter std::string name; /// Value of the parameter std::string value; /// File for debug output std::string dbg_filename; /// Position of parameter in parsed file UInt dbg_line, dbg_column; }; /* ------------------------------------------------------------------------ */ /* Sections Class */ /* ------------------------------------------------------------------------ */ /// ParserSection represents a branch of the parsing tree. class ParserSection { public: using SubSections = std::multimap<ParserType, ParserSection>; using Parameters = std::map<std::string, ParserParameter>; private: using const_section_iterator_ = SubSections::const_iterator; public: /* ------------------------------------------------------------------------ */ /* SubSection iterator */ /* ------------------------------------------------------------------------ */ /// Iterator on sections class const_section_iterator { public: using iterator_category = std::forward_iterator_tag; using value_type = ParserSection; using pointer = ParserSection *; using reference = ParserSection &; const_section_iterator() = default; const_section_iterator(const const_section_iterator_ & it) : it(it) {} const_section_iterator(const const_section_iterator & other) = default; const_section_iterator & operator=(const const_section_iterator & other) = default; const ParserSection & operator*() const { return it->second; } const ParserSection * operator->() const { return &(it->second); } bool operator==(const const_section_iterator & other) const { return it == other.it; } bool operator!=(const const_section_iterator & other) const { return it != other.it; } const_section_iterator & operator++() { ++it; return *this; } const_section_iterator operator++(int) { const_section_iterator tmp = *this; operator++(); return tmp; } private: const_section_iterator_ it; }; /* ------------------------------------------------------------------------ */ /* Parameters iterator */ /* ------------------------------------------------------------------------ */ /// Iterator on parameters class const_parameter_iterator { public: const_parameter_iterator(const const_parameter_iterator & other) = default; const_parameter_iterator(const Parameters::const_iterator & it) : it(it) {} const_parameter_iterator & operator=(const const_parameter_iterator & other) { if (this != &other) { it = other.it; } return *this; } const ParserParameter & operator*() const { return it->second; } const ParserParameter * operator->() { return &(it->second); }; bool operator==(const const_parameter_iterator & other) const { return it == other.it; } bool operator!=(const const_parameter_iterator & other) const { return it != other.it; } const_parameter_iterator & operator++() { ++it; return *this; } const_parameter_iterator operator++(int) { const_parameter_iterator tmp = *this; operator++(); return tmp; } private: Parameters::const_iterator it; }; /* ---------------------------------------------------------------------- */ ParserSection() : name(std::string()) {} ParserSection(const std::string & name, ParserType type) : name(name), type(type) {} ParserSection(const std::string & name, ParserType type, const std::string & option, const ParserSection & parent_section) : parent_section(&parent_section), name(name), type(type), option(option) {} ParserSection(const ParserSection & section) : parent_section(section.parent_section), name(section.name), type(section.type), option(section.option), parameters(section.parameters), sub_sections_by_type(section.sub_sections_by_type) { setChldrenPointers(); } ParserSection & operator=(const ParserSection & other) { if (&other != this) { parent_section = other.parent_section; name = other.name; type = other.type; option = other.option; parameters = other.parameters; sub_sections_by_type = other.sub_sections_by_type; setChldrenPointers(); } return *this; } virtual ~ParserSection(); virtual void printself(std::ostream & stream, unsigned int indent = 0) const; /* ---------------------------------------------------------------------- */ /* Creation functions */ /* ---------------------------------------------------------------------- */ public: ParserParameter & addParameter(const ParserParameter & param); ParserSection & addSubSection(const ParserSection & section); protected: /// Clean ParserSection content void clean() { parameters.clear(); sub_sections_by_type.clear(); } private: void setChldrenPointers() { for (auto && param_pair : this->parameters) { param_pair.second.setParent(*this); } for (auto && sub_sect_pair : this->sub_sections_by_type) { sub_sect_pair.second.setParent(*this); } } /* ---------------------------------------------------------------------- */ /* Accessors */ /* ---------------------------------------------------------------------- */ public: class SubSectionsRange : public std::pair<const_section_iterator, const_section_iterator> { public: SubSectionsRange(const const_section_iterator & first, const const_section_iterator & second) : std::pair<const_section_iterator, const_section_iterator>(first, second) {} auto begin() { return this->first; } auto end() { return this->second; } }; /// Get begin and end iterators on subsections of certain type auto getSubSections(ParserType type = ParserType::_not_defined) const { if (type != ParserType::_not_defined) { auto range = sub_sections_by_type.equal_range(type); return SubSectionsRange(range.first, range.second); } return SubSectionsRange(sub_sections_by_type.begin(), sub_sections_by_type.end()); } /// Get number of subsections of certain type UInt getNbSubSections(ParserType type = ParserType::_not_defined) const { if (type != ParserType::_not_defined) { return this->sub_sections_by_type.count(type); } return this->sub_sections_by_type.size(); } /// Get begin and end iterators on parameters auto getParameters() const { return std::pair<const_parameter_iterator, const_parameter_iterator>( parameters.begin(), parameters.end()); } /* ---------------------------------------------------------------------- */ /// Get parameter within specified context const ParserParameter & getParameter( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if ((search_ctx & _ppsc_current_scope) != 0) { it = parameters.find(name); } if (it == parameters.end()) { if ((search_ctx & _ppsc_parent_scope) != 0 and parent_section != nullptr) { return parent_section->getParameter(name, search_ctx); } AKANTU_SILENT_EXCEPTION( "The parameter " << name << " has not been found in the specified context"); } return it->second; } /* ------------------------------------------------------------------------ */ /// Get parameter within specified context, with a default value in case the /// parameter does not exists template <class T> T getParameter( const std::string & name, const T & default_value, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { try { T tmp = this->getParameter(name, search_ctx); return tmp; } catch (debug::Exception &) { return default_value; } } /* ------------------------------------------------------------------------ */ /// Check if parameter exists within specified context bool hasParameter( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if ((search_ctx & _ppsc_current_scope) != 0) { it = parameters.find(name); } if (it == parameters.end()) { if ((search_ctx & _ppsc_parent_scope) != 0 and parent_section != nullptr) { return parent_section->hasParameter(name, search_ctx); } return false; } return true; } /* -------------------------------------------------------------------------- */ /// Get value of given parameter in context template <class T> T getParameterValue( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { const ParserParameter & tmp_param = getParameter(name, search_ctx); T t = tmp_param; return t; } /* -------------------------------------------------------------------------- */ /// Get section name std::string getName() const { return name; } /// Get section type ParserType getType() const { return type; } /// Get section option std::string getOption(const std::string & def = "") const { return (not option.empty()) ? option : def; } protected: void setParent(const ParserSection & sect) { parent_section = § } /* ---------------------------------------------------------------------- */ /* Members */ /* ---------------------------------------------------------------------- */ private: /// Pointer to the parent section const ParserSection * parent_section{nullptr}; /// Name of section std::string name; /// Type of section, see AKANTU_SECTION_TYPES ParserType type{ParserType::_not_defined}; /// Section option std::string option; /// Map of parameters in section Parameters parameters; /// Multi-map of subsections SubSections sub_sections_by_type; }; /* ------------------------------------------------------------------------ */ /* Parser Class */ /* ------------------------------------------------------------------------ */ /// Root of parsing tree, represents the global ParserSection class Parser : public ParserSection { public: Parser() : ParserSection("global", ParserType::_global) {} void parse(const std::string & filename); std::string getLastParsedFile() const; static bool isPermissive() { return permissive_parser; } public: /// Parse real scalar static Real parseReal(const std::string & value, const ParserSection & section); /// Parse real vector static Vector<Real> parseVector(const std::string & value, const ParserSection & section); /// Parse real matrix static Matrix<Real> parseMatrix(const std::string & value, const ParserSection & section); /// Parse real random parameter static RandomParameter<Real> parseRandomParameter(const std::string & value, const ParserSection & section); protected: /// General parse function template <class T, class Grammar> static T parseType(const std::string & value, Grammar & grammar); protected: // friend class Parsable; static bool permissive_parser; std::string last_parsed_file; }; inline std::ostream & operator<<(std::ostream & stream, const ParserParameter & _this) { _this.printself(stream); return stream; } inline std::ostream & operator<<(std::ostream & stream, const ParserSection & section) { section.printself(stream); return stream; } } // namespace akantu namespace std { template <> struct iterator_traits<::akantu::Parser::const_section_iterator> { using iterator_category = input_iterator_tag; using value_type = ::akantu::ParserParameter; using difference_type = ptrdiff_t; using pointer = const ::akantu::ParserParameter *; using reference = const ::akantu::ParserParameter &; }; } // namespace std #include "parser_tmpl.hh" #endif /* AKANTU_PARSER_HH_ */ diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/common/constitutive_laws/constitutive_law_non_local_interface.hh similarity index 67% copy from src/model/solid_mechanics/materials/material_non_local.hh copy to src/model/common/constitutive_laws/constitutive_law_non_local_interface.hh index a5239c9b4..c58d22320 100644 --- a/src/model/solid_mechanics/materials/material_non_local.hh +++ b/src/model/common/constitutive_laws/constitutive_law_non_local_interface.hh @@ -1,119 +1,119 @@ /** - * @file material_non_local.hh + * @file constitutive_law_non_local.hh * * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Fri Jun 18 2010 - * @date last modification: Fri Apr 09 2021 + * @date last modification: Mon Sep 11 2017 * - * @brief Material class that handle the non locality of a law for example - * damage. + * @brief ConstitutiveLaw class that handle the non locality of a law for + * example damage. * - * - * @section LICENSE - * - * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * 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 + * 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 + * 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 + * 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 + * 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 "material.hh" +#include "constitutive_law.hh" /* -------------------------------------------------------------------------- */ -#ifndef AKANTU_MATERIAL_NON_LOCAL_HH_ -#define AKANTU_MATERIAL_NON_LOCAL_HH_ +#ifndef AKANTU_CONSTITUTIVE_LAW_NON_LOCAL_HH_ +#define AKANTU_CONSTITUTIVE_LAW_NON_LOCAL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ -class MaterialNonLocalInterface { +class ConstitutiveLawNonLocalInterface { /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// initialize the material the non local parts of the material - void initMaterialNonLocal() { + /// initialize the constitutive_law the non local parts of the + /// constitutive_law + void initConstitutiveLawNonLocal() { this->registerNeighborhood(); this->registerNonLocalVariables(); }; /// insert the quadrature points in the neighborhoods of the non-local manager virtual void insertIntegrationPointsInNeighborhoods( GhostType ghost_type, const ElementTypeMapReal & quadrature_points_coordinates) = 0; /// update the values in the non-local internal fields virtual void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened, const ID & field_id, GhostType ghost_type, ElementKind kind) = 0; - /// constitutive law - virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0; protected: - /// get the name of the neighborhood for this material + /// get the name of the neighborhood for this constitutive_law virtual ID getNeighborhoodName() = 0; - /// register the neighborhoods for the material + /// register the neighborhoods for the constitutive_law virtual void registerNeighborhood() = 0; /// register the non local internal variable virtual void registerNonLocalVariables() = 0; virtual inline void onElementsAdded(const Array<Element> & /*unused*/, const NewElementsEvent & /*unused*/) {} }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ -template <UInt dim, class LocalParent> -class MaterialNonLocal : public MaterialNonLocalInterface, public LocalParent { +template <UInt dim, class ConstitutiveLawNonLocalInterface, class ConstitutiveLawParent> +class ConstitutiveLawNonLocal : public ConstitutiveLawNonLocalInterface, + public ConstitutiveLawParent { + /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - explicit MaterialNonLocal(SolidMechanicsModel & model, const ID & id); + explicit ConstitutiveLawNonLocal( + typename ConstitutiveLawParent::ConstitutiveLawsHandler & handler, + const ID & id); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// insert the quadrature points in the neighborhoods of the non-local manager void insertIntegrationPointsInNeighborhoods( GhostType ghost_type, const ElementTypeMapReal & quadrature_points_coordinates) override; /// update the values in the non-local internal fields void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened, const ID & field_id, GhostType ghost_type, ElementKind kind) override; - /// register the neighborhoods for the material + /// register the neighborhoods for the constitutive_law void registerNeighborhood() override; protected: - /// get the name of the neighborhood for this material + /// get the name of the neighborhood for this constitutive_law ID getNeighborhoodName() override { return this->name; } }; } // namespace akantu /* -------------------------------------------------------------------------- */ -#include "material_non_local_tmpl.hh" +#include "constitutive_law_non_local_interface_tmpl.hh" -#endif /* AKANTU_MATERIAL_NON_LOCAL_HH_ */ +#endif /* AKANTU_CONSTITUTIVE_LAW_NON_LOCAL_HH_ */ diff --git a/src/model/solid_mechanics/materials/material_non_local_tmpl.hh b/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh similarity index 77% rename from src/model/solid_mechanics/materials/material_non_local_tmpl.hh rename to src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh index 6c033166f..a28b6fa41 100644 --- a/src/model/solid_mechanics/materials/material_non_local_tmpl.hh +++ b/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh @@ -1,129 +1,128 @@ /** - * @file material_non_local_tmpl.hh + * @file constitutive_law_non_local_tmpl.hh * * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Thu Jul 06 2017 * @date last modification: Fri Mar 26 2021 * - * @brief Implementation of material non-local + * @brief Implementation of constitutive_law non-local * * * @section LICENSE * * Copyright (©) 2016-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 "material.hh" -#include "material_non_local.hh" +#include "constitutive_law.hh" +#include "constitutive_law_non_local_interface.hh" #include "non_local_neighborhood.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template <UInt dim, class LocalParent> -MaterialNonLocal<dim, LocalParent>::MaterialNonLocal( - SolidMechanicsModel & model, const ID & id) - : LocalParent(model, id) { - AKANTU_DEBUG_IN(); - - AKANTU_DEBUG_OUT(); +template <UInt dim, class ConstitutiveLawParent> +ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::ConstitutiveLawNonLocal( + typename ConstitutiveLawParent::ConstitutiveLawsHandler & handler, + const ID & id) + : ConstitutiveLawParent(handler, id) { } /* -------------------------------------------------------------------------- */ -template <UInt dim, class LocalParent> -void MaterialNonLocal<dim, LocalParent>::insertIntegrationPointsInNeighborhoods( - GhostType ghost_type, - const ElementTypeMapReal & quadrature_points_coordinates) { +template <UInt dim, class ConstitutiveLawParent> +void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>:: + insertIntegrationPointsInNeighborhoods( + GhostType ghost_type, + const ElementTypeMapReal & quadrature_points_coordinates) { IntegrationPoint q; q.ghost_type = ghost_type; auto & neighborhood = this->model.getNonLocalManager().getNeighborhood( this->getNeighborhoodName()); for (auto & type : this->element_filter.elementTypes(dim, ghost_type, _ek_regular)) { q.type = type; const auto & elem_filter = this->element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element != 0U) { UInt nb_quad = this->getFEEngine().getNbIntegrationPoints(type, ghost_type); const auto & quads = quadrature_points_coordinates(type, ghost_type); auto nb_total_element = this->model.getMesh().getNbElement(type, ghost_type); auto quads_it = quads.begin_reinterpret(dim, nb_quad, nb_total_element); for (auto & elem : elem_filter) { Matrix<Real> quads = quads_it[elem]; q.element = elem; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.global_num = q.element * nb_quad + nq; neighborhood.insertIntegrationPoint(q, quads(nq)); } } } } } /* -------------------------------------------------------------------------- */ -template <UInt dim, class LocalParent> -void MaterialNonLocal<dim, LocalParent>::updateNonLocalInternals( +template <UInt dim, class ConstitutiveLawParent> +void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::updateNonLocalInternals( ElementTypeMapReal & non_local_flattened, const ID & field_id, GhostType ghost_type, ElementKind kind) { - /// loop over all types in the material + /// loop over all types in the constitutive_law for (auto & el_type : this->element_filter.elementTypes(dim, ghost_type, kind)) { Array<Real> & internal = this->template getInternal<Real>(field_id)(el_type, ghost_type); auto & internal_flat = non_local_flattened(el_type, ghost_type); auto nb_component = internal_flat.getNbComponent(); auto internal_it = internal.begin(nb_component); auto internal_flat_it = internal_flat.begin(nb_component); /// loop all elements for the given type const auto & filter = this->element_filter(el_type, ghost_type); UInt nb_quads = this->getFEEngine().getNbIntegrationPoints(el_type, ghost_type); for (auto & elem : filter) { for (UInt q = 0; q < nb_quads; ++q, ++internal_it) { UInt global_quad = elem * nb_quads + q; *internal_it = internal_flat_it[global_quad]; } } } } /* -------------------------------------------------------------------------- */ -template <UInt dim, class LocalParent> -void MaterialNonLocal<dim, LocalParent>::registerNeighborhood() { +template <UInt dim, class ConstitutiveLawParent> +void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::registerNeighborhood() { ID name = this->getNeighborhoodName(); - this->model.getNonLocalManager().registerNeighborhood(name, name); + this->handler.getNonLocalManager().registerNeighborhood(name, name); } } // namespace akantu diff --git a/src/model/common/constitutive_laws/constitutive_laws_handler.hh b/src/model/common/constitutive_laws/constitutive_laws_handler.hh new file mode 100644 index 000000000..9c6b8ef4f --- /dev/null +++ b/src/model/common/constitutive_laws/constitutive_laws_handler.hh @@ -0,0 +1,291 @@ +/** + * @file constitutive_laws_handler.hh + * + * @author Nicolas Richart <nicolas.richart@epfl.ch> + * @author Mohit Pundir <mohit.pundir@epfl.ch> + * + * @date creation ven mar 26 2021 + * + * @brief A handler for multiple consistutive laws + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see <http://www.gnu.org/licenses/>. + * + */ +/* -------------------------------------------------------------------------- */ +#include "constitutive_law.hh" +#include "mesh.hh" +#include "mesh_events.hh" +#include "non_local_manager_callback.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH +#define AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH + +namespace akantu { +template <class ConstitutiveLawsHandler> class ConstitutiveLawSelector; +class NonLocalManager; +} // namespace akantu + +/* -------------------------------------------------------------------------- */ +namespace akantu { + +template <class ConsistutiveLawsType> +class ConstitutiveLawsHandler : public NonLocalManagerCallback { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +public: + ConstitutiveLawsHandler(const Mesh & mesh, UInt spatial_dimension, + const ID & parent_id) + : constitutive_law_index("constitutive_law index", id), + constitutive_law_local_numbering("constitutive_law local numbering", + id) { + constitutive_law_selector = + std::make_shared<DefaultConstitutiveLawSelector>( + constitutive_law_index); + + constitutive_law_index.initialize(mesh, _element_kind = _ek_not_defined, + _default_value = UInt(-1), + _with_nb_element = true); + constitutive_law_local_numbering.initialize( + mesh, _element_kind = _ek_not_defined, _with_nb_element = true); + } + + virtual ~ConstitutiveLawsHandler(); + + class NewConstitutiveLawElementsEvent : public NewElementsEvent { + public: + AKANTU_GET_MACRO_NOT_CONST(ConstitutiveLawList, constitutive_law, + Array<UInt> &); + AKANTU_GET_MACRO(ConstitutiveLawList, constitutive_law, + const Array<UInt> &); + + protected: + Array<UInt> constitutive_law; + }; + + /* ------------------------------------------------------------------------ */ + /* ConstitutiveLaws */ + /* ------------------------------------------------------------------------ */ +public: + /// register an empty constitutive_law of a given type + auto & registerNewConstitutiveLaw(const ID & cl_name, const ID & cl_type, + const ID & opt_param); + + /// reassigns constitutive_laws depending on the constitutive_law selector + virtual void reassignConstitutiveLaw(); + + template <class Func> void for_each_constitutive_law(Func && func) { + for (auto && constitutive_law : constitutive_laws) { + std::forward<Func>(func)( + aka::as_type<ConstitutiveLawType>(*constitutive_laws)); + } + } + +protected: + /// register a constitutive_law in the dynamic database + auto & registerNewConstitutiveLaw(const ParserSection & cl_section); + + /// read the constitutive_law files to instantiate all the constitutive_laws + void instantiateConstitutiveLaws(); + + /// set the element_id_by_constitutive_law and add the elements to the good + /// constitutive_laws + virtual void assignConstitutiveLawToElements( + const ElementTypeMapArray<UInt> * filter = nullptr); + +protected: + void splitElementByConstitutiveLaw( + const Array<Element> & elements, + std::vector<Array<Element>> & elements_per_cl) const; + + template <typename Operation> + void splitByConstitutiveLaw(const Array<Element> & elements, + Operation && op) const; + + /* ------------------------------------------------------------------------ */ + /* Dumpable interface (kept for convenience) and dumper relative functions */ + /* ------------------------------------------------------------------------ */ +public: + //! decide wether a field is a constitutive_law internal or not + bool isInternal(const std::string & field_name, ElementKind element_kind); + + //! give the amount of data per element + virtual ElementTypeMap<UInt> + getInternalDataPerElem(const std::string & field_name, ElementKind kind); + + //! flatten a given constitutive_law internal field + ElementTypeMapArray<Real> & + flattenInternal(const std::string & field_name, ElementKind kind, + GhostType ghost_type = _not_ghost); + + //! flatten all the registered constitutive_law internals + void flattenAllRegisteredInternals(ElementKind kind); + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + /// get an iterable on the constitutive_laws + inline decltype(auto) getConstitutiveLaws(); + + /// get an iterable on the constitutive_laws + inline decltype(auto) getConstitutiveLaws() const; + + /// get a particular constitutive_law (by numerical constitutive_law index) + inline auto & getConstitutiveLaw(UInt cl_index); + + /// get a particular constitutive_law (by numerical constitutive_law index) + inline const auto & getConstitutiveLaw(UInt cl_index) const; + + /// get a particular constitutive_law (by constitutive_law name) + inline auto & getConstitutiveLaw(const std::string & name); + + /// get a particular constitutive_law (by constitutive_law name) + inline const auto & getConstitutiveLaw(const std::string & name) const; + + /// get a particular constitutive_law id from is name + inline UInt getConstitutiveLawIndex(const std::string & name) const; + + /// give the number of constitutive_laws + inline UInt getNbConstitutiveLaws() const { return constitutive_laws.size(); } + + /// give the constitutive_law internal index from its id + Int getInternalIndexFromID(const ID & id) const; + + AKANTU_GET_MACRO(ConstitutiveLawByElement, constitutive_law_index, + const ElementTypeMapArray<UInt> &); + AKANTU_GET_MACRO(ConstitutiveLawLocalNumbering, + constitutive_law_local_numbering, + const ElementTypeMapArray<UInt> &); + + /// vectors containing local constitutive_law element index for each global + /// element index + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConstitutiveLawByElement, + constitutive_law_index, UInt); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConstitutiveLawLocalNumbering, + constitutive_law_local_numbering, + UInt); + + AKANTU_GET_MACRO_NOT_CONST(ConstitutiveLawSelector, + *constitutive_law_selector, + ConstitutiveLawSelector &); + + /// Access the non_local_manager interface + AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); + + void setConstitutiveLawSelector( + std::shared_ptr<ConstitutiveLawSelector> constitutive_law_selector) { + this->constitutive_law_selector = std::move(constitutive_law_selector); + } + + /* ------------------------------------------------------------------------ */ + /* 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; + + /* ------------------------------------------------------------------------ */ + /* Mesh Event Handler inherited members */ + /* ------------------------------------------------------------------------ */ +protected: + void onElementsAdded(const Array<Element> & element_list, + const NewElementsEvent & event) override; + void onElementsRemoved(const Array<Element> & element_list, + const ElementTypeMapArray<UInt> & new_numbering, + const RemovedElementsEvent & event) override; + + /* ------------------------------------------------------------------------ */ + /* NonLocalManager inherited members */ + /* ------------------------------------------------------------------------ */ +protected: + void initializeNonLocal() override; + + void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; + + void computeNonLocalStresses(GhostType ghost_type) override; + + void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; + + /// update the values of the non local internal + void updateLocalInternal(ElementTypeMapReal & internal_flat, + GhostType ghost_type, ElementKind kind) override; + + /// copy the results of the averaging in the constitutive_laws + void updateNonLocalInternal(ElementTypeMapReal & internal_flat, + GhostType ghost_type, ElementKind kind) override; + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +private: + /// id of the derived class + const ID & parent_id; + + /// underlying mesh + const Mesh & _mesh; + + /// spatial_dimension of the derived model + UInt _spatial_dimension{_all_dimensions}; + + /// mapping between constitutive_law name and constitutive_law internal id + std::map<std::string, UInt> constitutive_laws_names_to_id; + + /// Arrays containing the constitutive_law index for each element + ElementTypeMapArray<UInt> constitutive_law_index; + + /// Arrays containing the position in the element filter of the + /// constitutive_law (constitutive_law's local numbering) + ElementTypeMapArray<UInt> constitutive_law_local_numbering; + + /// list of used constitutive_laws + std::vector<std::unique_ptr<ConstitutiveLaw>> constitutive_laws; + + /// class defining of to choose a constitutive_law + std::shared_ptr<ConstitutiveLawSelector> constitutive_law_selector; + + using flatten_internal_map = + std::map<std::pair<std::string, ElementKind>, + std::unique_ptr<ElementTypeMapArray<Real>>>; + + /// map a registered internals to be flattened for dump purposes + flatten_internal_map registered_internals; + + /// tells if the constitutive_law are instantiated + bool are_constitutive_laws_instantiated{false}; + + /// non local manager + std::unique_ptr<NonLocalManager> non_local_manager; + + template <class ConstitutiveLawsHandler> friend class ConstitutiveLaw; +}; + +} // namespace akantu +#endif /* AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH */ diff --git a/src/model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh b/src/model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh new file mode 100644 index 000000000..b0b6052e9 --- /dev/null +++ b/src/model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh @@ -0,0 +1,549 @@ +/** + * @file constitutive_laws_handler_tmpl.hh + * + * @author Nicolas Richart <nicolas.richart@epfl.ch> + * @author Mohit Pundir <mohit.pundir@epfl.ch> + * + * @date creation ven mar 26 2021 + * + * @brief Implementation of the consistutive laws handler. + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see <http://www.gnu.org/licenses/>. + * + */ +/* -------------------------------------------------------------------------- */ +#include "constitutive_law_non_local_interface.hh" +#include "constitutive_laws_handler.hh" +#include "non_local_manager.hh" +#include "parsable.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef AKANTU_CONSTITUTIVE_LAWS_HANDLER_TMPL_HH +#define AKANTU_CONSTITUTIVE_LAWS_HANDLER_TMPL_HH + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +inline decltype(auto) +ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaws() { + return make_transform_adaptor( + constitutive_laws, [](auto && value) -> decltype(auto) { + return aka::as_type<ConstitutiveLawType>(*value); + }); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +inline decltype(auto) +ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaws() const { + return make_transform_adaptor( + constitutive_laws, [](auto && value) -> decltype(auto) { + return aka::as_type<ConstitutiveLawType>(*value); + }); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +inline UInt +ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLawIndex( + const std::string & name) const { + auto it = constitutive_laws_names_to_id.find(name); + if (it == constitutive_laws_names_to_id.end()) { + AKANTU_SILENT_EXCEPTION( + "The model " << parent_id << " has no constitutive_law named " << name); + } + + return it->second; +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +inline auto & ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw( + UInt cl_index) { + AKANTU_DEBUG_ASSERT(cl_index < constitutive_laws.size(), + "The model " << parent_id + << " has no constitutive law number " + << cl_index); + return aka::as_type<ConstitutiveLawType>(*constitutive_laws.at(cl_index)); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +inline const auto & +ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw( + UInt cl_index) const { + AKANTU_DEBUG_ASSERT(cl_index < constitutive_laws.size(), + "The model " << parent_id + << " has no constitutive law number " + << cl_index); + return aka::as_type<ConstitutiveLawType>(*constitutive_laws.at(cl_index)); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +inline auto & ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw( + const std::string & name) { + std::map<std::string, UInt>::const_iterator it = + constitutive_laws_names_to_id.find(name); + if (it == constitutive_laws_names_to_id.end()) { + AKANTU_SILENT_EXCEPTION( + "The model " << parent_id << " has no constitutive_law named " << name); + } + + return aka::as_type<ConstitutiveLawType>(*constitutive_laws[it->second]); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +inline const auto & +ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw( + const std::string & name) const { + auto it = constitutive_laws_names_to_id.find(name); + if (it == constitutive_laws_names_to_id.end()) { + AKANTU_SILENT_EXCEPTION( + "The model " << parent_id << " has no constitutive_law named " << name); + } + return aka::as_type<ConstitutiveLawType>(*constitutive_laws[it->second]); +} + +/* -------------------------------------------------------------------------- */ +#define AKA_FWD_CL(...) ::std::forward<decltype(__VA_ARGS__)>(__VA_ARGS__) + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +template <typename Operation> +void ConstitutiveLawsHandler<ConstitutiveLawType>::splitByConstitutiveLaw( + const Array<Element> & elements, Operation && op) const { + std::vector<Array<Element>> elements_per_cl(constitutive_laws.size()); + this->splitElementByConstitutiveLaw(elements, elements_per_cl); + + for (auto && cl : zip(constitutive_laws, elements_per_cl)) { + AKA_FWD_CL(op)(AKA_FWD_CL(*std::get<0>(cl)), AKA_FWD_CL(std::get<1>(cl))); + } +} + +#undef AKA_FWD_CL + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +auto & ConstitutiveLawsHandler<ConstitutiveLawType>::registerNewConstitutiveLaw( + const ParserSection & section) { + std::string cl_name; + std::string cl_type = section.getName(); + std::string opt_param = section.getOption(); + + try { + std::string tmp = section.getParameter("name"); + cl_name = tmp; /** this can seam 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 constitutive_law of type \'" + << cl_type + << "\' in the input file has been defined without a name!"); + } + + auto && cl = this->registerNewConstitutiveLaw(cl_name, cl_type, opt_param); + + cl.parseSection(section); + + return cl; +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +auto & ConstitutiveLawsHandler<ConstitutiveLawType>::registerNewConstitutiveLaw( + const ID & cl_name, const ID & cl_type, const ID & opt_param) { + AKANTU_DEBUG_ASSERT(constitutive_laws_names_to_id.find(cl_name) == + constitutive_laws_names_to_id.end(), + "A constitutive_law with this name '" + << cl_name << "' has already been registered. " + << "Please use unique names for constitutive_laws"); + + UInt cl_count = constitutive_laws.size(); + constitutive_laws_names_to_id[cl_name] = cl_count; + + std::stringstream sstr_cl; + sstr_cl << this->id << ":" << cl_count << ":" << cl_type; + ID cl_id = sstr_cl.str(); + + std::unique_ptr<ConstitutiveLaw> constitutive_law = + ConstitutiveLawType::getFactory().allocate(cl_type, _spatial_dimension, + opt_param, *this, cl_id); + + constitutive_laws.push_back(std::move(constitutive_law)); + + return aka::as_type<ConstitutiveLawType>(*(constitutive_laws.back())); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler< + ConstitutiveLawType>::instantiateConstitutiveLaws() { + ParserSection model_section; + bool is_empty; + std::tie(model_section, is_empty) = this->getParserSection(); + + if (not is_empty) { + auto model_constitutive_laws = + model_section.getSubSections(ParserType::_constitutive_law); + for (const auto & section : model_constitutive_laws) { + this->registerNewConstitutiveLaw(section); + } + } + + auto sub_sections = + this->parser.getSubSections(ParserType::_constitutive_law); + for (const auto & section : sub_sections) { + this->registerNewConstitutiveLaw(section); + } + +#ifdef AKANTU_DAMAGE_NON_LOCAL + for (auto & constitutive_law : constitutive_laws) { + if (dynamic_cast<ConstitutiveLawNonLocalInterface *>( + constitutive_law.get()) == nullptr) { + continue; + } + + this->non_local_manager = std::make_unique<NonLocalManager>( + *this, *this, parent_id + ":non_local_manager"); + break; + } +#endif + + if (constitutive_laws.empty()) { + AKANTU_EXCEPTION("No constitutive_laws where instantiated for the model" + << parent_id); + } + are_constitutive_laws_instantiated = true; +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>:: + assignConstitutiveLawToElements(const ElementTypeMapArray<UInt> * filter) { + + for_each_element( + _mesh, + [&](auto && element) { + UInt cl_index = (*constitutive_law_selector)(element); + AKANTU_DEBUG_ASSERT(cl_index < constitutive_laws.size(), + "The constitutive_law selector returned an index " + "that does not exists"); + constitutive_law_index(element) = cl_index; + }, + _element_filter = filter, _ghost_type = _not_ghost); + + if (non_local_manager) { + non_local_manager->synchronize(*this, + SynchronizationTag::_constitutive_law_id); + } + + for_each_element( + _mesh, + [&](auto && element) { + auto cl_index = constitutive_law_index(element); + auto index = constitutive_laws[cl_index]->addElement(element); + constitutive_law_local_numbering(element) = index; + }, + _element_filter = filter, _ghost_type = _not_ghost); + + // synchronize the element constitutive_law arrays + this->synchronize(SynchronizationTag::_constitutive_law_id); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::initConstitutiveLaws() { + AKANTU_DEBUG_ASSERT(not constitutive_laws.empty(), + "No constitutive_law to initialize !"); + + this->assignConstitutiveLawToElements(); + + for (auto & constitutive_law : constitutive_laws) { + /// init internals properties + constitutive_law->initConstitutiveLaw(); + } + + this->synchronize(SynchronizationTag::_smm_init_cl); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +Int ConstitutiveLawsHandler<ConstitutiveLawType>::getInternalIndexFromID( + const ID & id) const { + AKANTU_DEBUG_IN(); + + auto it = constitutive_laws.begin(); + auto end = constitutive_laws.end(); + + for (; it != end; ++it) { + if ((*it)->getID() == id) { + AKANTU_DEBUG_OUT(); + return (it - constitutive_laws.begin()); + } + } + + AKANTU_DEBUG_OUT(); + return -1; +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::reassignConstitutiveLaw() { + AKANTU_DEBUG_IN(); + + std::vector<Array<Element>> element_to_add(constitutive_laws.size()); + std::vector<Array<Element>> element_to_remove(constitutive_laws.size()); + + for_each_element(mesh, [&](auto && element) { + auto old_constitutive_law = constitutive_law_index(element); + auto new_constitutive_law = (*constitutive_law_selector)(element); + if (old_constitutive_law != new_constitutive_law) { + element_to_add[new_constitutive_law].push_back(element); + element_to_remove[old_constitutive_law].push_back(element); + } + }); + + for (auto && data : enumerate(constitutive_laws)) { + auto cl_index = std::get<0>(data); + auto & cl = *std::get<1>(data); + + cl.removeElements(element_to_remove[cl_index]); + cl.addElements(element_to_add[cl_index]); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>:: + updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) { + const ID field_name = criterion.getName(); + for (auto & constitutive_law : constitutive_laws) { + if (!constitutive_law->isInternal<Real>(field_name, _ek_regular)) { + continue; + } + + for (auto ghost_type : ghost_types) { + constitutive_law->flattenInternal(field_name, criterion, ghost_type, + _ek_regular); + } + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::initializeNonLocal() { + this->non_local_manager->synchronize( + *this, SynchronizationTag::_constitutive_law_id); +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>:: + insertIntegrationPointsInNeighborhoods(GhostType ghost_type) { + for (auto & cl : constitutive_laws) { + ConstitutiveLawNonLocalInterface * cl_non_local; + if ((cl_non_local = dynamic_cast<ConstitutiveLawNonLocalInterface *>( + cl.get())) == nullptr) { + continue; + } + + ElementTypeMapArray<Real> quadrature_points_coordinates( + "quadrature_points_coordinates_tmp_nl", this->id); + quadrature_points_coordinates.initialize(this->getFEEngine(), + _nb_component = spatial_dimension, + _ghost_type = ghost_type); + + for (const auto & type : quadrature_points_coordinates.elementTypes( + Model::spatial_dimension, ghost_type)) { + this->getFEEngine().computeIntegrationPointsCoordinates( + quadrature_points_coordinates(type, ghost_type), type, ghost_type); + } + + cl_non_local->initConstitutiveLawNonLocal(); + + cl_non_local->insertIntegrationPointsInNeighborhoods( + ghost_type, quadrature_points_coordinates); + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::updateLocalInternal( + ElementTypeMapReal & internal_flat, GhostType ghost_type, + ElementKind kind) { + const ID field_name = internal_flat.getName(); + for (auto & constitutive_law : constitutive_laws) { + if (constitutive_law->isInternal<Real>(field_name, kind)) { + constitutive_law->flattenInternal(field_name, internal_flat, ghost_type, + kind); + } + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::updateNonLocalInternal( + ElementTypeMapReal & internal_flat, GhostType ghost_type, + ElementKind kind) { + + const ID field_name = internal_flat.getName(); + + for (auto & mat : constitutive_laws) { + if (not aka::is_of_type<ConstitutiveLawNonLocalInterface>(*mat)) { + continue; + } + + auto & mat_non_local = + dynamic_cast<ConstitutiveLawNonLocalInterface &>(*mat); + mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, + kind); + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>:: + splitElementByConstitutiveLaw( + const Array<Element> & elements, + std::vector<Array<Element>> & elements_per_cl) const { + for (const auto & el : elements) { + Element cl_el = el; + cl_el.element = this->constitutive_law_local_numbering(el); + elements_per_cl[this->constitutive_law_index(el)].push_back(cl_el); + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +UInt ConstitutiveLawsHandler<ConstitutiveLawType>::getNbData( + const Array<Element> & elements, const SynchronizationTag & tag) const { + UInt size{0}; + + if(tag == SynchronizationTag::_constitutive_law_id) { + size += elements.size() * sizeof(UInt); + break; + } + + + if (tag != SynchronizationTag::_constitutive_law_id) { + splitByConstitutiveLaw(elements, [&](auto && cl, auto && elements) { + size += cl.getNbData(elements, tag); + }); + } + return size; +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::packData( + CommunicationBuffer & buffer, const Array<Element> & elements, + const SynchronizationTag & tag) const { + if (tag == SynchronizationTag::_constitutive_law_id) { + DataAccessor<Element>::packElementalDataHelper( + constitutive_law_index, buffer, elements, false, getFEEngine()); + } + + if (tag != SynchronizationTag::_constitutive_law_id) { + splitByConstitutiveLaw(elements, [&](auto && cl, auto && elements) { + cl.packData(buffer, elements, tag); + }); + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::unpackData( + CommunicationBuffer & buffer, const Array<Element> & elements, + const SynchronizationTag & tag) { + if (tag == SynchronizationTag::_constitutive_law_id) { + for (auto && element : elements) { + UInt recv_cl_index; + buffer >> recv_cl_index; + UInt & cl_index = constitutive_law_index(element); + if (cl_index != UInt(-1)) { + continue; + } + + // add ghosts element to the correct constitutive_law + cl_index = recv_cl_index; + UInt index = constitutive_laws[cl_index]->addElement(element); + constitutive_law_local_numbering(element) = index; + } + } + + if (tag != SynchronizationTag::_constitutive_law_id) { + splitByConstitutiveLaw(elements, [&](auto && cl, auto && elements) { + cl.unpackData(buffer, elements, tag); + }); + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::onElementsAdded( + const Array<Element> & element_list, const NewElementsEvent & event) { + this->constitutive_law_index.initialize(mesh, _element_kind = _ek_not_defined, + _with_nb_element = true, + _default_value = UInt(-1)); + this->constitutive_law_local_numbering.initialize( + mesh, _element_kind = _ek_not_defined, _with_nb_element = true, + _default_value = UInt(-1)); + + ElementTypeMapArray<UInt> filter("new_element_filter", this->getID()); + + for (const auto & elem : element_list) { + if (mesh.getSpatialDimension(elem.type) != spatial_dimension) { + continue; + } + + if (!filter.exists(elem.type, elem.ghost_type)) { + filter.alloc(0, 1, elem.type, elem.ghost_type); + } + filter(elem.type, elem.ghost_type).push_back(elem.element); + } + + // this fails in parallel if the event is sent on facet between constructor + // and initFull \todo: to debug... + this->assignConstitutiveLawToElements(&filter); + + for (auto & constitutive_law : constitutive_laws) { + constitutive_law->onElementsAdded(element_list, event); + } +} + +/* -------------------------------------------------------------------------- */ +template <class ConstitutiveLawType> +void ConstitutiveLawsHandler<ConstitutiveLawType>::onElementsRemoved( + const Array<Element> & element_list, + const ElementTypeMapArray<UInt> & new_numbering, + const RemovedElementsEvent & event) { + for (auto & constitutive_law : constitutive_laws) { + constitutive_law->onElementsRemoved(element_list, new_numbering, event); + } +} + +} // namespace akantu + +#endif /* AKANTU_CONSTITUTIVE_LAWS_HANDLER_TMPL_HH */ diff --git a/src/model/common/consitutive_laws/internal_field.hh b/src/model/common/constitutive_laws/internal_field.hh similarity index 100% rename from src/model/common/consitutive_laws/internal_field.hh rename to src/model/common/constitutive_laws/internal_field.hh diff --git a/src/model/common/consitutive_laws/internal_field_tmpl.hh b/src/model/common/constitutive_laws/internal_field_tmpl.hh similarity index 100% rename from src/model/common/consitutive_laws/internal_field_tmpl.hh rename to src/model/common/constitutive_laws/internal_field_tmpl.hh diff --git a/src/model/common/consitutive_laws/random_internal_field.hh b/src/model/common/constitutive_laws/random_internal_field.hh similarity index 100% rename from src/model/common/consitutive_laws/random_internal_field.hh rename to src/model/common/constitutive_laws/random_internal_field.hh diff --git a/src/model/common/consitutive_laws/random_internal_field_tmpl.hh b/src/model/common/constitutive_laws/random_internal_field_tmpl.hh similarity index 100% rename from src/model/common/consitutive_laws/random_internal_field_tmpl.hh rename to src/model/common/constitutive_laws/random_internal_field_tmpl.hh diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/solid_mechanics/materials/material_non_local.hh index a5239c9b4..cca6e20ec 100644 --- a/src/model/solid_mechanics/materials/material_non_local.hh +++ b/src/model/solid_mechanics/materials/material_non_local.hh @@ -1,119 +1,59 @@ /** * @file material_non_local.hh * * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Apr 09 2021 * * @brief Material class that handle the non locality of a law for example * damage. * * * @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 "material.hh" +#include "constitutive_law_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_NON_LOCAL_HH_ #define AKANTU_MATERIAL_NON_LOCAL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ -class MaterialNonLocalInterface { +class MaterialNonLocalInterface : public ConstitutiveLawNonLocalInterface { /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// initialize the material the non local parts of the material - void initMaterialNonLocal() { - this->registerNeighborhood(); - this->registerNonLocalVariables(); - }; - - /// insert the quadrature points in the neighborhoods of the non-local manager - virtual void insertIntegrationPointsInNeighborhoods( - GhostType ghost_type, - const ElementTypeMapReal & quadrature_points_coordinates) = 0; - - /// update the values in the non-local internal fields - virtual void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened, - const ID & field_id, - GhostType ghost_type, - ElementKind kind) = 0; /// constitutive law virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0; - -protected: - /// get the name of the neighborhood for this material - virtual ID getNeighborhoodName() = 0; - - /// register the neighborhoods for the material - virtual void registerNeighborhood() = 0; - - /// register the non local internal variable - virtual void registerNonLocalVariables() = 0; - - virtual inline void onElementsAdded(const Array<Element> & /*unused*/, - const NewElementsEvent & /*unused*/) {} }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template <UInt dim, class LocalParent> -class MaterialNonLocal : public MaterialNonLocalInterface, public LocalParent { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - explicit MaterialNonLocal(SolidMechanicsModel & model, const ID & id); - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ -public: - /// insert the quadrature points in the neighborhoods of the non-local manager - void insertIntegrationPointsInNeighborhoods( - GhostType ghost_type, - const ElementTypeMapReal & quadrature_points_coordinates) override; - - /// update the values in the non-local internal fields - void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened, - const ID & field_id, GhostType ghost_type, - ElementKind kind) override; - - /// register the neighborhoods for the material - void registerNeighborhood() override; - -protected: - /// get the name of the neighborhood for this material - ID getNeighborhoodName() override { return this->name; } -}; - -} // namespace akantu - -/* -------------------------------------------------------------------------- */ -#include "material_non_local_tmpl.hh" +using MaterialNonLocal = + ConstitutiveLawNonLocal<dim, MaterialNonLocalInterface, LocalParent>; #endif /* AKANTU_MATERIAL_NON_LOCAL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 7136890bb..d83fd5211 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1259 +1,1108 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei <ramin.aghababaei@epfl.ch> * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Mauro Corrado <mauro.corrado@epfl.ch> * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@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 Clement Roux <clement.roux@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Apr 09 2021 * * @brief Implementation of the SolidMechanicsModel class * * * @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 "solid_mechanics_model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model_tmpl.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model * @param model_type this is an internal parameter for inheritance purposes */ - SolidMechanicsModel::SolidMechanicsModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr<DOFManager> dof_manager, const ModelType model_type) - : Model(mesh, model_type, dim, id), material_index("material index", id), - material_local_numbering("material local numbering", id) { + : Model(mesh, model_type, dim, id) { AKANTU_DEBUG_IN(); this->initDOFManager(dof_manager); this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); this->mesh.registerDumper<DumperParaview>("solid_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); - material_selector = std::make_shared<DefaultMaterialSelector>(material_index); + // material_selector = + // std::make_shared<DefaultMaterialSelector>(material_index); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); this->mesh.getDumper().setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param options * \parblock * contains the different options to initialize the model * \li \c analysis_method specify the type of solver to use * \endparblock */ void SolidMechanicsModel::initFullImpl(const ModelOptions & options) { - material_index.initialize(mesh, _element_kind = _ek_not_defined, - _default_value = UInt(-1), _with_nb_element = true); - material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, - _with_nb_element = true); - Model::initFullImpl(options); // initialize the materials if (not this->parser.getLastParsedFile().empty()) { this->instantiateMaterials(); this->initMaterials(); } this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { 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; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ std::tuple<ID, TimeStepSolverType> SolidMechanicsModel::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); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { auto & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->previous_displacement, spatial_dimension, "previous_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->blocked_dofs, spatial_dimension, "blocked_dofs"); this->allocNodalField(this->current_position, spatial_dimension, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->velocity, spatial_dimension, "velocity"); this->allocNodalField(this->acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForces(); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) const { // \TODO check the materials to know what is the correct answer if (matrix_id == "C") { return _mt_not_defined; } if (matrix_id == "K") { auto matrix_type = _unsymmetric; - for (auto & material : materials) { - matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id)); - } + for_each_constitutive_law([&](auto && material) { + matrix_type = std::max(matrix_type, material.getMatrixType(matrix_id)); + }); } - return _symmetric; +} +return _symmetric; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { this->assembleMass(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { this->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::beforeSolveStep() { for (auto & material : materials) { material->beforeSolveStep(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::afterSolveStep(bool converged) { for (auto & material : materials) { material->afterSolveStep(converged); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as \f$F_{int} = \int_{\Omega} N * \sigma d\Omega@\f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->zero(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); - for (auto & material : materials) { - material->computeAllStresses(_not_ghost); - } + for_each_constitutive_law( + [](auto && material) { material.computeAllStresses(_not_ghost); }); /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if (this->non_local_manager) { this->non_local_manager->computeAllNonLocalStresses(); } // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); - for (auto & material : materials) { - material->assembleInternalForces(_not_ghost); - } + for_each_constitutive_law( + [](auto && material) { material.assembleInternalForces(_not_ghost); }); // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); - for (auto & material : materials) { - material->assembleInternalForces(_ghost); - } + for_each_constitutive_law( + [](auto && material) { material.assembleInternalForces(_ghost); }); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix(bool need_to_reassemble) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); if (not this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", this->getMatrixType("K")); } // Check if materials need to recompute the matrix - for (auto & material : materials) { + for_each_constitutive_law([&](auto && material) { need_to_reassemble |= material->hasMatrixChanged("K"); - } + }); if (need_to_reassemble) { this->getDOFManager().getMatrix("K").zero(); // call compute stiffness matrix on each local elements - for (auto & material : materials) { - material->assembleStiffnessMatrix(_not_ghost); - } + for_each_constitutive_law( + [](auto && material) { material.assembleStiffnessMatrix(_not_ghost); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) { return; } this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array<Real> & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::updateDataForNonLocalCriterion( - ElementTypeMapReal & criterion) { - const ID field_name = criterion.getName(); - for (auto & material : materials) { - if (!material->isInternal<Real>(field_name, _ek_regular)) { - continue; - } - - for (auto ghost_type : ghost_types) { - material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); - } - } -} - /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(GhostType ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits<Real>::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); - auto mat_indexes = material_index(type, ghost_type).begin(); - auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); + auto mat_indexes = + this->getConstitutiveLawByElement(type, ghost_type).begin(); + auto mat_loc_num = + this->getConstitutiveLawLocalNumbering(type, ghost_type).begin(); Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { this->assembleLumpedMatrix("M"); auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const auto & v = *v_it; const auto & m = *m_it; Real mv2 = 0.; auto is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); auto count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits<Real>::epsilon()) { mv2 += v(i) * v(i) * m(i); } } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { this->assembleMatrix("M"); Array<Real> Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().assembleMatMulVectToArray("displacement", "M", *this->velocity, Mv); for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension), make_view(*this->velocity, spatial_dimension))) { ekin += std::get<2>(data).dot(std::get<1>(data)) * static_cast<Real>(mesh.isLocalOrMasterNode(std::get<0>(data))); } } else { AKANTU_ERROR("No function called to assemble the mass matrix."); } mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(ElementType type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array<UInt> filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); auto vit = vel_on_quad.begin(Model::spatial_dimension); auto vend = vel_on_quad.end(Model::spatial_dimension); Vector<Real> rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); decltype(ext_force_it) incr_or_velo_it; if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } else { incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); } Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) { work -= int_force(i) * incr_or_velo(i); } else { work += ext_force(i) * incr_or_velo(i); } } } } mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum); if (this->method != _static) { work *= this->getTimeStep(); } AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; - for (auto & material : materials) { - energy += material->getEnergy(energy_id); - } + for_each_constitutive_law( + [&](auto && material) { energy += material->getEnergy(energy_id); }); /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const ID & energy_id, const ID & group_id) { auto && group = mesh.getElementGroup(group_id); auto energy = 0.; for (auto && type : group.elementTypes()) { for (auto el : group.getElementsIterable(type)) { energy += getEnergy(energy_id, el); } } /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list, const NewElementsEvent & event) { - AKANTU_DEBUG_IN(); - - this->material_index.initialize(mesh, _element_kind = _ek_not_defined, - _with_nb_element = true, - _default_value = UInt(-1)); - this->material_local_numbering.initialize( - mesh, _element_kind = _ek_not_defined, _with_nb_element = true, - _default_value = UInt(-1)); - - ElementTypeMapArray<UInt> filter("new_element_filter", this->getID()); - - for (const auto & elem : element_list) { - if (mesh.getSpatialDimension(elem.type) != spatial_dimension) { - continue; - } - - if (!filter.exists(elem.type, elem.ghost_type)) { - filter.alloc(0, 1, elem.type, elem.ghost_type); - } - filter(elem.type, elem.ghost_type).push_back(elem.element); - } - - // this fails in parallel if the event is sent on facet between constructor - // and initFull \todo: to debug... - this->assignMaterialToElements(&filter); - - for (auto & material : materials) { - material->onElementsAdded(element_list, event); - } - - AKANTU_DEBUG_OUT(); + ConstitutiveLawsHandler<Material>::onElementsAdded(element_list, event); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( const Array<Element> & element_list, const ElementTypeMapArray<UInt> & new_numbering, const RemovedElementsEvent & event) { - for (auto & material : materials) { - material->onElementsRemoved(element_list, new_numbering, event); - } + ConstitutiveLawsHandler<Material>::onElementsRemoved(element_list, + new_numbering, event); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (displacement) { displacement->resize(nb_nodes, 0.); ++displacement_release; } if (mass) { mass->resize(nb_nodes, 0.); } if (velocity) { velocity->resize(nb_nodes, 0.); } if (acceleration) { acceleration->resize(nb_nodes, 0.); } if (external_force) { external_force->resize(nb_nodes, 0.); } if (internal_force) { internal_force->resize(nb_nodes, 0.); } if (blocked_dofs) { blocked_dofs->resize(nb_nodes, false); } if (current_position) { current_position->resize(nb_nodes, 0.); } if (previous_displacement) { previous_displacement->resize(nb_nodes, 0.); } if (displacement_increment) { displacement_increment->resize(nb_nodes, 0.); } for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } need_to_reassemble_lumped_mass = true; need_to_reassemble_mass = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(const Array<UInt> & /*element_list*/, const Array<UInt> & new_numbering, const RemovedNodesEvent & /*event*/) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) { mesh.removeNodesFromArray(*mass, new_numbering); } if (velocity) { mesh.removeNodesFromArray(*velocity, new_numbering); } if (acceleration) { mesh.removeNodesFromArray(*acceleration, new_numbering); } if (internal_force) { mesh.removeNodesFromArray(*internal_force, new_numbering); } if (external_force) { mesh.removeNodesFromArray(*external_force, new_numbering); } if (blocked_dofs) { mesh.removeNodesFromArray(*blocked_dofs, new_numbering); } // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) { mesh.removeNodesFromArray(*displacement_increment, new_numbering); } if (previous_displacement) { mesh.removeNodesFromArray(*previous_displacement, new_numbering); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Solid 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 << " ]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); if (velocity) { velocity->printself(stream, indent + 2); } if (acceleration) { acceleration->printself(stream, indent + 2); } if (mass) { mass->printself(stream, indent + 2); } external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + materials [" << std::endl; - for (const auto & material : materials) { - material->printself(stream, indent + 2); - } + for_each_constitutive_law( + [&](auto && material) { material.printself(stream, indent + 2); }); stream << space << " ]" << std::endl; stream << space << "]" << std::endl; } -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::initializeNonLocal() { - this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id); -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( - GhostType ghost_type) { - for (auto & mat : materials) { - MaterialNonLocalInterface * mat_non_local; - if ((mat_non_local = - dynamic_cast<MaterialNonLocalInterface *>(mat.get())) == nullptr) { - continue; - } - - ElementTypeMapArray<Real> quadrature_points_coordinates( - "quadrature_points_coordinates_tmp_nl", this->id); - quadrature_points_coordinates.initialize(this->getFEEngine(), - _nb_component = spatial_dimension, - _ghost_type = ghost_type); - - for (const auto & type : quadrature_points_coordinates.elementTypes( - Model::spatial_dimension, ghost_type)) { - this->getFEEngine().computeIntegrationPointsCoordinates( - quadrature_points_coordinates(type, ghost_type), type, ghost_type); - } - - mat_non_local->initMaterialNonLocal(); - - mat_non_local->insertIntegrationPointsInNeighborhoods( - ghost_type, quadrature_points_coordinates); - } -} - /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses(GhostType ghost_type) { - for (auto & mat : materials) { - if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) { + for_each_constitutive_law([&](auto && material) { + if (not aka::is_of_type<MaterialNonLocalInterface>(mat)) { continue; } - auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat); + auto & mat_non_local = aka::as_type<MaterialNonLocalInterface>(mat); mat_non_local.computeNonLocalStresses(ghost_type); - } -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::updateLocalInternal( - ElementTypeMapReal & internal_flat, GhostType ghost_type, - ElementKind kind) { - const ID field_name = internal_flat.getName(); - for (auto & material : materials) { - if (material->isInternal<Real>(field_name, kind)) { - material->flattenInternal(field_name, internal_flat, ghost_type, kind); - } - } -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::updateNonLocalInternal( - ElementTypeMapReal & internal_flat, GhostType ghost_type, - ElementKind kind) { - - const ID field_name = internal_flat.getName(); - - for (auto & mat : materials) { - if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) { - continue; - } - - auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat); - mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, - kind); - } + }); } /* -------------------------------------------------------------------------- */ FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { return getFEEngineClassBoundary<MyFEEngineType>(name); } -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::splitElementByMaterial( - const Array<Element> & elements, - std::vector<Array<Element>> & elements_per_mat) const { - for (const auto & el : elements) { - Element mat_el = el; - mat_el.element = this->material_local_numbering(el); - elements_per_mat[this->material_index(el)].push_back(mat_el); - } -} - /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::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); } switch (tag) { - case SynchronizationTag::_material_id: { - size += elements.size() * sizeof(UInt); - break; - } case SynchronizationTag::_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * Model::spatial_dimension; // mass vector break; } case SynchronizationTag::_smm_for_gradu: { size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real); // displacement break; } case SynchronizationTag::_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * Model::spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case SynchronizationTag::_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5; break; } default: { } } - if (tag != SynchronizationTag::_material_id) { - splitByMaterial(elements, [&](auto && mat, auto && elements) { - size += mat.getNbData(elements, tag); - }); - } + ConstitutiveLawsHandler<Material>::getNbData(elements, tag); AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array<Element> & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { - case SynchronizationTag::_material_id: { - packElementalDataHelper(material_index, buffer, elements, false, - getFEEngine()); - break; - } case SynchronizationTag::_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { packNodalDataHelper(*displacement, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*acceleration, buffer, elements, mesh); packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } - if (tag != SynchronizationTag::_material_id) { - splitByMaterial(elements, [&](auto && mat, auto && elements) { - mat.packData(buffer, elements, tag); - }); - } + ConstitutiveLawsHandler<Material>::packData(buffer, elements, tag); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array<Element> & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { - case SynchronizationTag::_material_id: { - for (auto && element : elements) { - UInt recv_mat_index; - buffer >> recv_mat_index; - UInt & mat_index = material_index(element); - if (mat_index != UInt(-1)) { - continue; - } - - // add ghosts element to the correct material - mat_index = recv_mat_index; - UInt index = materials[mat_index]->addElement(element); - material_local_numbering(element) = index; - } - break; - } case SynchronizationTag::_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*acceleration, buffer, elements, mesh); unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } - if (tag != SynchronizationTag::_material_id) { - splitByMaterial(elements, [&](auto && mat, auto && elements) { - mat.unpackData(buffer, elements, tag); - }); - } - + ConstitutiveLawsHandler<Material>::unpackData(buffer, elements, tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch (tag) { case SynchronizationTag::_smm_uv: { size += sizeof(Real) * Model::spatial_dimension * 2; break; } case SynchronizationTag::_smm_res: /* FALLTHRU */ case SynchronizationTag::_smm_mass: { size += sizeof(Real) * Model::spatial_dimension; break; } case SynchronizationTag::_for_dump: { size += sizeof(Real) * Model::spatial_dimension * 5; break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array<UInt> & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { packDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { packDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); packDOFDataHelper(*acceleration, buffer, dofs); packDOFDataHelper(*internal_force, buffer, dofs); packDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { unpackDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { unpackDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); unpackDOFDataHelper(*acceleration, buffer, dofs); unpackDOFDataHelper(*internal_force, buffer, dofs); unpackDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::applyEigenGradU( + const Matrix<Real> & prescribed_eigen_grad_u, const ID & material_name, + const GhostType ghost_type) { + AKANTU_DEBUG_ASSERT(prescribed_eigen_grad_u.size() == + spatial_dimension * spatial_dimension, + "The prescribed grad_u is not of the good size"); + for (auto & material : this->getConstitutiveLaws()) { + if (material->getName() == material_name) { + material->applyEigenGradU(prescribed_eigen_grad_u, ghost_type); + } + } +} + +/* -------------------------------------------------------------------------- */ + } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index c59136979..290fee36b 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,601 +1,497 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Apr 09 2021 * * @brief Model of Solid Mechanics * * * @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 "boundary_condition.hh" +#include "constitutive_laws_handler.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" -#include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; -class NonLocalManager; template <ElementKind kind, class IntegrationOrderFunctor> class IntegratorGauss; template <ElementKind kind> class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor<Element>, public DataAccessor<UInt>, public BoundaryCondition<SolidMechanicsModel>, - public NonLocalManagerCallback, + public ConstitutiveLawsHandler<Material>, public EventHandlerManager<SolidMechanicsModelEventHandler> { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - class NewMaterialElementsEvent : public NewElementsEvent { - public: - AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array<UInt> &); - AKANTU_GET_MACRO(MaterialList, material, const Array<UInt> &); - - protected: - Array<UInt> material; - }; - using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>; protected: using EventManager = EventHandlerManager<SolidMechanicsModelEventHandler>; public: SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", std::shared_ptr<DOFManager> dof_manager = nullptr, ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple<ID, TimeStepSolverType> getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(bool need_to_reassemble = false); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() const override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) const override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: - /// register an empty material of a given type - Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, - const ID & opt_param); - - /// reassigns materials depending on the material selector - virtual void reassignMaterial(); - /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u, const ID & material_name, GhostType ghost_type = _not_ghost); -protected: - /// register a material in the dynamic database - Material & registerNewMaterial(const ParserSection & mat_section); - - /// read the material files to instantiate all the materials - void instantiateMaterials(); - - /// set the element_id_by_material and add the elements to the good materials - virtual void - assignMaterialToElements(const ElementTypeMapArray<UInt> * filter = nullptr); - /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); public: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); protected: /// fill a vector of rho void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); - /* ------------------------------------------------------------------------ */ - /* NonLocalManager inherited members */ - /* ------------------------------------------------------------------------ */ -protected: - void initializeNonLocal() override; - - void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; - - void computeNonLocalStresses(GhostType ghost_type) override; - - void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; - - /// update the values of the non local internal - void updateLocalInternal(ElementTypeMapReal & internal_flat, - GhostType ghost_type, ElementKind kind) override; - - /// copy the results of the averaging in the materials - void updateNonLocalInternal(ElementTypeMapReal & internal_flat, - GhostType ghost_type, ElementKind kind) override; - /* ------------------------------------------------------------------------ */ /* 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: - void - splitElementByMaterial(const Array<Element> & elements, - std::vector<Array<Element>> & elements_per_mat) const; - - template <typename Operation> - void splitByMaterial(const Array<Element> & elements, Operation && op) const; - /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array<UInt> & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array<UInt> & element_list, const Array<UInt> & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array<Element> & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array<Element> & element_list, const ElementTypeMapArray<UInt> & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array<Element> & /*unused*/, const Array<Element> & /*unused*/, const ElementTypeMapArray<UInt> & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); - //! decide wether a field is a material internal or not - bool isInternal(const std::string & field_name, ElementKind element_kind); - //! give the amount of data per element - virtual ElementTypeMap<UInt> - getInternalDataPerElem(const std::string & field_name, ElementKind kind); - - //! flatten a given material internal field - ElementTypeMapArray<Real> & - flattenInternal(const std::string & field_name, ElementKind kind, - GhostType ghost_type = _not_ghost); - //! flatten all the registered material internals - void flattenAllRegisteredInternals(ElementKind kind); - - //! inverse operation of the flatten - void inflateInternal(const std::string & field_name, - const ElementTypeMapArray<Real> & field, - ElementKind kind, GhostType ghost_type = _not_ghost); - 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> createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr<dumpers::Field> createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; void dump(const std::string & dumper_name) override; void dump(const std::string & dumper_name, UInt step) override; void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array<Real> & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array<Real> & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials inline decltype(auto) getMaterials(); /// get an iterable on the materials inline decltype(auto) getMaterials() const; /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by numerical material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material (by material name) inline const Material & getMaterial(const Element & element) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials - inline UInt getNbMaterials() const { return materials.size(); } - - /// give the material internal index from its id - Int getInternalIndexFromID(const ID & id) const; + inline UInt getNbMaterials() const { return getNbConstitutiveLaws(); } /// compute the stable time step Real getStableTimeStep(); /** * @brief Returns the total energy for a given energy type * * Energy types of SolidMechanicsModel expected as argument are: * - `kinetic` * - `external work` * * Other energy types are passed on to the materials. All materials should * define a `potential` energy type. For additional energy types, see material * documentation. */ Real getEnergy(const std::string & energy_id); /// Compute energy for an element type and material index Real getEnergy(const std::string & energy_id, ElementType type, UInt index); /// Compute energy for an individual element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } /// Compute energy for an element group Real getEnergy(const ID & energy_id, const ID & group_id); - AKANTU_GET_MACRO(MaterialByElement, material_index, - const ElementTypeMapArray<UInt> &); - AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, - const ElementTypeMapArray<UInt> &); - - /// vectors containing local material element index for each global element - /// index - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, - UInt); - // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, - material_local_numbering, UInt); - // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, - // material_local_numbering, UInt); - - AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, material_selector, - std::shared_ptr<MaterialSelector>); + // this function is kept for backward compatinility + decltype(auto) getMaterialByElement() const { + return this->getConstitutiveLawByElement(); + } + + // this function is kept for backward compatinility + decltype(auto) getMaterialLocalNumbering() const { + return this->getConstitutiveLawLocalNumbering(); + } + + // this function is kept for backward compatinility + decltype(auto) getMaterialByElement(ElementType type, + GhostType ghost_type = _not_ghost) const { + return this->getConstitutiveLawByElement(type, ghost_type); + } + + // this function is kept for backward compatinility + decltype(auto) + getMaterialLocalNumbering(ElementType type, + GhostType ghost_type = _not_ghost) const { + return this->getConstitutiveLawLocalNumbering(type, ghost_type); + } + + // this function is kept for backward compatinility + decltype(auto) getMaterialSelector() { + return this->getConstitutiveLawSelector(); + } + + // this function is kept for backward compatinility void setMaterialSelector(std::shared_ptr<MaterialSelector> material_selector) { - this->material_selector = std::move(material_selector); + this->setConstitutiveLawSelector(std::move(material_selector)); } - /// Access the non_local_manager interface - AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); - /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array UInt displacement_release{0}; /// release version of the current_position array UInt current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; - /// mapping between material name and material internal id - std::map<std::string, UInt> materials_names_to_id; - protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr<Array<Real>> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr<Array<Real>> previous_displacement; /// increment of displacement std::unique_ptr<Array<Real>> displacement_increment; /// lumped mass array std::unique_ptr<Array<Real>> mass; /// velocities array std::unique_ptr<Array<Real>> velocity; /// accelerations array std::unique_ptr<Array<Real>> acceleration; /// external forces array std::unique_ptr<Array<Real>> external_force; /// internal forces array std::unique_ptr<Array<Real>> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr<Array<bool>> blocked_dofs; /// array of current position used during update residual std::unique_ptr<Array<Real>> current_position; - - /// Arrays containing the material index for each element - ElementTypeMapArray<UInt> material_index; - - /// Arrays containing the position in the element filter of the material - /// (material's local numbering) - ElementTypeMapArray<UInt> material_local_numbering; - - /// list of used materials - std::vector<std::unique_ptr<Material>> materials; - - /// class defining of to choose a material - std::shared_ptr<MaterialSelector> material_selector; - - using flatten_internal_map = - std::map<std::pair<std::string, ElementKind>, - std::unique_ptr<ElementTypeMapArray<Real>>>; - - /// tells if the material are instantiated - flatten_internal_map registered_internals; - - /// non local manager - std::unique_ptr<NonLocalManager> non_local_manager; - - /// tells if the material are instantiated - bool are_materials_instantiated{false}; - - friend class Material; - - template <class Model_> friend class CouplerSolidContactTemplate; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" - -#include "solid_mechanics_model_inline_impl.hh" -#include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.hh deleted file mode 100644 index d3dfed74f..000000000 --- a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.hh +++ /dev/null @@ -1,118 +0,0 @@ -/** - * @file solid_mechanics_model_inline_impl.hh - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch> - * @author Nicolas Richart <nicolas.richart@epfl.ch> - * - * @date creation: Wed Aug 04 2010 - * @date last modification: Fri Mar 26 2021 - * - * @brief Implementation of the inline functions of the SolidMechanicsModel - * class - * - * - * @section LICENSE - * - * Copyright (©) 2015-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_named_argument.hh" -#include "material_selector.hh" -#include "material_selector_tmpl.hh" -#include "solid_mechanics_model.hh" -/* -------------------------------------------------------------------------- */ - -#ifndef AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_HH_ -#define AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_HH_ - -namespace akantu { - -/* -------------------------------------------------------------------------- */ -inline decltype(auto) SolidMechanicsModel::getMaterials() { - return make_dereference_adaptor(materials); -} - -/* -------------------------------------------------------------------------- */ -inline decltype(auto) SolidMechanicsModel::getMaterials() const { - return make_dereference_adaptor(materials); -} - -/* -------------------------------------------------------------------------- */ -inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) { - AKANTU_DEBUG_ASSERT(mat_index < materials.size(), - "The model " << id << " has no material no " - << mat_index); - return *materials.at(mat_index); -} - -/* -------------------------------------------------------------------------- */ -inline const Material & SolidMechanicsModel::getMaterial(UInt mat_index) const { - AKANTU_DEBUG_ASSERT(mat_index < materials.size(), - "The model " << id << " has no material no " - << mat_index); - return *materials.at(mat_index); -} - -/* -------------------------------------------------------------------------- */ -inline Material & SolidMechanicsModel::getMaterial(const std::string & name) { - std::map<std::string, UInt>::const_iterator it = - materials_names_to_id.find(name); - if (it == materials_names_to_id.end()) { - AKANTU_SILENT_EXCEPTION("The model " << id << " has no material named " - << name); - } - - return *materials[it->second]; -} - -/* -------------------------------------------------------------------------- */ -inline const Material & -SolidMechanicsModel::getMaterial(const Element & element) const { - auto mat_id = material_index(element); - return *materials[mat_id]; -} - -/* -------------------------------------------------------------------------- */ -inline UInt -SolidMechanicsModel::getMaterialIndex(const std::string & name) const { - auto it = materials_names_to_id.find(name); - if (it == materials_names_to_id.end()) { - AKANTU_SILENT_EXCEPTION("The model " << id << " has no material named " - << name); - } - - return it->second; -} - -/* -------------------------------------------------------------------------- */ -inline const Material & -SolidMechanicsModel::getMaterial(const std::string & name) const { - auto it = materials_names_to_id.find(name); - if (it == materials_names_to_id.end()) { - AKANTU_SILENT_EXCEPTION("The model " << id << " has no material named " - << name); - } - return *materials[it->second]; -} - -/* -------------------------------------------------------------------------- */ -} // namespace akantu - -#endif /* AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_material.cc b/src/model/solid_mechanics/solid_mechanics_model_material.cc deleted file mode 100644 index c06fea251..000000000 --- a/src/model/solid_mechanics/solid_mechanics_model_material.cc +++ /dev/null @@ -1,247 +0,0 @@ -/** - * @file solid_mechanics_model_material.cc - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> - * @author Nicolas Richart <nicolas.richart@epfl.ch> - * - * @date creation: Fri Nov 26 2010 - * @date last modification: Fri Mar 26 2021 - * - * @brief instatiation of 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_math.hh" -#include "material_non_local.hh" -#include "mesh_iterators.hh" -#include "non_local_manager.hh" -#include "solid_mechanics_model.hh" -/* -------------------------------------------------------------------------- */ - -namespace akantu { - -/* -------------------------------------------------------------------------- */ -Material & -SolidMechanicsModel::registerNewMaterial(const ParserSection & section) { - std::string mat_name; - std::string mat_type = section.getName(); - std::string opt_param = section.getOption(); - - try { - std::string tmp = section.getParameter("name"); - mat_name = tmp; /** this can seam 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 material of type \'" - << mat_type - << "\' in the input file has been defined without a name!"); - } - Material & mat = this->registerNewMaterial(mat_name, mat_type, opt_param); - - mat.parseSection(section); - - return mat; -} - -/* -------------------------------------------------------------------------- */ -Material & SolidMechanicsModel::registerNewMaterial(const ID & mat_name, - const ID & mat_type, - const ID & opt_param) { - AKANTU_DEBUG_ASSERT(materials_names_to_id.find(mat_name) == - materials_names_to_id.end(), - "A material with this name '" - << mat_name << "' has already been registered. " - << "Please use unique names for materials"); - - UInt mat_count = materials.size(); - materials_names_to_id[mat_name] = mat_count; - - ID mat_id = this->id + ":" + std::to_string(mat_count) + ":" + mat_type; - - std::unique_ptr<Material> material = MaterialFactory::getInstance().allocate( - mat_type, spatial_dimension, opt_param, *this, mat_id); - material->setName(mat_name); - materials.push_back(std::move(material)); - - return *(materials.back()); -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::instantiateMaterials() { - ParserSection model_section; - bool is_empty; - std::tie(model_section, is_empty) = this->getParserSection(); - - if (not is_empty) { - auto model_materials = model_section.getSubSections(ParserType::_material); - for (const auto & section : model_materials) { - this->registerNewMaterial(section); - } - } - - auto sub_sections = this->parser.getSubSections(ParserType::_material); - for (const auto & section : sub_sections) { - this->registerNewMaterial(section); - } - -#ifdef AKANTU_DAMAGE_NON_LOCAL - for (auto & material : materials) { - if (dynamic_cast<MaterialNonLocalInterface *>(material.get()) == nullptr) { - continue; - } - - this->non_local_manager = std::make_unique<NonLocalManager>( - *this, *this, id + ":non_local_manager"); - break; - } -#endif - - if (materials.empty()) { - AKANTU_EXCEPTION("No materials where instantiated for the model" - << getID()); - } - are_materials_instantiated = true; -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::assignMaterialToElements( - const ElementTypeMapArray<UInt> * filter) { - - for_each_element( - mesh, - [&](auto && element) { - UInt mat_index = (*material_selector)(element); - AKANTU_DEBUG_ASSERT( - mat_index < materials.size(), - "The material selector returned an index that does not exists"); - material_index(element) = mat_index; - }, - _element_filter = filter, _ghost_type = _not_ghost); - - if (non_local_manager) { - non_local_manager->synchronize(*this, SynchronizationTag::_material_id); - } - - for_each_element( - mesh, - [&](auto && element) { - auto mat_index = material_index(element); - auto index = materials[mat_index]->addElement(element); - material_local_numbering(element) = index; - }, - _element_filter = filter, _ghost_type = _not_ghost); - - // synchronize the element material arrays - this->synchronize(SynchronizationTag::_material_id); -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::initMaterials() { - AKANTU_DEBUG_ASSERT(not materials.empty(), "No material to initialize !"); - - // if (!are_materials_instantiated) - // instantiateMaterials(); - - this->assignMaterialToElements(); - - for (auto & material : materials) { - /// init internals properties - material->initMaterial(); - } - - this->synchronize(SynchronizationTag::_smm_init_mat); - - if (this->non_local_manager) { - this->non_local_manager->initialize(); - } -} - -/* -------------------------------------------------------------------------- */ -Int SolidMechanicsModel::getInternalIndexFromID(const ID & id) const { - AKANTU_DEBUG_IN(); - - auto it = materials.begin(); - auto end = materials.end(); - - for (; it != end; ++it) { - if ((*it)->getID() == id) { - AKANTU_DEBUG_OUT(); - return (it - materials.begin()); - } - } - - AKANTU_DEBUG_OUT(); - return -1; -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::reassignMaterial() { - AKANTU_DEBUG_IN(); - - std::vector<Array<Element>> element_to_add(materials.size()); - std::vector<Array<Element>> element_to_remove(materials.size()); - - for_each_element( - mesh, - [&](auto && element) { - auto old_material = material_index(element); - auto new_material = (*material_selector)(element); - if (old_material != new_material) { - element_to_add[new_material].push_back(element); - element_to_remove[old_material].push_back(element); - } - }, - _spatial_dimension = spatial_dimension); - - for (auto && data : enumerate(materials)) { - auto mat_index = std::get<0>(data); - auto & mat = *std::get<1>(data); - - mat.removeElements(element_to_remove[mat_index]); - mat.addElements(element_to_add[mat_index]); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::applyEigenGradU( - const Matrix<Real> & prescribed_eigen_grad_u, const ID & material_name, - const GhostType ghost_type) { - AKANTU_DEBUG_ASSERT(prescribed_eigen_grad_u.size() == - spatial_dimension * spatial_dimension, - "The prescribed grad_u is not of the good size"); - for (auto & material : materials) { - if (material->getName() == material_name) { - material->applyEigenGradU(prescribed_eigen_grad_u, ghost_type); - } - } -} - -/* -------------------------------------------------------------------------- */ - -} // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh b/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh deleted file mode 100644 index 441ed43fd..000000000 --- a/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh +++ /dev/null @@ -1,63 +0,0 @@ -/** - * @file solid_mechanics_model_tmpl.hh - * - * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> - * @author Dana Christen <dana.christen@epfl.ch> - * @author Nicolas Richart <nicolas.richart@epfl.ch> - * - * @date creation: Fri Jun 18 2010 - * @date last modification: Fri Mar 26 2021 - * - * @brief template part of solid mechanics model - * - * - * @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 "material.hh" -#include "solid_mechanics_model.hh" -/* -------------------------------------------------------------------------- */ - -#ifndef AKANTU_SOLID_MECHANICS_MODEL_TMPL_HH_ -#define AKANTU_SOLID_MECHANICS_MODEL_TMPL_HH_ - -namespace akantu { - -#define FWD(...) ::std::forward<decltype(__VA_ARGS__)>(__VA_ARGS__) - -/* -------------------------------------------------------------------------- */ -template <typename Operation> -void SolidMechanicsModel::splitByMaterial(const Array<Element> & elements, - Operation && op) const { - std::vector<Array<Element>> elements_per_mat(materials.size()); - this->splitElementByMaterial(elements, elements_per_mat); - - for (auto && mat : zip(materials, elements_per_mat)) { - FWD(op)(FWD(*std::get<0>(mat)), FWD(std::get<1>(mat))); - } -} - -#undef FWD -/* -------------------------------------------------------------------------- */ - -} // namespace akantu - -#endif /* AKANTU_SOLID_MECHANICS_MODEL_TMPL_HH_ */