diff --git a/packages/core.cmake b/packages/core.cmake index 146fa2d57..66bfebd83 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,438 +1,439 @@ #=============================================================================== # @file core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @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 . # #=============================================================================== 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 cxx_std_17 FEATURES_PRIVATE cxx_lambdas cxx_nullptr cxx_range_for cxx_delegating_constructors DEPENDS INTERFACE akantu_iterators Boost Eigen3 ) 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_types_eigen_matrix_plugin.hh common/aka_types_eigen_matrix_base_plugin.hh common/aka_types_eigen_plain_object_base_plugin.hh common/aka_view_iterators.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/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.cc 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/damage_non_local.cmake b/packages/damage_non_local.cmake index 3de0a77d4..0b9b38300 100644 --- a/packages/damage_non_local.cmake +++ b/packages/damage_non_local.cmake @@ -1,62 +1,64 @@ #=============================================================================== # @file damage_non_local.cmake # # @author Aurelia Isabel Cuba Ramos # @author Nicolas Richart # # @date creation: Fri Jun 15 2012 # @date last modification: Thu Dec 17 2020 # # @brief package description for non-local 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 . # #=============================================================================== package_declare(damage_non_local DESCRIPTION "Package for Non-local damage constitutives laws Akantu" DEPENDS lapack) package_declare_sources(damage_non_local model/solid_mechanics/materials/material_damage/material_damage_non_local.hh model/solid_mechanics/materials/material_damage/material_marigo_non_local.cc model/solid_mechanics/materials/material_damage/material_marigo_non_local.hh model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh model/solid_mechanics/materials/material_damage/material_mazars_non_local_tmpl.hh model/solid_mechanics/materials/material_damage/material_von_mises_mazars_non_local.cc model/solid_mechanics/materials/material_damage/material_von_mises_mazars_non_local.hh model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh + model/solid_mechanics/materials/weight_functions/damaged_weight_function.cc model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.hh model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh + model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.cc model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.hh model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.hh model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh model/solid_mechanics/materials/weight_functions/stress_based_weight_function.cc model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.hh ) package_declare_material_infos(damage_non_local LIST AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST INCLUDE material_non_local_includes.hh ) diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh index ef7b989fc..4e593de7d 100644 --- a/src/common/aka_array_tmpl.hh +++ b/src/common/aka_array_tmpl.hh @@ -1,1167 +1,1167 @@ /** * @file aka_array_tmpl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Thu Jul 15 2010 * @date last modification: Fri Feb 26 2021 * * @brief Inline functions of the classes Array and ArrayBase * * * @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 . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Array */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" // NOLINT /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ //#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__ //#define __AKANTU_AKA_ARRAY_TMPL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { namespace debug { struct ArrayException : public Exception {}; } // namespace debug /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(Int size, Int nb_component, const ID & id) : ArrayBase(id) { allocate(size, nb_component); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(Int size, Int nb_component, const_reference value, const ID & id) : ArrayBase(id) { allocate(size, nb_component, value); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(const ArrayDataLayer & vect, const ID & id) : ArrayBase(vect, id) { this->data_storage = vect.data_storage; this->size_ = vect.size_; this->nb_component = vect.nb_component; this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer( const std::vector & vect) { this->data_storage = vect; this->size_ = vect.size(); this->nb_component = 1; this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer & ArrayDataLayer::operator=(const ArrayDataLayer & other) { if (this != &other) { this->data_storage = other.data_storage; this->nb_component = other.nb_component; this->size_ = other.size_; this->values = this->data_storage.data(); } return *this; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::allocate(Int new_size, Int nb_component) { this->nb_component = nb_component; this->resize(new_size); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::allocate(Int new_size, Int nb_component, const T & val) { this->nb_component = nb_component; this->resize(new_size, val); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(Int new_size) { this->data_storage.resize(new_size * this->nb_component); this->values = this->data_storage.data(); this->size_ = new_size; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(Int new_size, const T & value) { this->data_storage.resize(new_size * this->nb_component, value); this->values = this->data_storage.data(); this->size_ = new_size; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::reserve(Int size, Int new_size) { if (new_size != -1) { this->data_storage.resize(new_size * this->nb_component); } this->data_storage.reserve(size * this->nb_component); this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ /** * append a tuple to the array with the value value for all components * @param value the new last tuple or the array will contain nb_component copies * of value */ template inline void ArrayDataLayer::push_back(const T & value) { this->data_storage.push_back(value); this->values = this->data_storage.data(); this->size_ += 1; } /* -------------------------------------------------------------------------- */ /** * append a matrix or a vector to the array * @param new_elem a reference to a Matrix or Vector */ template template inline void ArrayDataLayer::push_back( const Eigen::MatrixBase & new_elem) { AKANTU_DEBUG_ASSERT( nb_component == new_elem.size(), "The vector(" << new_elem.size() << ") as not a size compatible with the Array (nb_component=" << nb_component << ")."); for (Idx i = 0; i < new_elem.size(); ++i) { this->data_storage.push_back(new_elem.array()[i]); } this->values = this->data_storage.data(); this->size_ += 1; } /* -------------------------------------------------------------------------- */ template inline Int ArrayDataLayer::getAllocatedSize() const { return this->data_storage.capacity() / this->nb_component; } /* -------------------------------------------------------------------------- */ template inline Int ArrayDataLayer::getMemorySize() const { return this->data_storage.capacity() * sizeof(T); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class ArrayDataLayer : public ArrayBase { public: using value_type = T; using reference = value_type &; using pointer_type = value_type *; using size_type = typename ArrayBase::size_type; using const_reference = const value_type &; public: ~ArrayDataLayer() override { deallocate(); } /// Allocation of a new vector ArrayDataLayer(Int size = 0, Int nb_component = 1, const ID & id = "") : ArrayBase(id) { allocate(size, nb_component); } /// Allocation of a new vector with a default value ArrayDataLayer(Int size, Int nb_component, const_reference value, const ID & id = "") : ArrayBase(id) { allocate(size, nb_component, value); } /// Copy constructor (deep copy) ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = "") : ArrayBase(vect, id) { allocate(vect.size(), vect.getNbComponent()); std::copy_n(vect.data(), this->size_ * this->nb_component, values); } /// Copy constructor (deep copy) explicit ArrayDataLayer(const std::vector & vect) { allocate(vect.size(), 1); std::copy_n(vect.data(), this->size_ * this->nb_component, values); } // copy operator inline ArrayDataLayer & operator=(const ArrayDataLayer & other) { if (this != &other) { allocate(other.size(), other.getNbComponent()); std::copy_n(other.data(), this->size_ * this->nb_component, values); } return *this; } // move constructor inline ArrayDataLayer(ArrayDataLayer && other) noexcept = default; // move assign inline ArrayDataLayer & operator=(ArrayDataLayer && other) noexcept = default; protected: // deallocate the memory virtual void deallocate() { // NOLINTNEXTLINE(cppcoreguidelines-owning-memory, // cppcoreguidelines-no-malloc) free(this->values); } // allocate the memory virtual inline void allocate(Int size, Int nb_component) { if (size != 0) { // malloc can return a non NULL pointer in case size is 0 this->values = static_cast( // NOLINT std::malloc(nb_component * size * sizeof(T))); // NOLINT } if (this->values == nullptr and size != 0) { throw std::bad_alloc(); } this->nb_component = nb_component; this->allocated_size = this->size_ = size; } // allocate and initialize the memory virtual inline void allocate(Int size, Int nb_component, const T & value) { allocate(size, nb_component); std::fill_n(values, size * nb_component, value); } public: /// append a tuple of size nb_component containing value inline void push_back(const_reference value) { resize(this->size_ + 1, value); } /// append a Vector or a Matrix template inline void push_back(const Eigen::MatrixBase & new_elem) { AKANTU_DEBUG_ASSERT( nb_component == new_elem.size(), "The vector(" << new_elem.size() << ") as not a size compatible with the Array (nb_component=" << nb_component << ")."); this->resize(this->size_ + 1); - make_view(*this, nb_component).begin()[this->size_].array() = + make_view(*this, nb_component).begin()[this->size_ - 1].array() = new_elem.array(); } /// changes the allocated size but not the size virtual void reserve(Int size, Int new_size = Int(-1)) { auto tmp_size = this->size_; if (new_size != Int(-1)) { tmp_size = new_size; } this->resize(size); this->size_ = std::min(this->size_, tmp_size); } /// change the size of the Array virtual void resize(Int size) { if (size * this->nb_component == 0) { free(values); // NOLINT: cppcoreguidelines-no-malloc values = nullptr; this->allocated_size = 0; } else { if (this->values == nullptr) { this->allocate(size, this->nb_component); return; } Int diff = size - allocated_size; Int size_to_allocate = (std::abs(diff) > AKANTU_MIN_ALLOCATION) ? size : (diff > 0) ? allocated_size + AKANTU_MIN_ALLOCATION : allocated_size; if (size_to_allocate == allocated_size) { // otherwhy the reserve + push_back might fail... this->size_ = size; return; } auto * tmp_ptr = reinterpret_cast( // NOLINT realloc(this->values, size_to_allocate * this->nb_component * sizeof(T))); if (tmp_ptr == nullptr) { throw std::bad_alloc(); } this->values = tmp_ptr; this->allocated_size = size_to_allocate; } this->size_ = size; } /// change the size of the Array and initialize the values virtual void resize(Int size, const T & val) { Int tmp_size = this->size_; this->resize(size); if (size > tmp_size) { // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) std::fill_n(values + this->nb_component * tmp_size, (size - tmp_size) * this->nb_component, val); } } /// get the amount of space allocated in bytes inline size_type getMemorySize() const final { return this->allocated_size * this->nb_component * sizeof(T); } /// Get the real size allocated in memory inline Int getAllocatedSize() const { return this->allocated_size; } /// give the address of the memory allocated for this vector [[deprecated("use data instead to be stl compatible")]] T * storage() const { return values; }; T * data() const { return values; }; protected: /// allocation type agnostic data access T * values{nullptr}; Int allocated_size{0}; }; /* -------------------------------------------------------------------------- */ template inline auto Array::operator()(Int i, Int j) -> reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component), "The value at position [" << i << "," << j << "] is out of range in array \"" << this->id << "\""); return this->values[i * this->nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline auto Array::operator()(Int i, Int j) const -> const_reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component), "The value at position [" << i << "," << j << "] is out of range in array \"" << this->id << "\""); // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) return this->values[i * this->nb_component + j]; } template inline auto Array::operator[](Int i) -> reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component), "The value at position [" << i << "] is out of range in array \"" << this->id << "\""); return this->values[i]; } /* -------------------------------------------------------------------------- */ template inline auto Array::operator[](Int i) const -> const_reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component), "The value at position [" << i << "] is out of range in array \"" << this->id << "\""); return this->values[i]; } /* -------------------------------------------------------------------------- */ /** * erase an element. If the erased element is not the last of the array, the * last element is moved into the hole in order to maintain contiguity. This * may invalidate existing iterators (For instance an iterator obtained by * Array::end() is no longer correct) and will change the order of the * elements. * @param i index of element to erase */ template inline void Array::erase(Idx i) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT((this->size_ > 0), "The array is empty"); AKANTU_DEBUG_ASSERT((i < this->size_), "The element at position [" << i << "] is out of range (" << i << ">=" << this->size_ << ")"); if (i != (this->size_ - 1)) { for (Idx j = 0; j < this->nb_component; ++j) { this->values[i * this->nb_component + j] = this->values[(this->size_ - 1) * this->nb_component + j]; } } this->resize(this->size_ - 1); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ template template inline auto Array::erase(const view_iterator & it) { auto && curr = it.data(); Idx pos = (curr - this->values) / this->nb_component; erase(pos); view_iterator rit = it; return --rit; } /* -------------------------------------------------------------------------- */ /** * Subtract another array entry by entry from this array in place. Both arrays * must * have the same size and nb_component. If the arrays have different shapes, * code compiled in debug mode will throw an expeption and optimised code * will behave in an unpredicted manner * @param other array to subtract from this * @return reference to modified this */ template Array & Array::operator-=(const Array & vect) { AKANTU_DEBUG_ASSERT((this->size_ == vect.size_) && (this->nb_component == vect.nb_component), "The too array don't have the same sizes"); T * a = this->values; T * b = vect.data(); for (Idx i = 0; i < this->size_ * this->nb_component; ++i) { *a -= *b; ++a; ++b; } return *this; } /* -------------------------------------------------------------------------- */ /** * Add another array entry by entry to this array in * place. Both arrays must have the same size and * nb_component. If the arrays have different shapes, code * compiled in debug mode will throw an expeption and * optimised code will behave in an unpredicted manner * @param other array to add to this * @return reference to modified this */ template Array & Array::operator+=(const Array & vect) { AKANTU_DEBUG_ASSERT((this->size_ == vect.size()) && (this->nb_component == vect.nb_component), "The too array don't have the same sizes"); T * a = this->values; T * b = vect.data(); for (Idx i = 0; i < this->size_ * this->nb_component; ++i) { *a++ += *b++; } return *this; } /* -------------------------------------------------------------------------- */ /** * Multiply all entries of this array by a scalar in place * @param alpha scalar multiplicant * @return reference to modified this */ template auto Array::operator*=(const T & alpha) -> Array & { T * a = this->values; for (Idx i = 0; i < this->size_ * this->nb_component; ++i) { *a++ *= alpha; } return *this; } /* ------------------------------------------------------------------------- */ /** * Compare this array element by element to another. * @param other array to compare to * @return true it all element are equal and arrays have * the same shape, else false */ template bool Array::operator==(const Array & other) const { bool equal = this->nb_component == other.nb_component && this->size_ == other.size_ && this->id == other.id; if (not equal) { return false; } if (this->values == other.data()) { return true; } return std::equal(this->values, this->values + this->size_ * this->nb_component, other.data()); } /* ------------------------------------------------------------------------ */ template bool Array::operator!=(const Array & other) const { return !operator==(other); } /* ------------------------------------------------------------------------ */ /** * set all tuples of the array to a given vector or matrix * @param vm Matrix or Vector to fill the array with */ template template ::value> *> inline void Array::set(const C & elem) { AKANTU_DEBUG_ASSERT( this->nb_component == elem.array().size(), "The size of the object does not match the number of components"); for (auto && v : make_view(*this, this->nb_component)) { v = elem.array(); } } /* ------------------------------------------------------------------------ */ template void Array::append(const Array & other) { AKANTU_DEBUG_ASSERT( this->nb_component == other.nb_component, "Cannot append an array with a different number of component"); Idx old_size = this->size_; this->resize(this->size_ + other.size()); T * tmp = this->values + this->nb_component * old_size; std::copy_n(other.data(), other.size() * this->nb_component, tmp); } /* ------------------------------------------------------------------------ */ /* Functions Array */ /* ------------------------------------------------------------------------ */ template Array::Array(Int size, Int nb_component, const ID & id) : parent(size, nb_component, id) {} template <> inline Array::Array(Int size, Int nb_component, const ID & id) : parent(size, nb_component, "", id) {} /* ------------------------------------------------------------------------ */ template Array::Array(Int size, Int nb_component, const_reference value, const ID & id) : parent(size, nb_component, value, id) {} /* ------------------------------------------------------------------------ */ template Array::Array(const Array & vect, const ID & id) : parent(vect, id) {} /* ------------------------------------------------------------------------ */ template auto Array::operator=(const Array & other) -> Array & { AKANTU_DEBUG_WARNING("You are copying the array " << this->id << " are you sure it is on purpose"); if (&other == this) { return *this; } parent::operator=(other); return *this; } /* ------------------------------------------------------------------------ */ template Array::Array(const std::vector & vect) : parent(vect) {} /* ------------------------------------------------------------------------ */ template Array::~Array() = default; /* ------------------------------------------------------------------------ */ /** * search elem in the array, return the position of the * first occurrence or -1 if not found * @param elem the element to look for * @return index of the first occurrence of elem or -1 if * elem is not present */ template Idx Array::find(const_reference elem) const { AKANTU_DEBUG_IN(); auto begin = this->begin(); auto end = this->end(); auto it = std::find(begin, end, elem); AKANTU_DEBUG_OUT(); return (it != end) ? it - begin : Idx(-1); } /* ------------------------------------------------------------------------ */ template template ::value> *> inline Idx Array::find(const V & elem) { AKANTU_DEBUG_ASSERT(elem.size() == this->nb_component, "Cannot find an element with a wrong size (" << elem.size() << ") != " << this->nb_component); auto && view = make_view(*this, elem.size()); auto begin = view.begin(); auto end = view.end(); auto it = std::find(begin, end, elem); return (it != end) ? it - begin : Idx(-1); } /* ------------------------------------------------------------------------ */ /** * copy the content of another array. This overwrites the * current content. * @param other Array to copy into this array. It has to * have the same nb_component as this. If compiled in * debug mode, an incorrect other will result in an * exception being thrown. Optimised code may result in * unpredicted behaviour. * @param no_sanity_check turns off all checkes */ template void Array::copy(const Array & other, bool no_sanity_check) { AKANTU_DEBUG_IN(); if (not no_sanity_check and (other.nb_component != this->nb_component)) { AKANTU_ERROR("The two arrays do not have the same " "number of components"); } this->resize((other.size_ * other.nb_component) / this->nb_component); std::copy_n(other.data(), this->size_ * this->nb_component, this->values); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ template class ArrayPrintHelper { public: template static void print_content(const Array & vect, std::ostream & stream, int indent) { std::string space(indent, AKANTU_INDENT); stream << space << " + values : {"; for (Idx i = 0; i < vect.size(); ++i) { stream << "{"; for (Idx j = 0; j < vect.getNbComponent(); ++j) { stream << vect(i, j); if (j != vect.getNbComponent() - 1) { stream << ", "; } } stream << "}"; if (i != vect.size() - 1) { stream << ", "; } } stream << "}" << std::endl; } }; template <> class ArrayPrintHelper { public: template static void print_content(__attribute__((unused)) const Array & vect, __attribute__((unused)) std::ostream & stream, __attribute__((unused)) int indent) {} }; /* ------------------------------------------------------------------------ */ template void Array::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); std::streamsize prec = stream.precision(); std::ios_base::fmtflags ff = stream.flags(); stream.setf(std::ios_base::showbase); stream.precision(2); stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; stream << space << " + id : " << this->id << std::endl; stream << space << " + size : " << this->size_ << std::endl; stream << space << " + nb_component : " << this->nb_component << std::endl; stream << space << " + allocated size : " << this->getAllocatedSize() << std::endl; stream << space << " + memory size : " << printMemorySize(this->getMemorySize()) << std::endl; if (not AKANTU_DEBUG_LEVEL_IS_TEST()) { stream << space << " + address : " << std::hex << this->values << std::dec << std::endl; } stream.precision(prec); stream.flags(ff); if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) { ArrayPrintHelper::value>::print_content( *this, stream, indent); } stream << space << "]" << std::endl; } /* ------------------------------------------------------------------------ */ /* ArrayFilter */ /* ------------------------------------------------------------------------ */ template class ArrayFilter { public: /// const iterator template class const_iterator { public: Idx getCurrentIndex() { return (*this->filter_it * this->nb_item_per_elem + this->sub_element_counter); } inline const_iterator() = default; inline const_iterator(original_iterator origin_it, filter_iterator filter_it, Int nb_item_per_elem) : origin_it(std::move(origin_it)), filter_it(std::move(filter_it)), nb_item_per_elem(nb_item_per_elem), sub_element_counter(0){}; inline bool operator!=(const_iterator & other) const { return !((*this) == other); } inline bool operator==(const_iterator & other) const { return (this->origin_it == other.origin_it && this->filter_it == other.filter_it && this->sub_element_counter == other.sub_element_counter); } inline bool operator!=(const const_iterator & other) const { return !((*this) == other); } inline bool operator==(const const_iterator & other) const { return (this->origin_it == other.origin_it && this->filter_it == other.filter_it && this->sub_element_counter == other.sub_element_counter); } inline const_iterator & operator++() { ++sub_element_counter; if (sub_element_counter == nb_item_per_elem) { sub_element_counter = 0; ++filter_it; } return *this; }; inline decltype(auto) operator*() { return origin_it[nb_item_per_elem * (*filter_it) + sub_element_counter]; }; private: original_iterator origin_it; filter_iterator filter_it; /// the number of item per element Int nb_item_per_elem; /// counter for every sub element group Int sub_element_counter; }; using vector_iterator = void; // iterator>; using array_type = Array; using const_vector_iterator = const_iterator::const_scalar_iterator>; using value_type = typename array_type::value_type; private: /* ---------------------------------------------------------------------- */ /* Constructors/Destructors */ /* ---------------------------------------------------------------------- */ public: ArrayFilter(const Array & array, const Array & filter, Int nb_item_per_elem) : array(array), filter(filter), nb_item_per_elem(nb_item_per_elem){}; decltype(auto) begin_reinterpret(Int n, Int new_size) const { Int new_nb_item_per_elem = this->nb_item_per_elem; if (new_size != 0 && n != 0) new_nb_item_per_elem = this->array.getNbComponent() * this->filter.size() * this->nb_item_per_elem / (n * new_size); return const_vector_iterator(make_view(this->array, n).begin(), this->filter.begin(), new_nb_item_per_elem); }; decltype(auto) end_reinterpret(Int n, Int new_size) const { Int new_nb_item_per_elem = this->nb_item_per_elem; if (new_size != 0 && n != 0) new_nb_item_per_elem = this->array.getNbComponent() * this->filter.size() * this->nb_item_per_elem / (n * new_size); return const_vector_iterator(make_view(this->array, n).begin(), this->filter.end(), new_nb_item_per_elem); }; // vector_iterator begin_reinterpret(Int, Int) { throw; }; // vector_iterator end_reinterpret(Int, Int) { throw; }; /// return the size of the filtered array which is the filter size Int size() const { return this->filter.size() * this->nb_item_per_elem; }; /// the number of components of the filtered array Int getNbComponent() const { return this->array.getNbComponent(); }; /// tells if the container is empty bool empty() const [[gnu::warn_unused_result]] { return (size() == 0); } /* ---------------------------------------------------------------------- */ /* Class Members */ /* ---------------------------------------------------------------------- */ private: /// reference to array of data const Array & array; /// reference to the filter used to select elements const Array & filter; /// the number of item per element Int nb_item_per_elem; }; /* ------------------------------------------------------------------------ */ /* Begin/End functions implementation */ /* ------------------------------------------------------------------------ */ namespace detail { template constexpr auto take_front_impl(Tuple && t, std::index_sequence /*idxs*/) { return std::make_tuple(std::get(std::forward(t))...); } template constexpr auto take_front(Tuple && t) { return take_front_impl(std::forward(t), std::make_index_sequence{}); } template std::string to_string_all(T &&... t) { if (sizeof...(T) == 0) { return ""; } std::stringstream ss; bool noComma = true; ss << "("; (void)std::initializer_list{ (ss << (noComma ? "" : ", ") << t, noComma = false)...}; ss << ")"; return ss.str(); } template struct InstantiationHelper { template static auto instantiate(T && data, Ns... ns) { return std::make_unique(data, ns...); } }; template <> struct InstantiationHelper<0> { template static auto instantiate(T && data) { return data; } }; template decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) { const auto is_const_arr = std::is_const>::value; using type = ViewIteratorHelper_t; using iterator = std::conditional_t, view_iterator>; static_assert(sizeof...(Ns), "You should provide a least one size"); if (array.getNbComponent() * array.size() != Int(product_all(std::forward(ns)...))) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::ArrayException(), "The iterator on " << debug::demangle(typeid(Arr).name()) << to_string_all(array.size(), array.getNbComponent()) << "is not compatible with the type " << debug::demangle(typeid(type).name()) << to_string_all(ns...)); } return aka::apply([&](auto... n) { return iterator(data, n...); }, take_front(std::make_tuple(ns...))); } } // namespace detail /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ template template inline auto Array::begin(Ns &&... ns) { return detail::get_iterator(*this, this->values, std::forward(ns)..., this->size_); } template template inline auto Array::end(Ns &&... ns) { return detail::get_iterator(*this, this->values + this->nb_component * this->size_, std::forward(ns)..., this->size_); } template template inline auto Array::begin(Ns &&... ns) const { return detail::get_iterator(*this, this->values, std::forward(ns)..., this->size_); } template template inline auto Array::end(Ns &&... ns) const { return detail::get_iterator(*this, this->values + this->nb_component * this->size_, std::forward(ns)..., this->size_); } template template inline auto Array::begin_reinterpret(Ns &&... ns) { return detail::get_iterator(*this, this->values, std::forward(ns)...); } template template inline auto Array::end_reinterpret(Ns &&... ns) { return detail::get_iterator( *this, this->values + detail::product_all(std::forward(ns)...), std::forward(ns)...); } template template inline auto Array::begin_reinterpret(Ns &&... ns) const { return detail::get_iterator(*this, this->values, std::forward(ns)...); } template template inline auto Array::end_reinterpret(Ns &&... ns) const { return detail::get_iterator( *this, this->values + detail::product_all(std::forward(ns)...), std::forward(ns)...); } /* ------------------------------------------------------------------------ */ /* Views */ /* ------------------------------------------------------------------------ */ namespace detail { template class ArrayView { using tuple = std::tuple; public: using size_type = Idx; ~ArrayView() = default; constexpr ArrayView(Array && array, Ns... ns) noexcept : array(array), sizes(std::move(ns)...) {} constexpr ArrayView(const ArrayView & array_view) = default; constexpr ArrayView(ArrayView && array_view) noexcept = default; constexpr ArrayView & operator=(const ArrayView & array_view) = default; constexpr ArrayView & operator=(ArrayView && array_view) noexcept = default; auto begin() { return aka::apply( [&](auto &&... ns) { return detail::get_iterator(array.get(), array.get().data(), std::forward(ns)...); }, sizes); } auto begin() const { return aka::apply( [&](auto &&... ns) { return detail::get_iterator(array.get(), array.get().data(), std::forward(ns)...); }, sizes); } auto end() { return aka::apply( [&](auto &&... ns) { return detail::get_iterator( array.get(), array.get().data() + detail::product_all(std::forward(ns)...), std::forward(ns)...); }, sizes); } auto end() const { return aka::apply( [&](auto &&... ns) { return detail::get_iterator( array.get(), array.get().data() + detail::product_all(std::forward(ns)...), std::forward(ns)...); }, sizes); } auto cbegin() const { return this->begin(); } auto cend() const { return this->end(); } constexpr auto size() const { return std::get::value - 1>(sizes); } constexpr auto dims() const { return std::tuple_size::value - 1; } private: std::reference_wrapper> array; tuple sizes; }; /* ---------------------------------------------------------------------- */ template class ArrayView &, Ns...> { using tuple = std::tuple; public: constexpr ArrayView(const ArrayFilter & array, Ns... ns) : array(array), sizes(std::move(ns)...) {} constexpr ArrayView(const ArrayView & array_view) = default; constexpr ArrayView(ArrayView && array_view) = default; constexpr ArrayView & operator=(const ArrayView & array_view) = default; constexpr ArrayView & operator=(ArrayView && array_view) = default; auto begin() const { return aka::apply( [&](auto &&... ns) { return array.get().begin_reinterpret( std::forward(ns)...); }, sizes); } auto end() const { return aka::apply( [&](auto &&... ns) { return array.get().end_reinterpret( std::forward(ns)...); }, sizes); } auto cbegin() const { return this->begin(); } auto cend() const { return this->end(); } constexpr auto size() const { return std::get::value - 1>(sizes); } constexpr auto dims() const { return std::tuple_size::value - 1; } private: std::reference_wrapper> array; tuple sizes; }; } // namespace detail /* ------------------------------------------------------------------------ */ template decltype(auto) make_view(Array && array, const Ns... ns) { static_assert(aka::conjunction>...>::value, "Ns should be integral types"); AKANTU_DEBUG_ASSERT((detail::product_all(ns...) != 0), "You must specify non zero dimensions"); auto size = std::forward(array).size() * std::forward(array).getNbComponent() / detail::product_all(ns...); return detail::ArrayView..., std::common_type_t>( std::forward(array), std::move(ns)..., size); } } // namespace akantu //#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */ diff --git a/src/common/aka_tensor.hh b/src/common/aka_tensor.hh new file mode 100644 index 000000000..61746d601 --- /dev/null +++ b/src/common/aka_tensor.hh @@ -0,0 +1,462 @@ +/** + * @file aka_tensor.hh + * + * @author Nicolas Richart + * + * @date creation dim jan 23 2022 + * + * @brief A Documented file. + * + * @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 . + * + */ + +/* -------------------------------------------------------------------------- */ + +#ifndef AKA_TENSOR_HH_ +#define AKA_TENSOR_HH_ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +/* TensorBase */ +/* -------------------------------------------------------------------------- */ +template class TensorBase : public TensorTrait { + using RetType = TensorBase; + static_assert(ndim > 2, "TensorBase cannot by used for dimensions < 3"); + +protected: + using idx_t = Idx; + + template + using valid_args_t = typename std::enable_if< + aka::conjunction, + std::is_enum>...>::value and + ndim == sizeof...(Args), + int>::type; + +public: + using proxy = TensorBase; + using size_type = Idx; + template = 0> + TensorBase() { + n.fill(0); + } + + TensorBase() { n.fill(0); } + + template = 0> + constexpr TensorBase(Args... args) + : n{idx_t(args)...}, _size(detail::product_all(args...)) {} + + constexpr TensorBase(const TensorBase & other) + : n(other.n), _size(other._size), values(other.values) {} + + constexpr TensorBase(TensorBase && other) noexcept + : n(std::move(other.n)), _size(std::exchange(other._size, 0)), + values(std::exchange(other.values, nullptr)) {} + +protected: + template + constexpr auto check_indices( + const Array & idx, + std::integer_sequence /* for_template_deduction */) const { + bool result = true; + (void)std::initializer_list{(result &= idx[I] < n[I], 0)...}; + return result; + } + + template constexpr auto compute_index(Args... args) const { + std::array idx{idx_t(args)...}; + AKANTU_DEBUG_ASSERT( + check_indices(idx, + std::make_integer_sequence{}), + "The indexes are out of bound"); + idx_t index = 0, i = (sizeof...(Args)) - 1; + for (; i > 0; i--) { + index += idx[i]; + if (i > 0) { + index *= n[i - 1]; + } + } + return index + idx[0]; + } + + template + constexpr auto get_slice(idx_t s, std::index_sequence) { + return S(this->values + s * detail::product_all(n[I]...), n[I]...); + } + + template + constexpr auto get_slice(idx_t s, std::index_sequence) const { + return S(this->values + s * detail::product_all(n[I]...), n[I]...); + } + +public: + template = 0> + inline auto operator()(Args... args) -> T & { + return *(this->values + compute_index(std::move(args)...)); + } + + template = 0> + inline auto operator()(Args... args) const -> const T & { + return *(this->values + compute_index(std::move(args)...)); + } + + template < + class R = T, idx_t _ndim = ndim, + std::enable_if_t<(_ndim > 3) and not std::is_const::value> * = nullptr> + inline auto operator()(idx_t s) { + return get_slice>( + s, std::make_index_sequence()); + } + + template 3)> * = nullptr> + inline auto operator()(idx_t s) const { + return get_slice>( + s, std::make_index_sequence()); + } + + template ::value> * = + nullptr> + inline auto operator()(idx_t s) { + return get_slice< + Eigen::Map>>( + s, std::make_index_sequence()); + } + + template * = nullptr> + inline auto operator()(idx_t s) const { + return get_slice, Eigen::Dynamic, Eigen::Dynamic>>>( + s, std::make_index_sequence()); + } + +protected: + template auto transform(Operator && op) -> RetType & { + std::transform(this->values, this->values + this->_size, this->values, + std::forward(op)); + return *(static_cast(this)); + } + + template + auto transform(Other && other, Operator && op) -> RetType & { + AKANTU_DEBUG_ASSERT(_size == other.size(), + "The two tensors do not have the same size " + << this->_size << " != " << other._size); + + std::transform(this->values, this->values + this->_size, other.values, + this->values, std::forward(op)); + return *(static_cast(this)); + } + + template auto accumulate(T init, Operator && op) -> T { + return std::accumulate(this->values, this->values + this->_size, + std::move(init), std::forward(op)); + } + + template + auto transform_reduce(Other && other, T init, Accumulate && acc, + Operator && op) -> T { + return std::inner_product( + this->values, this->values + this->_size, other.data(), std::move(init), + std::forward(acc), std::forward(op)); + } + + // element wise arithmetic operators ----------------------------------------- +public: + inline decltype(auto) operator+=(const TensorBase & other) { + return transform(other, [](auto && a, auto && b) { return a + b; }); + } + + /* ------------------------------------------------------------------------ */ + inline auto operator-=(const TensorBase & other) -> TensorBase & { + return transform(other, [](auto && a, auto && b) { return a - b; }); + } + + /* ------------------------------------------------------------------------ */ + inline auto operator+=(const T & x) -> TensorBase & { + return transform([&x](auto && a) { return a + x; }); + } + + /* ------------------------------------------------------------------------ */ + inline auto operator-=(const T & x) -> TensorBase & { + return transform([&x](auto && a) { return a - x; }); + } + + /* ------------------------------------------------------------------------ */ + inline auto operator*=(const T & x) -> TensorBase & { + return transform([&x](auto && a) { return a * x; }); + } + + /* ---------------------------------------------------------------------- */ + inline auto operator/=(const T & x) -> TensorBase & { + return transform([&x](auto && a) { return a / x; }); + } + + /// Y = \alpha X + Y + inline auto aXplusY(const TensorBase & other, const T alpha = 1.) + -> TensorBase & { + return transform(other, + [&alpha](auto && a, auto && b) { return alpha * a + b; }); + } + + /* ------------------------------------------------------------------------ */ + auto data() -> T * { return values; } + auto data() const -> const T * { return values; } + + // clang-format off + [[deprecated("use data instead to be stl compatible")]] + auto storage() -> T*{ + return values; + } + + [[deprecated("use data instead to be stl compatible")]] + auto storage() const -> const T * { + return values; + } + // clang-format on + + auto size() const { return _size; } + auto size(idx_t i) const { + AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim + << " dimensions, not " + << (i + 1)); + return n[i]; + }; + + inline void set(const T & t) { std::fill_n(values, _size, t); }; + inline void clear() { set(T()); }; + +public: + /// "Entrywise" norm norm @f[ \|\boldsymbol{T}\|_p = \left( + /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}} + /// @f] + template * = nullptr> + auto lpNorm() const -> T { + return accumulate( + T(), [](auto && init, auto && a) { return init + std::abs(a); }); + } + + template * = nullptr> + auto lpNorm() const -> T { + return accumulate(T(), [](auto && init, auto && a) { + return std::max(init, std::abs(a)); + }); + } + + template * = nullptr> + auto norm() const -> T { + return std::sqrt( + accumulate(T(), [](auto && init, auto && a) { return init + a * a; })); + } + + template 2)> * = nullptr> + auto norm() const -> T { + return std::pow(accumulate(T(), + [](auto && init, auto && a) { + return init + std::pow(a, norm_type); + }), + 1. / norm_type); + } + + auto norm() const -> T { return lpNorm<2>(); } + +protected: + template = 0> + void serialize(std::ostream & stream, Args... args) const { + stream << this->operator()(std::move(args)...); + } + + template = 0> + void serialize(std::ostream & stream, Args... args) const { + stream << "["; + for (idx_t i = 0; i < n[N]; ++i) { + if (i != 0) { + stream << ","; + } + serialize(stream, std::move(args)..., i); + } + stream << "]"; + } + +public: + void printself(std::ostream & stream) const { serialize<0>(stream); }; + +protected: + template + constexpr decltype(auto) begin(std::index_sequence) { + return view_iterator>(values, + n[I]...); + } + + template + constexpr decltype(auto) end(std::index_sequence) { + return view_iterator>(values + _size, + n[I]...); + } + + template + constexpr decltype(auto) begin(std::index_sequence) const { + return const_view_iterator>(values, + n[I]...); + } + + template + constexpr decltype(auto) end(std::index_sequence) const { + return const_view_iterator>( + values + _size, n[I]...); + } + +public: + decltype(auto) begin() { return begin(std::make_index_sequence{}); } + decltype(auto) end() { return end(std::make_index_sequence{}); } + + decltype(auto) begin() const { + return begin(std::make_index_sequence{}); + } + decltype(auto) end() const { + return end(std::make_index_sequence{}); + } + +protected: + // size per dimension + std::array n; + + // total storage size + idx_t _size{0}; + + // actual data location + T * values{nullptr}; +}; + +/* -------------------------------------------------------------------------- */ +/* TensorProxy */ +/* -------------------------------------------------------------------------- */ +template class TensorProxy : public TensorBase { +private: + using parent = TensorBase; + +public: + // proxy constructor + template + constexpr TensorProxy(T * data, Args... args) : parent(args...) { + this->values = data; + } + + constexpr TensorProxy(const TensorProxy & other) : parent(other) { + this->values = other.values; + } + + constexpr TensorProxy(const Tensor & other) : parent(other) { + this->values = other.values; + } + + // move constructors --------------------------------------------------------- + // proxy -> proxy + TensorProxy(TensorProxy && other) noexcept : parent(other) {} + + auto operator=(const TensorBase & other) -> TensorProxy & { + AKANTU_DEBUG_ASSERT( + other.size() == this->size(), + "You are trying to copy too a tensors proxy with the wrong size " + << this->_size << " != " << other._size); + + static_assert(std::is_trivially_copyable{}, + "Cannot copy a tensor on non trivial types"); + + std::copy(other.values, other.values + this->_size, this->values); + return *this; + } +}; + +/* -------------------------------------------------------------------------- */ +/* Tensor */ +/* -------------------------------------------------------------------------- */ +template class Tensor : public TensorBase { +private: + using parent = TensorBase; + +public: + template constexpr Tensor(Args... args) : parent(args...) { + static_assert( + std::is_trivially_constructible{}, + "Cannot create a tensor on non trivially constructible types"); + this->values = new T[this->_size]; + } + + /* ------------------------------------------------------------------------ */ + virtual ~Tensor() { delete[] this->values; } + + // copy constructors --------------------------------------------------------- + constexpr Tensor(const Tensor & other) : parent(other) { + this->values = new T[this->_size]; + std::copy(other.values, other.values + this->_size, this->values); + } + + constexpr explicit Tensor(const TensorProxy & other) + : parent(other) { + // static_assert(false, "Copying data are you sure"); + this->values = new T[this->_size]; + std::copy(other.values, other.values + this->_size, this->values); + } + + // move constructors --------------------------------------------------------- + // proxy -> proxy, non proxy -> non proxy + Tensor(Tensor && other) noexcept : parent(other) {} + + // copy operator ------------------------------------------------------------- + /// operator= copy-and-swap + auto operator=(const TensorBase & other) -> Tensor & { + if (&other == this) + return *this; + + std::cout << "Warning: operator= delete data" << std::endl; + delete[] this->values; + this->n = other.n; + this->_size = other._size; + + static_assert( + std::is_trivially_constructible{}, + "Cannot create a tensor on non trivially constructible types"); + + this->values = new T[this->_size]; + + static_assert(std::is_trivially_copyable{}, + "Cannot copy a tensor on non trivial types"); + + std::copy(other.values, other.values + this->_size, this->values); + + return *this; + } +}; + +/* -------------------------------------------------------------------------- */ +template using Tensor3 = Tensor; +template using Tensor3Proxy = TensorProxy; +template using Tensor3Base = TensorBase; + +} // namespace akantu + +#endif // AKA_TENSOR_HH_ diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh index 4db44bf67..50d3d144e 100644 --- a/src/common/aka_types.hh +++ b/src/common/aka_types.hh @@ -1,871 +1,500 @@ /** * @file aka_types.hh * * @author Nicolas Richart * * @date creation: Thu Feb 17 2011 * @date last modification: Wed Dec 09 2020 * * @brief description of the "simple" types * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_compatibilty_with_cpp_standard.hh" #include "aka_error.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #ifndef AKANTU_AKA_TYPES_HH #define AKANTU_AKA_TYPES_HH /* -------------------------------------------------------------------------- */ namespace aka { template struct is_eigen_map : public std::false_type {}; } // namespace aka #define EIGEN_DEFAULT_DENSE_INDEX_TYPE akantu::Idx #define EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION Eigen::ColMajor /* -------------------------------------------------------------------------- */ #define EIGEN_MATRIXBASE_PLUGIN "aka_types_eigen_matrix_base_plugin.hh" #define EIGEN_MATRIX_PLUGIN "aka_types_eigen_matrix_plugin.hh" #define EIGEN_PLAINOBJECTBASE_PLUGIN \ "aka_types_eigen_plain_object_base_plugin.hh" #include #include /* -------------------------------------------------------------------------- */ namespace akantu { using Eigen::Ref; template using Vector = Eigen::Matrix; template using Matrix = Eigen::Matrix; template using VectorProxy = Eigen::Map::value, const Vector, n>, Vector, n>>>; template using MatrixProxy = Eigen::Map::value, const Matrix, m, n>, Matrix, m, n>>>; using VectorXr = Vector; using MatrixXr = Matrix; enum NormType : int8_t { L_1 = 1, L_2 = 2, L_inf = -1 }; struct TensorTraitBase {}; template struct TensorTrait : public TensorTraitBase {}; } // namespace akantu namespace aka { template using is_matrix = aka::bool_constant; template using is_vector = aka::bool_constant; template using are_vectors = aka::conjunction...>; template using are_matrices = aka::conjunction...>; template using enable_if_matrices_t = std::enable_if_t::value>; template using enable_if_vectors_t = std::enable_if_t::value>; -// template struct is_eigen_map : public std::false_type {}; - template struct is_eigen_map> : public std::true_type {}; -} // namespace aka - /* -------------------------------------------------------------------------- */ -namespace aka { template struct is_tensor : public std::is_base_of {}; template struct is_tensor> : public std::true_type {}; template struct is_tensor> : public std::true_type {}; template using is_tensor_n = std::is_base_of, T>; -// template using is_vector = is_tensor_n; - -// template using is_matrix = is_tensor_n; - template using enable_if_tensors_n = std::enable_if< aka::conjunction< is_tensor_n..., std::is_same< std::common_type_t...>, std::decay_t>...>::value, T>; template using enable_if_tensors_n_t = typename enable_if_tensors_n::type; -// template -// using enable_if_vectors_t = typename enable_if_tensors_n<1, T, -// Ms...>::type; template - -// using enable_if_matricies_t = typename enable_if_tensors_n<2, T, -// Ms...>::type; } // namespace aka -namespace akantu { +namespace akantu { // fwd declaration template class TensorBase; template class TensorProxy; template class Tensor; } // namespace akantu -#include "aka_view_iterators.hh" - -namespace akantu { - /* -------------------------------------------------------------------------- */ -template class TensorBase : public TensorTrait { - using RetType = TensorBase; - static_assert(ndim > 2, "TensorBase cannot by used for dimensions < 3"); - -protected: - using idx_t = Idx; - - template - using valid_args_t = typename std::enable_if< - aka::conjunction, - std::is_enum>...>::value and - ndim == sizeof...(Args), - int>::type; - -public: - using proxy = TensorBase; - using size_type = Idx; - template = 0> - TensorBase() { - n.fill(0); - } - - TensorBase() { n.fill(0); } - - template = 0> - constexpr TensorBase(Args... args) - : n{idx_t(args)...}, _size(detail::product_all(args...)) {} - - constexpr TensorBase(const TensorBase & other) - : n(other.n), _size(other._size), values(other.values) {} - - constexpr TensorBase(TensorBase && other) noexcept - : n(std::move(other.n)), _size(std::exchange(other._size, 0)), - values(std::exchange(other.values, nullptr)) {} - -protected: - template - constexpr auto check_indices( - const Array & idx, - std::integer_sequence /* for_template_deduction */) const { - bool result = true; - (void)std::initializer_list{(result &= idx[I] < n[I], 0)...}; - return result; - } - - template constexpr auto compute_index(Args... args) const { - std::array idx{idx_t(args)...}; - AKANTU_DEBUG_ASSERT( - check_indices(idx, - std::make_integer_sequence{}), - "The indexes are out of bound"); - idx_t index = 0, i = (sizeof...(Args)) - 1; - for (; i > 0; i--) { - index += idx[i]; - if (i > 0) { - index *= n[i - 1]; - } - } - return index + idx[0]; - } - - template - constexpr auto get_slice(idx_t s, std::index_sequence) { - return S(this->values + s * detail::product_all(n[I]...), n[I]...); - } - - template - constexpr auto get_slice(idx_t s, std::index_sequence) const { - return S(this->values + s * detail::product_all(n[I]...), n[I]...); - } - -public: - template = 0> - inline auto operator()(Args... args) -> T & { - return *(this->values + compute_index(std::move(args)...)); - } - - template = 0> - inline auto operator()(Args... args) const -> const T & { - return *(this->values + compute_index(std::move(args)...)); - } - - template < - class R = T, idx_t _ndim = ndim, - std::enable_if_t<(_ndim > 3) and not std::is_const::value> * = nullptr> - inline auto operator()(idx_t s) { - return get_slice>( - s, std::make_index_sequence()); - } - - template 3)> * = nullptr> - inline auto operator()(idx_t s) const { - return get_slice>( - s, std::make_index_sequence()); - } - - template ::value> * = - nullptr> - inline auto operator()(idx_t s) { - return get_slice< - Eigen::Map>>( - s, std::make_index_sequence()); - } - - template * = nullptr> - inline auto operator()(idx_t s) const { - return get_slice, Eigen::Dynamic, Eigen::Dynamic>>>( - s, std::make_index_sequence()); - } - -protected: - template auto transform(Operator && op) -> RetType & { - std::transform(this->values, this->values + this->_size, this->values, - std::forward(op)); - return *(static_cast(this)); - } - - template - auto transform(Other && other, Operator && op) -> RetType & { - AKANTU_DEBUG_ASSERT(_size == other.size(), - "The two tensors do not have the same size " - << this->_size << " != " << other._size); - - std::transform(this->values, this->values + this->_size, other.values, - this->values, std::forward(op)); - return *(static_cast(this)); - } - - template auto accumulate(T init, Operator && op) -> T { - return std::accumulate(this->values, this->values + this->_size, - std::move(init), std::forward(op)); - } - - template - auto transform_reduce(Other && other, T init, Accumulate && acc, - Operator && op) -> T { - return std::inner_product( - this->values, this->values + this->_size, other.data(), std::move(init), - std::forward(acc), std::forward(op)); - } - - // arithmetic operators ------------------------------------------------------ - /* ------------------------------------------------------------------------ */ -public: - inline decltype(auto) operator+=(const TensorBase & other) { - return transform(other, [](auto && a, auto && b) { return a + b; }); - } - - /* ------------------------------------------------------------------------ */ - inline auto operator-=(const TensorBase & other) -> TensorBase & { - return transform(other, [](auto && a, auto && b) { return a - b; }); - } - - /* ------------------------------------------------------------------------ */ - inline auto operator+=(const T & x) -> TensorBase & { - return transform([&x](auto && a) { return a + x; }); - } - - /* ------------------------------------------------------------------------ */ - inline auto operator-=(const T & x) -> TensorBase & { - return transform([&x](auto && a) { return a - x; }); - } - - /* ------------------------------------------------------------------------ */ - inline auto operator*=(const T & x) -> TensorBase & { - return transform([&x](auto && a) { return a * x; }); - } - - /* ---------------------------------------------------------------------- */ - inline auto operator/=(const T & x) -> TensorBase & { - return transform([&x](auto && a) { return a / x; }); - } - - /// Y = \alpha X + Y - inline auto aXplusY(const TensorBase & other, const T alpha = 1.) - -> TensorBase & { - return transform(other, - [&alpha](auto && a, auto && b) { return alpha * a + b; }); - } - - /* ------------------------------------------------------------------------ */ - auto data() -> T * { return values; } - auto data() const -> const T * { return values; } - - // clang-format off - [[deprecated("use data instead to be stl compatible")]] - auto storage() -> T*{ - return values; - } - - [[deprecated("use data instead to be stl compatible")]] - auto storage() const -> const T * { - return values; - } - // clang-format on - - auto size() const { return _size; } - auto size(idx_t i) const { - AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim - << " dimensions, not " - << (i + 1)); - return n[i]; - }; - - inline void set(const T & t) { std::fill_n(values, _size, t); }; - inline void clear() { set(T()); }; - -public: - /// "Entrywise" norm norm @f[ \|\boldsymbol{T}\|_p = \left( - /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}} - /// @f] - template - auto norm() -> std::enable_if_t const { - return accumulate( - T(), [](auto && init, auto && a) { return init + std::abs(a); }); - } - - template - auto norm() -> std::enable_if_t const { - return accumulate(T(), [](auto && init, auto && a) { - return std::max(init, std::abs(a)); - }); - } - - template - auto norm() -> std::enable_if_t const { - return std::sqrt( - accumulate(T(), [](auto && init, auto && a) { return init + a * a; })); - } - - template - auto norm() -> std::enable_if_t<(norm_type > L_2), T> const { - return std::pow(accumulate(T(), - [](auto && init, auto && a) { - return init + std::pow(a, norm_type); - }), - 1. / norm_type); - } - - auto norm() const -> T { return norm(); } - -protected: - template = 0> - void serialize(std::ostream & stream, Args... args) const { - stream << this->operator()(std::move(args)...); - } - - template = 0> - void serialize(std::ostream & stream, Args... args) const { - stream << "["; - for (idx_t i = 0; i < n[N]; ++i) { - if (i != 0) - stream << ","; - serialize(stream, std::move(args)..., i); - } - stream << "]"; - } - -public: - void printself(std::ostream & stream) const { serialize<0>(stream); }; - -protected: - template - constexpr decltype(auto) begin(std::index_sequence) { - return view_iterator>(values, - n[I]...); - } - - template - constexpr decltype(auto) end(std::index_sequence) { - return view_iterator>(values + _size, - n[I]...); - } - - template - constexpr decltype(auto) begin(std::index_sequence) const { - return const_view_iterator>(values, - n[I]...); - } - - template - constexpr decltype(auto) end(std::index_sequence) const { - return const_view_iterator>( - values + _size, n[I]...); - } - -public: - decltype(auto) begin() { return begin(std::make_index_sequence{}); } - decltype(auto) end() { return end(std::make_index_sequence{}); } - - decltype(auto) begin() const { - return begin(std::make_index_sequence{}); - } - decltype(auto) end() const { - return end(std::make_index_sequence{}); - } - -protected: - // size per dimension - std::array n; - - // total storage size - idx_t _size{0}; - - // actual data location - T * values{nullptr}; -}; - +#include "aka_view_iterators.hh" /* -------------------------------------------------------------------------- */ -template class TensorProxy : public TensorBase { -private: - using parent = TensorBase; - -public: - // proxy constructor - template - constexpr TensorProxy(T * data, Args... args) : parent(args...) { - this->values = data; - } - - constexpr TensorProxy(const TensorProxy & other) : parent(other) { - this->values = other.values; - } - - constexpr TensorProxy(const Tensor & other) : parent(other) { - this->values = other.values; - } - - // move constructors --------------------------------------------------------- - // proxy -> proxy - TensorProxy(TensorProxy && other) noexcept : parent(other) {} - - auto operator=(const TensorBase & other) -> TensorProxy & { - AKANTU_DEBUG_ASSERT( - other.size() == this->size(), - "You are trying to copy too a tensors proxy with the wrong size " - << this->_size << " != " << other._size); - - static_assert(std::is_trivially_copyable{}, - "Cannot copy a tensor on non trivial types"); - - std::copy(other.values, other.values + this->_size, this->values); - return *this; - } -}; - +#include "aka_tensor.hh" /* -------------------------------------------------------------------------- */ -template class Tensor : public TensorBase { -private: - using parent = TensorBase; - -public: - template constexpr Tensor(Args... args) : parent(args...) { - static_assert( - std::is_trivially_constructible{}, - "Cannot create a tensor on non trivially constructible types"); - this->values = new T[this->_size]; - } - /* ------------------------------------------------------------------------ */ - virtual ~Tensor() { delete[] this->values; } - - // copy constructors --------------------------------------------------------- - constexpr Tensor(const Tensor & other) : parent(other) { - this->values = new T[this->_size]; - std::copy(other.values, other.values + this->_size, this->values); - } - - constexpr explicit Tensor(const TensorProxy & other) - : parent(other) { - // static_assert(false, "Copying data are you sure"); - this->values = new T[this->_size]; - std::copy(other.values, other.values + this->_size, this->values); - } - - // move constructors --------------------------------------------------------- - // proxy -> proxy, non proxy -> non proxy - Tensor(Tensor && other) noexcept : parent(other) {} - - // copy operator ------------------------------------------------------------- - /// operator= copy-and-swap - auto operator=(const TensorBase & other) -> Tensor & { - if (&other == this) - return *this; - - std::cout << "Warning: operator= delete data" << std::endl; - delete[] this->values; - this->n = other.n; - this->_size = other._size; - - static_assert( - std::is_trivially_constructible{}, - "Cannot create a tensor on non trivially constructible types"); - - this->values = new T[this->_size]; - - static_assert(std::is_trivially_copyable{}, - "Cannot copy a tensor on non trivial types"); - - std::copy(other.values, other.values + this->_size, this->values); - - return *this; - } -}; - -/* -------------------------------------------------------------------------- */ -template using Tensor3 = Tensor; -template using Tensor3Proxy = TensorProxy; -template using Tensor3Base = TensorBase; +namespace akantu { class ArrayBase; /* -------------------------------------------------------------------------- */ namespace details { template struct MapPlainObjectType { using type = T; }; template struct MapPlainObjectType< Eigen::Map> { using type = PlainObjectType; }; template using MapPlainObjectType_t = typename MapPlainObjectType::type; template struct EigenMatrixViewHelper {}; template struct EigenMatrixViewHelper { using type = Eigen::Matrix; }; template struct EigenMatrixViewHelper { using type = Eigen::Matrix; }; template using EigenMatrixViewHelper_t = typename EigenMatrixViewHelper::type; template class EigenView { static_assert(sizeof...(sizes) == 1 or sizeof...(sizes) == 2, "Eigen only supports Vector and Matrices"); public: using size_type = typename std::decay_t::size_type; using value_type = typename std::decay_t::value_type; EigenView(Array && array, decltype(sizes)... sizes_) : array(array), sizes_(sizes_...) {} EigenView(Array && array) : array(array), sizes_(sizes...) {} EigenView(const EigenView & other) = default; EigenView(EigenView && other) noexcept = default; auto operator=(const EigenView & other) -> EigenView & = default; auto operator=(EigenView && other) noexcept -> EigenView & = default; template ().data())>, std::enable_if_t::value> * = nullptr> decltype(auto) begin() { return aka::make_from_tuple<::akantu::view_iterator< Eigen::Map>>>( std::tuple_cat(std::make_tuple(array.get().data()), sizes_)); } template ().data())>, std::enable_if_t::value> * = nullptr> decltype(auto) end() { return aka::make_from_tuple<::akantu::view_iterator< Eigen::Map>>>( std::tuple_cat(std::make_tuple(array.get().data() + array_size()), sizes_)); } decltype(auto) begin() const { return aka::make_from_tuple<::akantu::view_iterator< Eigen::Map>>>( std::tuple_cat(std::make_tuple(array.get().data()), sizes_)); } decltype(auto) end() const { return aka::make_from_tuple<::akantu::view_iterator< Eigen::Map>>>( std::tuple_cat(std::make_tuple(array.get().data() + array_size()), sizes_)); } private: template < class A = Array, std::enable_if_t>::value> * = nullptr> auto array_size() { return array.get().size() * array.get().getNbComponent(); } template >::value> * = nullptr> auto array_size() { return array.get().size(); } private: std::reference_wrapper> array; std::tuple sizes_; }; } // namespace details template decltype(auto) make_view(Array && array, Idx rows = RowsAtCompileTime) { return details::EigenView( std::forward(array), rows); } template decltype(auto) make_view(Array && array, Idx rows = RowsAtCompileTime, Idx cols = ColsAtCompileTime) { return details::EigenView( std::forward(array), rows, cols); } } // namespace akantu namespace Eigen { template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase::zero() { return this->fill(0.); } /* -------------------------------------------------------------------------- */ /* Vector */ /* -------------------------------------------------------------------------- */ template template ::value and ED::IsVectorAtCompileTime> *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::begin() { return ::akantu::view_iterator( this->derived().data()); } template template ::value and ED::IsVectorAtCompileTime> *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::end() { return ::akantu::view_iterator( this->derived().data() + this->size()); } template template *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::begin() const { using Scalar = typename Derived::Scalar; return ::akantu::const_view_iterator(this->derived().data()); } template template *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::end() const { using Scalar = typename Derived::Scalar; return ::akantu::const_view_iterator(this->derived().data() + this->size()); } /* -------------------------------------------------------------------------- */ /* Matrix */ /* -------------------------------------------------------------------------- */ template template ::value and not ED::IsVectorAtCompileTime> *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::begin() { return ::akantu::view_iterator< Map>>( this->derived().data(), this->rows()); } template template ::value and not ED::IsVectorAtCompileTime> *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::end() { return ::akantu::view_iterator< Map>>( this->derived().data() + this->size(), this->rows()); } template template *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::begin() const { using Scalar = typename Derived::Scalar; return ::akantu::const_view_iterator< Map>>( const_cast(this->derived().data()), this->rows()); } template template *> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) MatrixBase::end() const { using Scalar = typename Derived::Scalar; return ::akantu::const_view_iterator< Map>>( const_cast(this->derived().data()) + this->size(), this->rows()); } template template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase::eig(const MatrixBase & values_) const { EigenSolver>> solver(*this, false); - Eigen::MatrixBase & values = - const_cast &>( - values_); // as advised by the Eigen developers - - values = solver.eigenvalues(); + using OtherScalar = typename OtherDerived::Scalar; + + // as advised by the Eigen developers even though this is a UB + auto & values = const_cast &>(values_); + akantu::static_if(std::is_floating_point{}) + .then([&](auto && solver) { values = solver.eigenvalues().real(); }) + .else_([&](auto && solver) { values = solver.eigenvalues(); })( + std::forward(solver)); } template template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase::eig(const MatrixBase & values_, - const MatrixBase & vectors_) const { + const MatrixBase & vectors_, bool sort) const { EigenSolver>> - solver(*this, false); + solver(*this, true); - Eigen::MatrixBase & values = const_cast &>( - values_); // as advised by the Eigen developers - Eigen::MatrixBase & vectors = const_cast &>( - vectors_); // as advised by the Eigen developers + // as advised by the Eigen developers even though this is a UB + auto & values = const_cast &>(values_); + auto & vectors = const_cast &>(vectors_); - values = solver.eigenvalues(); - vectors = solver.eigenvectors(); + auto norm = this->norm(); + + using OtherScalar = typename D1::Scalar; + + if ((solver.eigenvectors().imag().template lpNorm() > + 1e-15 * norm) and + std::is_floating_point::value) { + AKANTU_EXCEPTION("This matrix has complex eigenvectors()"); + } + + if (not sort) { + akantu::static_if(std::is_floating_point{}) + .then([&](auto && solver) { + values = solver.eigenvalues().real(); + vectors = solver.eigenvectors().real(); + }) + .else_([&](auto && solver) { + values = solver.eigenvalues(); + vectors = solver.eigenvectors(); + })(std::forward(solver)); + return; + } + + if (not std::is_floating_point::value) { + AKANTU_EXCEPTION("Cannot sort complex eigen values"); + } + + values = solver.eigenvalues().real(); + + PermutationMatrix P(values.size()); + P.setIdentity(); + + std::sort(P.indices().data(), P.indices().data() + P.indices().size(), + [&values](const Index & a, const Index & b) { + return (values(a) - values(b)) > 0; + }); + + akantu::static_if(std::is_floating_point{}) + .then([&](auto && solver) { + values = P.transpose() * values; + vectors = solver.eigenvectors().real() * P; + })(std::forward(solver)); + return; } template template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase::eigh(const MatrixBase & values_) const { SelfAdjointEigenSolver< akantu::details::MapPlainObjectType_t>> - solver(*this, false); - Eigen::MatrixBase & values = - const_cast &>( - values_); // as advised by the Eigen developers + solver(*this, EigenvaluesOnly); + // as advised by the Eigen developers even though this is a UB + auto & values = const_cast &>(values_); values = solver.eigenvalues(); } template template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase::eigh(const MatrixBase & values_, - const MatrixBase & vectors_) const { + const MatrixBase & vectors_, bool sort) const { SelfAdjointEigenSolver< akantu::details::MapPlainObjectType_t>> - solver(*this, false); + solver(*this, ComputeEigenvectors); - Eigen::MatrixBase & values = const_cast &>( - values_); // as advised by the Eigen developers - Eigen::MatrixBase & vectors = const_cast &>( - vectors_); // as advised by the Eigen developers + // as advised by the Eigen developers, even though this is a UB + auto & values = const_cast &>(values_); + auto & vectors = const_cast &>(vectors_); + + if (not sort) { + values = solver.eigenvalues(); + vectors = solver.eigenvectors(); + return; + } values = solver.eigenvalues(); - vectors = solver.eigenvectors(); + PermutationMatrix P(values.size()); + P.setIdentity(); + + std::sort(P.indices().data(), P.indices().data() + P.indices().size(), + [&values](const Index & a, const Index & b) { + return (values(a) - values(b)) > 0; + }); + + values = P.transpose() * values(); + vectors = solver.eigenvectors() * P; } } // namespace Eigen namespace std { -template -struct is_convertible, - Eigen::Map> - : aka::bool_constant::value> {}; +template +struct is_convertible, + Eigen::Map> + : aka::bool_constant::value> {}; } // namespace std #endif /* AKANTU_AKA_TYPES_HH */ diff --git a/src/common/aka_types_eigen_matrix_base_plugin.hh b/src/common/aka_types_eigen_matrix_base_plugin.hh index d3f053016..8d1b7bd7c 100644 --- a/src/common/aka_types_eigen_matrix_base_plugin.hh +++ b/src/common/aka_types_eigen_matrix_base_plugin.hh @@ -1,263 +1,218 @@ using size_type = Index; template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE MatrixBase(std::initializer_list list) : MatrixBase() { Index i = 0; for (auto val : list) { this->operator()(i++) = val; } } template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE MatrixBase(std::initializer_list> list) : MatrixBase() { static_assert(std::is_trivially_copyable{}, "Cannot create a tensor on non trivial types"); Index m = list.size(); Index n = 0; for (auto row : list) { n = std::max(n, row.size()); } if (RowsAtCompileTime != -1 and RowsAtCompileTime != m) { AKANTU_EXCEPTION( "The size of the matrix does not correspond to the initializer_list"); } if (ColsAtCompileTime != -1 and ColsAtCompileTime != n) { AKANTU_EXCEPTION( "The size of the matrix does not correspond to the initializer_list"); } if (RowsAtCompileTime != -1 or ColsAtCompileTime != -1) { this->resize(m, n); } Index i = 0, j = 0; for (auto & row : list) { for (auto & val : row) { this->operator()(i, j++) = val; } ++i; j = 0; } } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void zero(); template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto operator()(Index c) { auto & d = this->derived(); return Map>( d.data() + c * d.outerStride(), d.outerStride()); } template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto operator()(Index c) const { const auto & d = this->derived(); return d.col(c); } template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index c) { return Base::operator()(c); } template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index c) const { return Base::operator()(c); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index i, Index j) { return Base::operator()(i, j); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) operator()(Index i, Index j) const { return Base::operator()(i, j); } template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin(); template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end(); template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin() const; template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end() const; template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and not ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin(); template < typename ED = Derived, typename T = std::remove_reference_t().data())>, std::enable_if_t::value and not ED::IsVectorAtCompileTime> * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end(); template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) begin() const; template * = nullptr> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto) end() const; // clang-format off [[deprecated("use data instead to be stl compatible")]] Scalar * storage() { return this->data(); } [[deprecated("use data instead to be stl compatible")]] const Scalar * storage() const { return this->data(); } // clang-format on -template * = nullptr> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto size() const { - return this->cols(); +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index size() const { + return this->rows() * this->cols(); } -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto size(Index i) const { +template * = nullptr> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index size(Index i) const { AKANTU_DEBUG_ASSERT(i < 2, "This tensor has only " << 2 << " dimensions, not " << (i + 1)); return (i == 0) ? this->rows() : this->cols(); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void set(const Scalar & t) { this->fill(t); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eye(const Scalar & t = Scalar()) { (*this) = t * Matrix::Identity(this->rows(), this->cols()); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void clear() { this->fill(Scalar()); }; template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto distance(const MatrixBase & other) const { return (*this - other).norm(); } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar doubleDot(const MatrixBase & other) const { Scalar sum = 0; eigen_assert(rows() == cols() and rows() == other.rows() and cols() == other.cols()); this->cwiseProduct(other).sum(); } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eig(const MatrixBase & other) const; template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void -eig(const MatrixBase & values, const MatrixBase & vectors) const; +eig(const MatrixBase & values, const MatrixBase & vectors, + bool sort = std::is_floating_point::value) const; template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eigh(const MatrixBase & other) const; template -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void -eigh(const MatrixBase & values, const MatrixBase & vectors) const; -/* -public: -*/ -// template ::value, int> = 0> -// inline bool equal(const VectorBase & v, -// R tolerance = Math::getTolerance()) const { -// T * a = this->data(); -// T * b = v.data(); -// idx_t i = 0; -// while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance)) -// ++i; -// return i == this->_size; -// } - -// template ::value, int> = 0> -// inline short compare(const TensorBase & v, -// Real tolerance = Math::getTolerance()) const { -// T * a = this->data(); -// T * b = v.data(); -// for (UInt i(0); i < this->_size; ++i, ++a, ++b) { -// if (std::abs(*a - *b) > tolerance) -// return (((*a - *b) > tolerance) ? 1 : -1); -// } -// return 0; -// } - -// template ::value, int> = 0> -// inline bool equal(const TensorBase & v) const { -// return std::equal(this->values, this->values + this->_size, v.data()); -// } - -// template ::value, int> = 0> -// inline short compare(const TensorBase & v) const { -// T * a = this->data(); -// T * b = v.data(); -// for (idx_t i(0); i < this->_size; ++i, ++a, ++b) { -// if (*a > *b) -// return 1; -// else if (*a < *b) -// return -1; -// } -// return 0; -// } +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void eigh(const MatrixBase & values, + const MatrixBase & vectors, + bool sort = true) const; template inline bool operator<=(const MatrixBase & v) const { return this->isMuchSmallerThan(v); } template inline bool operator>=(const MatrixBase & v) const { return v.isMuchSmallerThan(*this); } template inline bool operator<(const MatrixBase & v) const { return (*this <= v) and (*this != v); } template inline bool operator>(const MatrixBase & v) const { return (*this >= v) and (*this != v); } diff --git a/src/io/mesh_io/mesh_io_msh.cc b/src/io/mesh_io/mesh_io_msh.cc index 30c577c07..a2b3d2124 100644 --- a/src/io/mesh_io/mesh_io_msh.cc +++ b/src/io/mesh_io/mesh_io_msh.cc @@ -1,1124 +1,1118 @@ /** * @file mesh_io_msh.cc * * @author Dana Christen * @author Mauro Corrado * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Oct 29 2020 * * @brief Read/Write for MSH files generated by gmsh * * * @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 . * */ /* ----------------------------------------------------------------------------- Version (Legacy) 1.0 $NOD number-of-nodes node-number x-coord y-coord z-coord ... $ENDNOD $ELM number-of-elements elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list ... $ENDELM ----------------------------------------------------------------------------- Version 2.1 $MeshFormat version-number file-type data-size $EndMeshFormat $Nodes number-of-nodes node-number x-coord y-coord z-coord ... $EndNodes $Elements number-of-elements elm-number elm-type number-of-tags < tag > ... node-number-list ... $EndElements $PhysicalNames number-of-names physical-dimension physical-number "physical-name" ... $EndPhysicalNames $NodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... node-number value ... ... $EndNodeData $ElementData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number value ... ... $EndElementData $ElementNodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number number-of-nodes-per-element value ... ... $ElementEndNodeData ----------------------------------------------------------------------------- elem-type 1: 2-node line. 2: 3-node triangle. 3: 4-node quadrangle. 4: 4-node tetrahedron. 5: 8-node hexahedron. 6: 6-node prism. 7: 5-node pyramid. 8: 3-node second order line 9: 6-node second order triangle 10: 9-node second order quadrangle 11: 10-node second order tetrahedron 12: 27-node second order hexahedron 13: 18-node second order prism 14: 14-node second order pyramid 15: 1-node point. 16: 8-node second order quadrangle 17: 20-node second order hexahedron 18: 15-node second order prism 19: 13-node second order pyramid 20: 9-node third order incomplete triangle 21: 10-node third order triangle 22: 12-node fourth order incomplete triangle 23: 15-node fourth order triangle 24: 15-node fifth order incomplete triangle 25: 21-node fifth order complete triangle 26: 4-node third order edge 27: 5-node fourth order edge 28: 6-node fifth order edge 29: 20-node third order tetrahedron 30: 35-node fourth order tetrahedron 31: 56-node fifth order tetrahedron -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIOMSH::MeshIOMSH() { canReadSurface = true; canReadExtendedData = true; _msh_nodes_per_elem[_msh_not_defined] = 0; _msh_nodes_per_elem[_msh_segment_2] = 2; _msh_nodes_per_elem[_msh_triangle_3] = 3; _msh_nodes_per_elem[_msh_quadrangle_4] = 4; _msh_nodes_per_elem[_msh_tetrahedron_4] = 4; _msh_nodes_per_elem[_msh_hexahedron_8] = 8; _msh_nodes_per_elem[_msh_prism_1] = 6; _msh_nodes_per_elem[_msh_pyramid_1] = 1; _msh_nodes_per_elem[_msh_segment_3] = 3; _msh_nodes_per_elem[_msh_triangle_6] = 6; _msh_nodes_per_elem[_msh_quadrangle_9] = 9; _msh_nodes_per_elem[_msh_tetrahedron_10] = 10; _msh_nodes_per_elem[_msh_hexahedron_27] = 27; _msh_nodes_per_elem[_msh_hexahedron_20] = 20; _msh_nodes_per_elem[_msh_prism_18] = 18; _msh_nodes_per_elem[_msh_prism_15] = 15; _msh_nodes_per_elem[_msh_pyramid_14] = 14; _msh_nodes_per_elem[_msh_point] = 1; _msh_nodes_per_elem[_msh_quadrangle_8] = 8; _msh_to_akantu_element_types[_msh_not_defined] = _not_defined; _msh_to_akantu_element_types[_msh_segment_2] = _segment_2; _msh_to_akantu_element_types[_msh_triangle_3] = _triangle_3; _msh_to_akantu_element_types[_msh_quadrangle_4] = _quadrangle_4; _msh_to_akantu_element_types[_msh_tetrahedron_4] = _tetrahedron_4; _msh_to_akantu_element_types[_msh_hexahedron_8] = _hexahedron_8; _msh_to_akantu_element_types[_msh_prism_1] = _pentahedron_6; _msh_to_akantu_element_types[_msh_pyramid_1] = _not_defined; _msh_to_akantu_element_types[_msh_segment_3] = _segment_3; _msh_to_akantu_element_types[_msh_triangle_6] = _triangle_6; _msh_to_akantu_element_types[_msh_quadrangle_9] = _not_defined; _msh_to_akantu_element_types[_msh_tetrahedron_10] = _tetrahedron_10; _msh_to_akantu_element_types[_msh_hexahedron_27] = _not_defined; _msh_to_akantu_element_types[_msh_hexahedron_20] = _hexahedron_20; _msh_to_akantu_element_types[_msh_prism_18] = _not_defined; _msh_to_akantu_element_types[_msh_prism_15] = _pentahedron_15; _msh_to_akantu_element_types[_msh_pyramid_14] = _not_defined; _msh_to_akantu_element_types[_msh_point] = _point_1; _msh_to_akantu_element_types[_msh_quadrangle_8] = _quadrangle_8; _akantu_to_msh_element_types[_not_defined] = _msh_not_defined; _akantu_to_msh_element_types[_segment_2] = _msh_segment_2; _akantu_to_msh_element_types[_segment_3] = _msh_segment_3; _akantu_to_msh_element_types[_triangle_3] = _msh_triangle_3; _akantu_to_msh_element_types[_triangle_6] = _msh_triangle_6; _akantu_to_msh_element_types[_tetrahedron_4] = _msh_tetrahedron_4; _akantu_to_msh_element_types[_tetrahedron_10] = _msh_tetrahedron_10; _akantu_to_msh_element_types[_quadrangle_4] = _msh_quadrangle_4; _akantu_to_msh_element_types[_quadrangle_8] = _msh_quadrangle_8; _akantu_to_msh_element_types[_hexahedron_8] = _msh_hexahedron_8; _akantu_to_msh_element_types[_hexahedron_20] = _msh_hexahedron_20; _akantu_to_msh_element_types[_pentahedron_6] = _msh_prism_1; _akantu_to_msh_element_types[_pentahedron_15] = _msh_prism_15; _akantu_to_msh_element_types[_point_1] = _msh_point; #if defined(AKANTU_STRUCTURAL_MECHANICS) _akantu_to_msh_element_types[_bernoulli_beam_2] = _msh_segment_2; _akantu_to_msh_element_types[_bernoulli_beam_3] = _msh_segment_2; _akantu_to_msh_element_types[_discrete_kirchhoff_triangle_18] = _msh_triangle_3; #endif std::map::iterator it; for (it = _akantu_to_msh_element_types.begin(); it != _akantu_to_msh_element_types.end(); ++it) { UInt nb_nodes = _msh_nodes_per_elem[it->second]; std::vector tmp(nb_nodes); for (UInt i = 0; i < nb_nodes; ++i) { tmp[i] = i; } switch (it->first) { case _tetrahedron_10: tmp[8] = 9; tmp[9] = 8; break; case _pentahedron_6: tmp[0] = 2; tmp[1] = 0; tmp[2] = 1; tmp[3] = 5; tmp[4] = 3; tmp[5] = 4; break; case _pentahedron_15: tmp[0] = 2; tmp[1] = 0; tmp[2] = 1; tmp[3] = 5; tmp[4] = 3; tmp[5] = 4; tmp[6] = 8; tmp[8] = 11; tmp[9] = 6; tmp[10] = 9; tmp[11] = 10; tmp[12] = 14; tmp[14] = 12; break; case _hexahedron_20: tmp[9] = 11; tmp[10] = 12; tmp[11] = 9; tmp[12] = 13; tmp[13] = 10; tmp[17] = 19; tmp[18] = 17; tmp[19] = 18; break; default: // nothing to change break; } _read_order[it->first] = tmp; } } /* -------------------------------------------------------------------------- */ MeshIOMSH::~MeshIOMSH() = default; /* -------------------------------------------------------------------------- */ namespace { struct File { std::string filename; std::ifstream infile; std::string line; size_t current_line{0}; size_t first_node_number{std::numeric_limits::max()}, last_node_number{0}; bool has_physical_names{false}; std::unordered_map node_tags; std::unordered_map element_tags; double version{0}; int size_of_size_t{0}; Mesh & mesh; MeshAccessor mesh_accessor; std::multimap, int> entity_tag_to_physical_tags; File(const std::string & filename, Mesh & mesh) : filename(filename), mesh(mesh), mesh_accessor(mesh) { infile.open(filename.c_str()); if (not infile.good()) { AKANTU_EXCEPTION("Cannot open file " << filename); } } ~File() { infile.close(); } auto good() { return infile.good(); } std::stringstream get_line() { std::string tmp_str; if (infile.eof()) { AKANTU_EXCEPTION("Reached the end of the file " << filename); } std::getline(infile, tmp_str); line = trim(tmp_str); ++current_line; return std::stringstream(line); } template void read_line(Ts &&... ts) { auto && sstr = get_line(); (void)std::initializer_list{ (sstr >> std::forward(ts), 0)...}; } }; } // namespace /* -------------------------------------------------------------------------- */ template void MeshIOMSH::populateReaders2(File & file, Readers & readers) { readers["$NOD"] = readers["$Nodes"] = [&](const std::string & /*unused*/) { UInt nb_nodes; file.read_line(nb_nodes); Array & nodes = file.mesh_accessor.getNodes(); nodes.resize(nb_nodes); file.mesh_accessor.setNbGlobalNodes(nb_nodes); size_t index; Vector coord(3); /// for each node, read the coordinates for (auto && data : enumerate(make_view(nodes, nodes.getNbComponent()))) { file.read_line(index, coord(0), coord(1), coord(2)); if (index > std::numeric_limits::max()) { AKANTU_EXCEPTION( "There are more nodes in this files than the index type of akantu " "can handle, consider recompiling with a bigger index type"); } file.first_node_number = std::min(file.first_node_number, index); file.last_node_number = std::max(file.last_node_number, index); for (auto && coord_data : zip(std::get<1>(data), coord)) { std::get<0>(coord_data) = std::get<1>(coord_data); } file.node_tags[index] = std::get<0>(data); } }; readers["$ELM"] = readers["$Elements"] = [&](const std::string & /*unused*/) { Int nb_elements; file.read_line(nb_elements); Int index; UInt msh_type; ElementType akantu_type; for (UInt i = 0; i < nb_elements; ++i) { auto && sstr_elem = file.get_line(); sstr_elem >> index; sstr_elem >> msh_type; /// get the connectivity vector depending on the element type akantu_type = this->_msh_to_akantu_element_types[MSHElementType(msh_type)]; if (akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " << msh_type << " at line " << file.current_line); continue; } Element elem{akantu_type, 0, _not_ghost}; auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type); auto node_per_element = connectivity.getNbComponent(); auto & read_order = this->_read_order[akantu_type]; /// read tags informations if (file.version < 2) { Int tag0; Int tag1; Int nb_nodes; // reg-phys, reg-elem, number-of-nodes sstr_elem >> tag0 >> tag1 >> nb_nodes; auto & data0 = file.mesh_accessor.template getData("tag_0", akantu_type); data0.push_back(tag0); auto & data1 = file.mesh_accessor.template getData("tag_1", akantu_type); data1.push_back(tag1); } else if (file.version < 4) { Int nb_tags; sstr_elem >> nb_tags; for (Int j = 0; j < nb_tags; ++j) { Int tag; sstr_elem >> tag; auto & data = file.mesh_accessor.template getData( "tag_" + std::to_string(j), akantu_type); data.push_back(tag); } } Vector local_connect(node_per_element); for (Int j = 0; j < node_per_element; ++j) { Int node_index; sstr_elem >> node_index; AKANTU_DEBUG_ASSERT(node_index <= file.last_node_number, "Node number not in range : line " << file.current_line); local_connect(read_order[j]) = file.node_tags[node_index]; } connectivity.push_back(local_connect); elem.element = connectivity.size() - 1; file.element_tags[index] = elem; } }; readers["$Periodic"] = [&](const std::string &) { Int nb_periodic_entities; file.read_line(nb_periodic_entities); file.mesh_accessor.getNodesFlags().resize(file.mesh.getNbNodes(), NodeFlag::_normal); for (Int p = 0; p < nb_periodic_entities; ++p) { // dimension slave-tag master-tag Int dimension; file.read_line(dimension); // transformation file.get_line(); // nb nodes Int nb_nodes; file.read_line(nb_nodes); for (Int n = 0; n < nb_nodes; ++n) { // slave master auto && sstr = file.get_line(); // The info in the mesh seem inconsistent so they are ignored for now. continue; if (dimension == file.mesh.getSpatialDimension() - 1) { Idx slave, master; sstr >> slave; sstr >> master; file.mesh_accessor.addPeriodicSlave(file.node_tags[slave], file.node_tags[master]); } } } // mesh_accessor.markMeshPeriodic(); }; } /* -------------------------------------------------------------------------- */ template void MeshIOMSH::populateReaders4(File & file, Readers & readers) { static std::map entity_type{ {0, "points"}, {1, "curve"}, {2, "surface"}, {3, "volume"}, }; readers["$Entities"] = [&](const std::string & /*unused*/) { size_t num_entity[4]; file.read_line(num_entity[0], num_entity[1], num_entity[2], num_entity[3]); for (auto entity_dim : arange(4)) { for (auto _ [[gnu::unused]] : arange(num_entity[entity_dim])) { auto && sstr = file.get_line(); int tag; double min_x; double min_y; double min_z; double max_x; double max_y; double max_z; size_t num_physical_tags; sstr >> tag >> min_x >> min_y >> min_z; if (entity_dim > 0 or file.version < 4.1) { sstr >> max_x >> max_y >> max_z; } sstr >> num_physical_tags; for (auto _ [[gnu::unused]] : arange(num_physical_tags)) { int phys_tag; sstr >> phys_tag; std::string physical_name; if (this->physical_names.find(phys_tag) == this->physical_names.end()) { physical_name = "msh_block_" + std::to_string(phys_tag); } else { physical_name = this->physical_names[phys_tag]; } if (not file.mesh.elementGroupExists(physical_name)) { file.mesh.createElementGroup(physical_name, entity_dim); } else { file.mesh.getElementGroup(physical_name).addDimension(entity_dim); } file.entity_tag_to_physical_tags.insert( std::make_pair(std::make_pair(tag, entity_dim), phys_tag)); } } } }; readers["$Nodes"] = [&](const std::string & /*unused*/) { size_t num_blocks; size_t num_nodes; if (file.version >= 4.1) { file.read_line(num_blocks, num_nodes, file.first_node_number, file.last_node_number); } else { file.read_line(num_blocks, num_nodes); } auto & nodes = file.mesh_accessor.getNodes(); nodes.reserve(num_nodes); file.mesh_accessor.setNbGlobalNodes(num_nodes); if (num_nodes > std::numeric_limits::max()) { AKANTU_EXCEPTION( "There are more nodes in this files than the index type of akantu " "can handle, consider recompiling with a bigger index type"); } size_t node_id{0}; + auto dim = nodes.getNbComponent(); + for (auto block [[gnu::unused]] : arange(num_blocks)) { int entity_dim; int entity_tag; int parametric; size_t num_nodes_in_block; Vector pos(3); - Vector real_pos(nodes.getNbComponent()); if (file.version >= 4.1) { file.read_line(entity_dim, entity_tag, parametric, num_nodes_in_block); if (parametric) { AKANTU_EXCEPTION( "Akantu does not support parametric nodes in msh files"); } for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { size_t tag; file.read_line(tag); file.node_tags[tag] = node_id; ++node_id; } for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { file.read_line(pos(_x), pos(_y), pos(_z)); - for (auto && data : zip(real_pos, pos)) { - std::get<0>(data) = std::get<1>(data); - } - - nodes.push_back(real_pos); + nodes.push_back(pos.block(0, 0, dim, 1)); } } else { file.read_line(entity_tag, entity_dim, parametric, num_nodes_in_block); for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { size_t tag; file.read_line(tag, pos(_x), pos(_y), pos(_z)); if (file.version < 4.1) { file.first_node_number = std::min(file.first_node_number, tag); file.last_node_number = std::max(file.last_node_number, tag); } - for (auto && data : zip(real_pos, pos)) { - std::get<0>(data) = std::get<1>(data); - } + nodes.push_back(pos.block(0, 0, dim, 1)); - nodes.push_back(real_pos); file.node_tags[tag] = node_id; ++node_id; } } } }; readers["$Elements"] = [&](const std::string & /*unused*/) { size_t num_blocks; size_t num_elements; file.read_line(num_blocks, num_elements); for (auto block [[gnu::unused]] : arange(num_blocks)) { int entity_dim; int entity_tag; int element_type; size_t num_elements_in_block; if (file.version >= 4.1) { file.read_line(entity_dim, entity_tag, element_type, num_elements_in_block); } else { file.read_line(entity_tag, entity_dim, element_type, num_elements_in_block); } /// get the connectivity vector depending on the element type auto && akantu_type = this->_msh_to_akantu_element_types[(MSHElementType)element_type]; if (akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " << element_type << " at line " << file.current_line); continue; } Element elem{akantu_type, 0, _not_ghost}; auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type); Vector local_connect(connectivity.getNbComponent()); auto && read_order = this->_read_order[akantu_type]; auto & data0 = file.mesh_accessor.template getData("tag_0", akantu_type); data0.resize(data0.size() + num_elements_in_block, 0); auto & physical_data = file.mesh_accessor.template getData( "physical_names", akantu_type); physical_data.resize(physical_data.size() + num_elements_in_block, ""); for (auto _ [[gnu::unused]] : arange(num_elements_in_block)) { auto && sstr_elem = file.get_line(); std::size_t elem_tag; sstr_elem >> elem_tag; for (auto && c : arange(connectivity.getNbComponent())) { std::size_t node_tag; sstr_elem >> node_tag; AKANTU_DEBUG_ASSERT(node_tag <= file.last_node_number, "Node number not in range : line " << file.current_line); node_tag = file.node_tags[node_tag]; local_connect(read_order[c]) = node_tag; } connectivity.push_back(local_connect); elem.element = connectivity.size() - 1; file.element_tags[elem_tag] = elem; auto range = file.entity_tag_to_physical_tags.equal_range( std::make_pair(entity_tag, entity_dim)); bool first = true; for (auto it = range.first; it != range.second; ++it) { auto phys_it = this->physical_names.find(it->second); if (first) { data0(elem.element) = it->second; // for compatibility with version 2 if (phys_it != this->physical_names.end()) { physical_data(elem.element) = phys_it->second; } first = false; } if (phys_it != this->physical_names.end()) { file.mesh.getElementGroup(phys_it->second).add(elem, true, false); } } } } for (auto && element_group : file.mesh.iterateElementGroups()) { element_group.getNodeGroup().optimize(); } }; } /* -------------------------------------------------------------------------- */ void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { File file(filename, mesh); std::map> readers; readers["$MeshFormat"] = [&](const std::string & /*unused*/) { auto && sstr = file.get_line(); int format; sstr >> file.version >> format; if (format != 0) { AKANTU_ERROR("This reader can only read ASCII files."); } if (file.version > 2) { sstr >> file.size_of_size_t; if (file.size_of_size_t > int(sizeof(UInt))) { AKANTU_DEBUG_INFO("The size of the indexes in akantu might be to small " "to read this file (akantu " << sizeof(UInt) << " vs. msh file " << file.size_of_size_t << ")"); } } if (file.version < 4) { this->populateReaders2(file, readers); } else { this->populateReaders4(file, readers); } }; auto && read_data = [&](auto && entity_tags, auto && get_data, auto && read_data) { auto read_data_tags = [&](auto x) { UInt number_of_tags{0}; file.read_line(number_of_tags); std::vector tags(number_of_tags); for (auto && tag : tags) { file.read_line(tag); } return tags; }; auto && string_tags = read_data_tags(std::string{}); auto && real_tags [[gnu::unused]] = read_data_tags(double{}); auto && int_tags = read_data_tags(int{}); for (auto & s : string_tags) { s = trim(s, '"'); } auto id = string_tags[0]; auto size = int_tags[2]; auto nb_component = int_tags[1]; auto & data = get_data(id, size, nb_component); for (auto n [[gnu::unused]] : arange(size)) { auto && sstr = file.get_line(); size_t tag; sstr >> tag; const auto & entity = entity_tags[tag]; read_data(entity, sstr, data, nb_component); } }; readers["$NodeData"] = [&](const std::string & /*unused*/) { /* $NodeData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... nodeTag(size_t) value(double) ... ... $EndNodeData */ read_data( file.node_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> Array & { auto & data = file.mesh.template getNodalData(id, nb_component); data.resize(size); return data; }, [&](auto && node, auto && sstr, auto && data, auto && nb_component [[gnu::unused]]) { for (auto c : arange(nb_component)) { sstr >> data(node, c); } }); }; readers["$ElementData"] = [&](const std::string & /*unused*/) { /* $ElementData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... elementTag(size_t) value(double) ... ... $EndElementData */ read_data( file.element_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> ElementTypeMapArray & { file.mesh.template getElementalData(id); return file.mesh.template getElementalData(id); }, [&](auto && element, auto && sstr, auto && data, auto && nb_component) { if (not data.exists(element.type)) { data.alloc(mesh.getNbElement(element.type), nb_component, element.type, element.ghost_type); } auto & data_array = data(element.type); for (auto c : arange(nb_component)) { sstr >> data_array(element.element, c); } }); }; readers["$ElementNodeData"] = [&](const std::string & /*unused*/) { /* $ElementNodeData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... elementTag(size_t) value(double) ... ... $EndElementNodeData */ read_data( file.element_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> ElementTypeMapArray & { file.mesh.template getElementalData(id); auto & data = file.mesh.template getElementalData(id); data.isNodal(true); return data; }, [&](auto && element, auto && sstr, auto && data, auto && nb_component) { int nb_nodes_per_element; sstr >> nb_nodes_per_element; if (not data.exists(element.type)) { data.alloc(mesh.getNbElement(element.type), nb_component * nb_nodes_per_element, element.type, element.ghost_type); } auto & data_array = data(element.type); for (auto c : arange(nb_component)) { sstr >> data_array(element.element, c); } }); }; readers["$PhysicalNames"] = [&](const std::string & /*unused*/) { file.has_physical_names = true; int num_of_phys_names; file.read_line(num_of_phys_names); /// the format line for (auto k [[gnu::unused]] : arange(num_of_phys_names)) { int phys_name_id; int phys_dim; std::string phys_name; file.read_line(phys_dim, phys_name_id, std::quoted(phys_name)); this->physical_names[phys_name_id] = phys_name; } }; readers["Unsupported"] = [&](const std::string & _block) { std::string block = _block.substr(1); AKANTU_DEBUG_WARNING("Unsupported block_kind " << block << " at line " << file.current_line); auto && end_block = "$End" + block; while (file.line != end_block) { file.get_line(); } }; while (file.good()) { std::string block; file.read_line(block); auto && it = readers.find(block); if (it != readers.end()) { it->second(block); std::string end_block; file.read_line(end_block); block = block.substr(1); if (end_block != "$End" + block) { AKANTU_EXCEPTION("The reader failed to properly read the block " << block << ". Expected a $End" << block << " at line " << file.current_line); } } else if (not block.empty()) { readers["Unsupported"](block); } } // mesh.updateTypesOffsets(_not_ghost); if (file.version < 4) { this->constructPhysicalNames("tag_0", mesh); if (file.has_physical_names) { mesh.createGroupsFromMeshData("physical_names"); } } MeshUtils::fillElementToSubElementsData(mesh); } /* -------------------------------------------------------------------------- */ void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) { std::ofstream outfile; const Array & nodes = mesh.getNodes(); outfile.open(filename.c_str()); outfile << "$MeshFormat" << "\n"; outfile << "2.2 0 8" << "\n"; outfile << "$EndMeshFormat" << "\n"; outfile << std::setprecision(std::numeric_limits::digits10); outfile << "$Nodes" << "\n"; outfile << nodes.size() << "\n"; outfile << std::uppercase; for (Int i = 0; i < nodes.size(); ++i) { auto offset = i * nodes.getNbComponent(); outfile << i + 1; for (Int j = 0; j < nodes.getNbComponent(); ++j) { outfile << " " << nodes.data()[offset + j]; } for (UInt p = nodes.getNbComponent(); p < 3; ++p) { outfile << " " << 0.; } outfile << "\n"; ; } outfile << std::nouppercase; outfile << "$EndNodes" << "\n"; outfile << "$Elements" << "\n"; Int nb_elements = 0; for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { const auto & connectivity = mesh.getConnectivity(type, _not_ghost); nb_elements += connectivity.size(); } outfile << nb_elements << "\n"; std::map element_to_msh_element; Idx element_idx = 1; auto element = ElementNull; for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { const auto & connectivity = mesh.getConnectivity(type, _not_ghost); element.type = type; UInt * tag[2] = {nullptr, nullptr}; if (mesh.hasData("tag_0", type, _not_ghost)) { const auto & data_tag_0 = mesh.getData("tag_0", type, _not_ghost); tag[0] = data_tag_0.data(); } if (mesh.hasData("tag_1", type, _not_ghost)) { const auto & data_tag_1 = mesh.getData("tag_1", type, _not_ghost); tag[1] = data_tag_1.data(); } for (auto && data : enumerate(make_view(connectivity, connectivity.getNbComponent()))) { element.element = std::get<0>(data); const auto & conn = std::get<1>(data); element_to_msh_element.insert(std::make_pair(element, element_idx)); outfile << element_idx << " " << _akantu_to_msh_element_types[type] << " 2"; /// \todo write the real data in the file for (UInt t = 0; t < 2; ++t) { if (tag[t] != nullptr) { outfile << " " << tag[t][element.element]; } else { outfile << " 0"; } } for (auto && c : conn) { outfile << " " << c + 1; } outfile << "\n"; element_idx++; } } outfile << "$EndElements" << "\n"; if (mesh.hasData(MeshDataType::_nodal)) { auto && tags = mesh.getTagNames(); for (auto && tag : tags) { auto type = mesh.getTypeCode(tag, MeshDataType::_nodal); if (type != MeshDataTypeCode::_real) { AKANTU_DEBUG_WARNING( "The field " << tag << " is ignored by the MSH writer, msh files do not support " << type << " data"); continue; } auto && data = mesh.getNodalData(tag); outfile << "$NodeData" << "\n"; outfile << "1" << "\n"; outfile << "\"" << tag << "\"\n"; outfile << "1\n0.0" << "\n"; outfile << "3\n0" << "\n"; outfile << data.getNbComponent() << "\n"; outfile << data.size() << "\n"; for (auto && d : enumerate(make_view(data, data.getNbComponent()))) { outfile << std::get<0>(d) + 1; for (auto && v : std::get<1>(d)) { outfile << " " << v; } outfile << "\n"; } outfile << "$EndNodeData" << "\n"; } } if (mesh.hasData(MeshDataType::_elemental)) { auto && tags = mesh.getTagNames(); for (auto && tag : tags) { auto && data = mesh.getElementalData(tag); auto type = mesh.getTypeCode(tag, MeshDataType::_elemental); if (type != MeshDataTypeCode::_real) { AKANTU_DEBUG_WARNING( "The field " << tag << " is ignored by the MSH writer, msh files do not support " << type << " data"); continue; } if (data.isNodal()) { continue; } auto size = data.size(); if (size == 0) { continue; } auto && nb_components = data.getNbComponents(); auto nb_component = nb_components(*(data.elementTypes().begin())); outfile << "$ElementData" << "\n"; outfile << "1" << "\n"; outfile << "\"" << tag << "\"\n"; outfile << "1\n0.0" << "\n"; outfile << "3\n0" << "\n"; outfile << nb_component << "\n"; outfile << size << "\n"; Element element; for (auto type : data.elementTypes()) { element.type = type; for (auto && _ : enumerate(make_view(data(type), nb_components(type)))) { element.element = std::get<0>(_); outfile << element_to_msh_element[element]; for (auto && v : std::get<1>(_)) { outfile << " " << v; } outfile << "\n"; } } outfile << "$EndElementData" << "\n"; } } outfile.close(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/common/non_local_toolbox/base_weight_function.cc b/src/model/common/non_local_toolbox/base_weight_function.cc new file mode 100644 index 000000000..9f4491b8c --- /dev/null +++ b/src/model/common/non_local_toolbox/base_weight_function.cc @@ -0,0 +1,36 @@ +/** + * @file base_weight_function.cc + * + * @author Nicolas Richart + * + * @date creation sam jan 22 2022 + * + * @brief A Documented file. + * + * @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 . + * + */ + +#include "base_weight_function.hh" + +namespace akantu { + +INSTANTIATE_NL_NEIGHBORHOOD(base_wf, BaseWeightFunction); + +} // namespace akantu diff --git a/src/model/common/non_local_toolbox/base_weight_function_inline_impl.hh b/src/model/common/non_local_toolbox/base_weight_function_inline_impl.hh index 983f784ed..b0b0a062a 100644 --- a/src/model/common/non_local_toolbox/base_weight_function_inline_impl.hh +++ b/src/model/common/non_local_toolbox/base_weight_function_inline_impl.hh @@ -1,76 +1,74 @@ /** * @file base_weight_function_inline_impl.hh * * @author Nicolas Richart * @author Cyprien Wolff * * @date creation: Wed Sep 01 2010 * @date last modification: Wed Sep 27 2017 * * @brief Implementation of inline function of base weight function * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "base_weight_function.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_BASE_WEIGHT_FUNCTION_INLINE_IMPL_HH_ #define AKANTU_BASE_WEIGHT_FUNCTION_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline void BaseWeightFunction::init() { /// compute R^2 for a given non-local radius this->R2 = this->R * this->R; } /* -------------------------------------------------------------------------- */ inline void BaseWeightFunction::setRadius(Real radius) { /// set the non-local radius and update R^2 accordingly this->R = radius; this->R2 = this->R * this->R; } /* -------------------------------------------------------------------------- */ inline Real BaseWeightFunction::operator()(Real r, const IntegrationPoint & /* q1 */, const IntegrationPoint & /* q2 */) const { /// initialize the weight Real w = 0; /// compute weight for given r if (r <= this->R) { Real alpha = (1. - r * r / this->R2); w = alpha * alpha; // *weight = 1 - sqrt(r / radius); } return w; } -INSTANTIATE_NL_NEIGHBORHOOD(base_wf, BaseWeightFunction); - } // namespace akantu #endif /* AKANTU_BASE_WEIGHT_FUNCTION_INLINE_IMPL_HH_ */ diff --git a/src/model/phase_field/phasefields/phasefield_exponential.hh b/src/model/phase_field/phasefields/phasefield_exponential.hh index 088cbe33a..da0fd55a5 100644 --- a/src/model/phase_field/phasefields/phasefield_exponential.hh +++ b/src/model/phase_field/phasefields/phasefield_exponential.hh @@ -1,133 +1,133 @@ /** * @file phasefield_exponential.hh * * @author Mohit Pundir * * @date creation: Fri Jun 19 2020 * @date last modification: Wed Jun 23 2021 * * @brief Phasefield law for approximating discrete crack as an exponential * * * @section LICENSE * * Copyright (©) 2018-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 . * */ /* -------------------------------------------------------------------------- */ #include "phasefield.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ #define __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ namespace akantu { class PhaseFieldExponential : public PhaseField { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PhaseFieldExponential(PhaseFieldModel & model, const ID & id = ""); ~PhaseFieldExponential() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: void computePhiOnQuad(const Matrix & /*strain_quad*/, Real & /*phi_quad*/, Real & /*phi_hist_quad*/); void computeDrivingForce(ElementType /*el_type*/, GhostType /*ghost_type*/) override; inline void computeDrivingForceOnQuad(const Real & /*phi_quad*/, Real & /*driving_force_quad*/); inline void computeDamageEnergyDensityOnQuad(const Real & /*phi_quad*/, Real & /*dam_energy_quad*/); public: void updateInternalParameters() override; }; /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computeDrivingForceOnQuad(const Real & phi_quad, Real & driving_force_quad) { driving_force_quad = 2.0 * phi_quad; } /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computeDamageEnergyDensityOnQuad( const Real & phi_quad, Real & dam_energy_quad) { dam_energy_quad = 2.0 * phi_quad + this->g_c / this->l0; } /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computePhiOnQuad(const Matrix & strain_quad, Real & phi_quad, Real & phi_hist_quad) { Matrix strain_plus(spatial_dimension, spatial_dimension); // Matrix strain_minus(spatial_dimension, spatial_dimension); Matrix strain_dir(spatial_dimension, spatial_dimension); Matrix strain_diag_plus(spatial_dimension, spatial_dimension); // Matrix strain_diag_minus(spatial_dimension, spatial_dimension); Vector strain_values(spatial_dimension); Real trace_plus, trace_minus; strain_plus.zero(); // strain_minus.zero(); strain_dir.zero(); strain_values.zero(); strain_diag_plus.zero(); // strain_diag_minus.zero(); - strain_quad.eigh(strain_values, strain_dir); + strain_quad.eig(strain_values, strain_dir); for (UInt i = 0; i < spatial_dimension; i++) { strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); // strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); } strain_plus = strain_dir * strain_diag_plus * strain_dir.transpose(); // Nicolas: @mohit is the second transpose really here ? // strain_minus = strain_dir * strain_diag_minus * strain_dir.transpose(); trace_plus = std::max(Real(0.), strain_quad.trace()); // trace_minus = std::min(Real(0.), strain_quad.trace()); auto I = Matrix::Identity(spatial_dimension, spatial_dimension); Matrix sigma_plus = I * lambda * trace_plus + 2 * mu * strain_plus; // Matrix sigma_minus = I * lambda * trace_minus + 2 * mu * // strain_minus; phi_quad = 0.5 * sigma_plus.doubleDot(strain_quad); if (phi_quad < phi_hist_quad) { phi_quad = phi_hist_quad; } } } // namespace akantu #endif diff --git a/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh b/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh index 04d0ecfa0..adb149b06 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh @@ -1,378 +1,378 @@ /** * @file material_anisotropic_damage_tmpl.hh * * @author Emil Gallyamov * @author Nicolas Richart * * @date creation: Wed Jul 03 2019 * @date last modification: Fri Jul 24 2020 * * @brief Base class for anisotropic damage materials * * * @section LICENSE * * Copyright (©) 2018-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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "material_anisotropic_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_ #define AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_ namespace akantu { struct EmptyIteratorContainer { using size_type = std::size_t; struct iterator { auto & operator++() { return *this; } Real operator*() { return 0; } bool operator!=(const iterator & /*unused*/) const { return true; } bool operator==(const iterator & /*unused*/) const { return false; } }; auto begin() const // NOLINT(readability-convert-member-functions-to-static) { return iterator(); } auto end() const // NOLINT(readability-convert-member-functions-to-static) { return iterator(); } }; } // namespace akantu namespace std { template <> struct iterator_traits<::akantu::EmptyIteratorContainer::iterator> { using iterator_category = forward_iterator_tag; using value_type = akantu::Real; using difference_type = std::ptrdiff_t; using pointer = akantu::Real *; using reference = akantu::Real &; }; } // namespace std namespace akantu { namespace { template void tensorPlus_(const Eigen::MatrixBase & A, Op && oper) { Vector A_eigs; - A.eigh(A_eigs); + A.eig(A_eigs); for (auto & ap : A_eigs) { oper(ap); } } template auto tensorPlus2(const Eigen::MatrixBase & A) { Real square = 0; tensorPlus_(A, [&](Real eig) { eig = std::max(eig, 0.); square += eig * eig; }); return square; } template auto tensorPlusTrace(const Eigen::MatrixBase & A) { Real trace_plus = 0.; Real trace_minus = 0.; tensorPlus_(A, [&](Real eig) { trace_plus += std::max(eig, 0.); trace_minus += std::min(eig, 0.); }); return std::make_pair(trace_plus, trace_minus); } template auto tensorPlusOp(const Eigen::MatrixBase & A, Eigen::MatrixBase & A_directions, Op && oper) { Vector A_eigs; Matrix A_diag; - A.eigh(A_eigs, A_directions); + A.eig(A_eigs, A_directions); for (auto && data : enumerate(A_eigs)) { auto i = std::get<0>(data); A_diag(i, i) = oper(std::max(std::get<1>(data), 0.), i); } return A_directions * A_diag * A_directions.transpose(); } template auto tensorPlus(const Eigen::MatrixBase & A, Eigen::MatrixBase & A_directions) { return tensorPlusOp(A, A_directions, [](Real x, Real) { return x; }); } template auto tensorPlusOp(const Eigen::MatrixBase & A, Op && oper) { Matrix A_directions; return tensorPlusOp(A, A_directions, std::forward(oper)); } template auto tensorPlus(const Eigen::MatrixBase & A) { return tensorPlusOp(A, [](Real x, Real) { return x; }); } template auto tensorSqrt(const Eigen::MatrixBase & A) { return tensorPlusOp(A, [](Real x, Int) { return std::sqrt(x); }); } } // namespace /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> MaterialAnisotropicDamage:: MaterialAnisotropicDamage(SolidMechanicsModel & model, const ID & id) : Parent(model, id), damage("damage_tensor", *this), elastic_stress("elastic_stress", *this), equivalent_strain("equivalent_strain", *this), trace_damage("trace_damage", *this), equivalent_strain_function(*this), damage_threshold_function(*this) { this->registerParam("Dc", Dc, _pat_parsable, "Critical damage"); this->damage.initialize(dim * dim); this->elastic_stress.initialize(dim * dim); this->equivalent_strain.initialize(1); this->trace_damage.initialize(1); this->trace_damage.initializeHistory(); } /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> template void MaterialAnisotropicDamage:: damageStress(Eigen::MatrixBase & sigma, const Eigen::MatrixBase & sigma_el, const Eigen::MatrixBase & D, Real TrD) { // σ_(n + 1) = (1 − D_(n + 1))^(1/2) σ~_(n + 1) (1 − D_(n + 1))^(1 / 2) // - ((1 − D_(n + 1)) : σ~_(n + 1))/ (3 - Tr(D_(n+1))) (1 − D_(n + 1)) // + 1/3 (1 - Tr(D_(n+1)) _+ + _-) I Matrix one_D = Matrix::Identity() - D; auto sqrt_one_D = tensorSqrt(one_D); Real Tr_sigma_plus; Real Tr_sigma_minus; std::tie(Tr_sigma_plus, Tr_sigma_minus) = tensorPlusTrace(sigma_el); auto I = Matrix::Identity(); sigma = sqrt_one_D * sigma_el * sqrt_one_D - (one_D.doubleDot(sigma_el) / (dim - TrD) * one_D) + 1. / dim * ((1 - TrD) * Tr_sigma_plus - Tr_sigma_minus) * I; } /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> template void MaterialAnisotropicDamage::computeStressOnQuad(Args && args) { auto & sigma = tuple::get<"sigma"_h>(args); auto & grad_u = tuple::get<"grad_u"_h>(args); auto & sigma_th = tuple::get<"sigma_th"_h>(args); auto & sigma_el = tuple::get<"sigma_el"_h>(args); auto & epsilon_hat = tuple::get<"epsilon_hat"_h>(args); auto & D = tuple::get<"damage"_h>(args); auto & TrD_n_1 = tuple::get<"TrD_n_1"_h>(args); auto & TrD = tuple::get<"TrD"_h>(args); auto & equivalent_strain_data = tuple::get<"equivalent_strain_data"_h>(args); auto & damage_threshold_data = tuple::get<"damage_threshold_data"_h>(args); Matrix Dtmp; Real TrD_n_1_tmp; Matrix epsilon; // yes you read properly this is a label for a goto auto computeDamage = [&]() { MaterialElastic::computeStressOnQuad(args); epsilon = this->template gradUToEpsilon(grad_u); // evaluate the damage criteria epsilon_hat = equivalent_strain_function(epsilon, equivalent_strain_data); // evolve the damage if needed auto K_TrD = damage_threshold_function.K(TrD, damage_threshold_data); auto f = epsilon_hat - K_TrD; // if test function > 0 evolve the damage if (f > 0) { TrD_n_1_tmp = damage_threshold_function.K_inv(epsilon_hat, damage_threshold_data); auto epsilon_plus = tensorPlus(epsilon); auto delta_lambda = (TrD_n_1_tmp - TrD) / (epsilon_hat * epsilon_hat); Dtmp = D + delta_lambda * epsilon_plus; return true; } return false; }; // compute a temporary version of the new damage auto is_damage_updated = computeDamage(); if (is_damage_updated) { /// Check and correct for broken case if (Dtmp.trace() > Dc) { if (epsilon.trace() > 0) { // tensile loading auto kpa = this->kpa; auto lambda = this->lambda; // change kappa to Kappa_broken = (1-Dc) Kappa kpa = (1 - Dc) * kpa; this->E = 9 * kpa * (kpa - lambda) / (3 * kpa - lambda); this->nu = lambda / (3 * kpa - lambda); this->updateInternalParameters(); computeDamage(); } else if (std::abs(epsilon.trace()) < 1e-10) { // deviatoric case Matrix n; std::vector ns; tensorPlusOp(Dtmp, n, [&](Real x, Int i) { if (x > this->Dc) { ns.push_back(i); return this->Dc; } return x; }); } } TrD_n_1 = TrD_n_1_tmp; D = Dtmp; } else { TrD_n_1 = TrD; } // apply the damage to the stress damageStress(sigma, sigma_el, D, TrD_n_1); } /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> void MaterialAnisotropicDamage::computeStress(ElementType type, GhostType ghost_type) { for (auto && args : getArguments(type, ghost_type)) { computeStressOnQuad(args); } } /* -------------------------------------------------------------------------- */ /* EquivalentStrain functions */ /* -------------------------------------------------------------------------- */ template class EquivalentStrainMazars : public EmptyIteratorContainer { public: EquivalentStrainMazars(Material & /*mat*/) {} template Real operator()(const Eigen::MatrixBase & epsilon, Other &&... /*other*/) { Real epsilon_hat = 0.; std::tie(epsilon_hat, std::ignore) = tensorPlusTrace(epsilon); return std::sqrt(epsilon_hat); } }; template class EquivalentStrainMazarsDruckerPrager : public EquivalentStrainMazars { public: EquivalentStrainMazarsDruckerPrager(Material & mat) : EquivalentStrainMazars(mat) { mat.registerParam("k", k, _pat_parsable, "k"); } template Real operator()(const Eigen::MatrixBase & epsilon, Real) { Real epsilon_hat = EquivalentStrainMazars::operator()(epsilon); epsilon_hat += k * epsilon.trace(); return epsilon_hat; } protected: Real k; }; /* -------------------------------------------------------------------------- */ /* DamageThreshold functions */ /* -------------------------------------------------------------------------- */ template class DamageThresholdLinear : public EmptyIteratorContainer { public: DamageThresholdLinear(Material & mat) : mat(mat) { mat.registerParam("A", A, _pat_parsable, "A"); mat.registerParam("K0", K0, _pat_parsable, "K0"); } template Real K(Real x, Other &&... /*other*/) { return 1. / A * x + K0; } template Real K_inv(Real x, Other &&... /*other*/) { return A * (x - K0); } private: Material & mat; Real A; Real K0; }; template class DamageThresholdTan : public EmptyIteratorContainer { public: DamageThresholdTan(Material & mat) : mat(mat) { mat.registerParam("a", a, _pat_parsable, "a"); mat.registerParam("A", A, _pat_parsable, "A"); mat.registerParam("K0", K0, _pat_parsable, "K0"); } template Real K(Real x, Other &&... /*other*/) { return a * std::tan(std::atan2(x, a) - std::atan2(K0, a)); } template Real K_inv(Real x, Other &&... /*other*/) { return a * A * (std::atan2(x, a) - std::atan2(K0, a)); } private: Material & mat; Real a{2.93e-4}; Real A{5e3}; Real K0{5e-5}; }; } // namespace akantu #endif /* AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_mazars.hh b/src/model/solid_mechanics/materials/material_damage/material_mazars.hh index 60e688969..f1f93f7fc 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_mazars.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_mazars.hh @@ -1,130 +1,130 @@ /** * @file material_mazars.hh * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Marion Estelle Chambart * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Jul 24 2020 * * @brief Material Following the Mazars law for damage evolution * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_MAZARS_HH_ #define AKANTU_MATERIAL_MAZARS_HH_ namespace akantu { /** * Material Mazars * * parameters in the material files : * - rho : density (default: 0) * - E : Young's modulus (default: 0) * - nu : Poisson's ratio (default: 1/2) * - K0 : Damage threshold * - At : Parameter damage traction 1 * - Bt : Parameter damage traction 2 * - Ac : Parameter damage compression 1 * - Bc : Parameter damage compression 2 * - beta : Parameter for shear */ template class Parent = MaterialElastic> class MaterialMazars : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ using parent_damage = MaterialDamage; public: MaterialMazars(SolidMechanicsModel & model, const ID & id = ""); ~MaterialMazars() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; protected: /// constitutive law for a given quadrature point - template inline void computeStressOnQuad(Args && arguments); + template inline void computeStressOnQuad(Args && args); template - inline void computeDamageAndStressOnQuad(Args && arguments); + inline void computeDamageAndStressOnQuad(Args && args); template inline void - computeDamageOnQuad(Args && arguments, + computeDamageOnQuad(Args && args, const Eigen::MatrixBase & epsilon_princ); public: decltype(auto) getArguments(ElementType el_type, GhostType ghost_type) { return zip_append( parent_damage::getArguments(el_type, ghost_type), tuple::get<"K0"_h>() = make_view(this->K0(el_type, ghost_type)), tuple::get<"Ehat"_h>() = broadcast(this->Ehat, this->damage(el_type, ghost_type).size())); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// damage threshold RandomInternalField K0; /// parameter damage traction 1 Real At; /// parameter damage traction 2 Real Bt; /// parameter damage compression 1 Real Ac; /// parameter damage compression 2 Real Bc; /// parameter for shear Real beta; /// specify the variable to average false = ehat, true = damage (only valid /// for non local version) bool damage_in_compute_stress; Real Ehat{0}; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_mazars_inline_impl.hh" #endif /* __AKANTU_MATERIAL_MAZARS_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh b/src/model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh index 468bdb510..ddac2dab8 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh @@ -1,190 +1,189 @@ /** * @file material_mazars_inline_impl.hh * * @author Marion Estelle Chambart * @author Marion Estelle Chambart * @author Nicolas Richart * * @date creation: Wed Apr 06 2011 * @date last modification: Thu Feb 20 2020 * * @brief Implementation of the inline functions of the material 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 . * */ /* -------------------------------------------------------------------------- */ #include "material_linear_isotropic_hardening.hh" #include "material_mazars.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class Parent> MaterialMazars::MaterialMazars(SolidMechanicsModel & model, const ID & id) : parent_damage(model, id), K0("K0", *this), damage_in_compute_stress(true) { this->registerParam("K0", this->K0, _pat_parsable, "K0"); this->registerParam("At", this->At, Real(0.8), _pat_parsable, "At"); this->registerParam("Ac", this->Ac, Real(1.4), _pat_parsable, "Ac"); this->registerParam("Bc", this->Bc, Real(1900.), _pat_parsable, "Bc"); this->registerParam("Bt", this->Bt, Real(12000.), _pat_parsable, "Bt"); this->registerParam("beta", this->beta, Real(1.06), _pat_parsable, "beta"); this->K0.initialize(1); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialMazars::computeStress(ElementType el_type, GhostType ghost_type) { auto & damage = this->damage(el_type, ghost_type); Real Ehat = 0; auto && arguments = getArguments(el_type, ghost_type); - for (auto && data : arguments) { - computeStressOnQuad(data); + for (auto && args : arguments) { + computeStressOnQuad(args); } } /* -------------------------------------------------------------------------- */ template class Parent> template -inline void -MaterialMazars::computeStressOnQuad(Args && arguments) { - Parent::computeStressOnQuad(arguments); +inline void MaterialMazars::computeStressOnQuad(Args && args) { + Parent::computeStressOnQuad(args); - auto && grad_u = tuple::get<"grad_u"_h>(arguments); + auto && grad_u = tuple::get<"grad_u"_h>(args); static_if(tuple::has_t<"inelastic_strain"_h, Args>()) .then([&grad_u](auto && args) { grad_u -= tuple::get<"inelastic_strain"_h>(args); - })(std::forward(arguments)); + })(std::forward(args)); Matrix epsilon = Matrix::Zero(); epsilon.block(0, 0) = Material::gradUToEpsilon(grad_u); Vector Fdiag; - epsilon.eigh(Fdiag); + epsilon.eig(Fdiag); - auto && Ehat = tuple::get<"Ehat"_h>(arguments); + auto && Ehat = tuple::get<"Ehat"_h>(args); Ehat = 0.; for (Int i = 0; i < 3; ++i) { Real epsilon_p = std::max(Real(0.), Fdiag(i)); Ehat += epsilon_p * epsilon_p; } Ehat = std::sqrt(Ehat); if (damage_in_compute_stress) { - computeDamageOnQuad(arguments, Fdiag); + computeDamageOnQuad(args, Fdiag); } if (not this->is_non_local) { - computeDamageAndStressOnQuad(arguments); + computeDamageAndStressOnQuad(args); } } /* -------------------------------------------------------------------------- */ template class Parent> template inline void -MaterialMazars::computeDamageAndStressOnQuad(Args && arguments) { - auto && grad_u = tuple::get<"grad_u"_h>(arguments); +MaterialMazars::computeDamageAndStressOnQuad(Args && args) { + auto && grad_u = tuple::get<"grad_u"_h>(args); if (not damage_in_compute_stress) { Vector Fdiag; Matrix epsilon = Matrix::Zero(); epsilon.block(0, 0) = Material::gradUToEpsilon(grad_u); - epsilon.eigh(Fdiag); - computeDamageOnQuad(arguments, Fdiag); + epsilon.eig(Fdiag); + computeDamageOnQuad(args, Fdiag); } - auto && sigma = tuple::get<"sigma"_h>(arguments); - auto && dam = tuple::get<"damage"_h>(arguments); + auto && sigma = tuple::get<"sigma"_h>(args); + auto && dam = tuple::get<"damage"_h>(args); sigma *= 1 - dam; static_if(tuple::has_t<"inelastic_strain"_h, Args>()) .then([&grad_u](auto && args) { grad_u += tuple::get<"inelastic_strain"_h>(args); - })(std::forward(arguments)); + })(std::forward(args)); } /* -------------------------------------------------------------------------- */ template class Parent> template inline void MaterialMazars::computeDamageOnQuad( - Args && arguments, const Eigen::MatrixBase & epsilon_princ) { - auto && dam = tuple::get<"damage"_h>(arguments); - auto && Ehat = tuple::get<"Ehat"_h>(arguments); + Args && args, const Eigen::MatrixBase & epsilon_princ) { + auto && dam = tuple::get<"damage"_h>(args); + auto && Ehat = tuple::get<"Ehat"_h>(args); auto Fs = Ehat - K0; if (Fs <= 0.) { return; } auto dam_t = 1 - K0 * (1 - At) / Ehat - At * (exp(-Bt * (Ehat - K0))); auto dam_c = 1 - K0 * (1 - Ac) / Ehat - Ac * (exp(-Bc * (Ehat - K0))); auto Cdiag = this->E * (1 - this->nu) / ((1 + this->nu) * (1 - 2 * this->nu)); Vector sigma_princ; sigma_princ(0) = Cdiag * epsilon_princ(0) + this->lambda * (epsilon_princ(1) + epsilon_princ(2)); sigma_princ(1) = Cdiag * epsilon_princ(1) + this->lambda * (epsilon_princ(0) + epsilon_princ(2)); sigma_princ(2) = Cdiag * epsilon_princ(2) + this->lambda * (epsilon_princ(1) + epsilon_princ(0)); Vector sigma_p; for (Int i = 0; i < 3; i++) { sigma_p(i) = std::max(Real(0.), sigma_princ(i)); } // sigma_p *= 1. - dam; auto trace_p = this->nu / this->E * (sigma_p(0) + sigma_p(1) + sigma_p(2)); Real alpha_t = 0; for (Int i = 0; i < 3; ++i) { auto epsilon_t = (1 + this->nu) / this->E * sigma_p(i) - trace_p; auto epsilon_p = std::max(Real(0.), epsilon_princ(i)); alpha_t += epsilon_t * epsilon_p; } alpha_t /= Ehat * Ehat; alpha_t = std::min(alpha_t, Real(1.)); auto alpha_c = 1. - alpha_t; alpha_t = std::pow(alpha_t, beta); alpha_c = std::pow(alpha_c, beta); auto damtemp = alpha_t * dam_t + alpha_c * dam_c; dam = std::max(damtemp, dam); dam = std::min(dam, Real(1.)); } } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh b/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh index 618a13ed8..78ef122f5 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh @@ -1,147 +1,147 @@ /** * @file material_von_mises_mazars_inline_impl.hh * * @author Mohit Pundir * * @date creation: Wed Apr 06 2011 * @date last modification: Wed Dec 23 2020 * * @brief Mazars damage with Von Misses criteria * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material_von_mises_mazars.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialVonMisesMazars::computeStressOnQuad( const Matrix & grad_u, Matrix & sigma, Real & dam, Real & Ehat) { Matrix epsilon(3, 3); epsilon.zero(); for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { epsilon(i, j) = .5 * (grad_u(i, j) + grad_u(j, i)); } } Vector Fdiag(3); - epsilon.eigh(Fdiag); + epsilon.eig(Fdiag); Ehat = 0.; for (UInt i = 0; i < 3; ++i) { Real epsilon_p = std::max(Real(0.), Fdiag(i)); Ehat += epsilon_p * epsilon_p; } Ehat = std::sqrt(Ehat); // MaterialElastic::computeStressOnQuad(grad_u, sigma); if (damage_in_compute_stress) { computeDamageOnQuad(Ehat, sigma, Fdiag, dam); } if (not this->is_non_local) { computeDamageAndStressOnQuad(grad_u, sigma, dam, Ehat); } } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialVonMisesMazars::computeDamageAndStressOnQuad( const Matrix & grad_u, Matrix & sigma, Real & dam, Real & Ehat) { if (!damage_in_compute_stress) { auto && Fdiag = Vector::Zero(); auto && epsilon = Matrix::Zero(); epsilon.block(0, 0, dim, dim) = Material::gradUToEpsilon(grad_u); - epsilon.eigh(Fdiag); + epsilon.eig(Fdiag); computeDamageOnQuad(Ehat, sigma, Fdiag, dam); } sigma *= 1 - dam; } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialVonMisesMazars::computeDamageOnQuad( const Real & epsilon_equ, __attribute__((unused)) const Matrix & sigma, const Vector & epsilon_princ, Real & dam) { Real Fs = epsilon_equ - K0; if (Fs > 0.) { Real dam_t; Real dam_c; dam_t = 1 - K0 * (1 - At) / epsilon_equ - At * (exp(-Bt * (epsilon_equ - K0))); dam_c = 1 - K0 * (1 - Ac) / epsilon_equ - Ac * (exp(-Bc * (epsilon_equ - K0))); Real Cdiag; Cdiag = this->E * (1 - this->nu) / ((1 + this->nu) * (1 - 2 * this->nu)); Vector sigma_princ(3); sigma_princ(0) = Cdiag * epsilon_princ(0) + this->lambda * (epsilon_princ(1) + epsilon_princ(2)); sigma_princ(1) = Cdiag * epsilon_princ(1) + this->lambda * (epsilon_princ(0) + epsilon_princ(2)); sigma_princ(2) = Cdiag * epsilon_princ(2) + this->lambda * (epsilon_princ(1) + epsilon_princ(0)); Vector sigma_p(3); for (UInt i = 0; i < 3; i++) { sigma_p(i) = std::max(Real(0.), sigma_princ(i)); } // sigma_p *= 1. - dam; Real trace_p = this->nu / this->E * (sigma_p(0) + sigma_p(1) + sigma_p(2)); Real alpha_t = 0; for (UInt i = 0; i < 3; ++i) { Real epsilon_t = (1 + this->nu) / this->E * sigma_p(i) - trace_p; Real epsilon_p = std::max(Real(0.), epsilon_princ(i)); alpha_t += epsilon_t * epsilon_p; } alpha_t /= epsilon_equ * epsilon_equ; alpha_t = std::min(alpha_t, Real(1.)); Real alpha_c = 1. - alpha_t; alpha_t = std::pow(alpha_t, beta); alpha_c = std::pow(alpha_c, beta); Real damtemp; damtemp = alpha_t * dam_t + alpha_c * dam_c; dam = std::max(damtemp, dam); dam = std::min(dam, Real(1.)); } } } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc index 50174ba67..4c37449bf 100644 --- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc +++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc @@ -1,269 +1,269 @@ /** * @file material_elastic_linear_anisotropic.cc * * @author Aurelia Isabel Cuba Ramos * @author Till Junge * @author Enrico Milanese * @author Nicolas Richart * * @date creation: Wed Sep 25 2013 * @date last modification: Fri Jul 24 2020 * * @brief Anisotropic elastic material * * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "material_elastic_linear_anisotropic.hh" #include "solid_mechanics_model.hh" #include #include namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialElasticLinearAnisotropic::MaterialElasticLinearAnisotropic( SolidMechanicsModel & model, const ID & id, bool symmetric) : Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim), C(voigt_h::size, voigt_h::size), eigC(voigt_h::size), symmetric(symmetric), was_stiffness_assembled(false) { AKANTU_DEBUG_IN(); this->dir_vecs.push_back(std::make_unique>()); (*this->dir_vecs.back())[0] = 1.; this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod, "Direction of main material axis"); if (dim > 1) { this->dir_vecs.push_back(std::make_unique>()); (*this->dir_vecs.back())[1] = 1.; this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod, "Direction of secondary material axis"); } if (dim > 2) { this->dir_vecs.push_back(std::make_unique>()); (*this->dir_vecs.back())[2] = 1.; this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod, "Direction of tertiary material axis"); } for (auto i : arange(voigt_h::size)) { decltype(i) start = 0; if (this->symmetric) { start = i; } for (auto j : arange(start, voigt_h::size)) { auto param = "C" + std::to_string(i + 1) + std::to_string(j + 1); this->registerParam(param, this->Cprime(i, j), Real(0.), _pat_parsmod, "Coefficient " + param); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::updateInternalParameters() { Material::updateInternalParameters(); if (this->symmetric) { for (auto i : arange(voigt_h::size)) { for (auto j : arange(i + 1, voigt_h::size)) { this->Cprime(j, i) = this->Cprime(i, j); } } } this->rotateCprime(); - this->C.eigh(this->eigC); + this->C.eig(this->eigC); this->was_stiffness_assembled = false; } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::rotateCprime() { // start by filling the empty parts fo Cprime auto diff = Dim * Dim - voigt_h::size; for (auto i : arange(voigt_h::size, Dim * Dim)) { for (auto j : arange(Dim * Dim)) { this->Cprime(i, j) = this->Cprime(i - diff, j); } } for (auto i : arange(Dim * Dim)) { for (auto j : arange(voigt_h::size, Dim * Dim)) { this->Cprime(i, j) = this->Cprime(i, j - diff); } } // construction of rotator tensor // normalise rotation matrix for (auto j : arange(Dim)) { auto && rot_vec = this->rot_mat(j); rot_vec = *this->dir_vecs[j]; rot_vec.normalize(); } // make sure the vectors form a right-handed base Vector v1, v2, v3; v3.zero(); if (Dim == 2) { for (auto i : arange(Dim)) { v1[i] = this->rot_mat(0, i); v2[i] = this->rot_mat(1, i); } v3 = v1.cross(v2); if (v3.norm() < 8 * std::numeric_limits::epsilon()) { AKANTU_ERROR("The axis vectors parallel."); } v3.normalize(); } else if (Dim == 3) { v1 = this->rot_mat(0); v2 = this->rot_mat(1); v3 = this->rot_mat(2); } auto test_axis = v1.cross(v2) - v3; if (test_axis.norm() > 8 * std::numeric_limits::epsilon()) { AKANTU_ERROR("The axis vectors do not form a right-handed coordinate " << "system. I. e., ||n1 x n2 - n3|| should be zero, but " << "it is " << test_axis.norm() << "."); } // create the rotator and the reverse rotator Matrix rotator; Matrix revrotor; for (auto i : arange(Dim)) { for (auto j : arange(Dim)) { for (auto k : arange(Dim)) { for (auto l : arange(Dim)) { auto I = voigt_h::mat[i][j]; auto J = voigt_h::mat[k][l]; rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j); revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l); } } } } // create the full rotated matrix Matrix Cfull(Dim * Dim, Dim * Dim); Cfull = rotator * Cprime * revrotor; for (auto i : arange(voigt_h::size)) { for (auto j : arange(voigt_h::size)) { this->C(i, j) = Cfull(i, j); } } } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeStress( ElementType el_type, GhostType ghost_type) { auto && arguments = getArguments(el_type, ghost_type); for (auto && data : arguments) { this->computeStressOnQuad(data); } } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeTangentModuli( ElementType el_type, Array & tangent_matrix, GhostType ghost_type) { auto && arguments = Material::getArgumentsTangent(tangent_matrix, el_type, ghost_type); for (auto && args : arguments) { this->computeTangentModuliOnQuad(args); } this->was_stiffness_assembled = true; } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computePotentialEnergy( ElementType el_type) { AKANTU_DEBUG_ASSERT(!this->finite_deformation, "finite deformation not possible in material anisotropic " "(TO BE IMPLEMENTED)"); auto && arguments = Material::getArguments(el_type, _not_ghost); for (auto && args : zip(arguments, this->potential_energy(el_type, _not_ghost))) { this->computePotentialEnergyOnQuad(std::get<0>(args), std::get<1>(args)); } } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computePotentialEnergyByElement( ElementType type, Int index, Vector & epot_on_quad_points) { auto gradu_view = make_view(this->gradu(type)); auto stress_view = make_view(this->stress(type)); auto nb_quadrature_points = this->fem.getNbIntegrationPoints(type); auto gradu_it = gradu_view.begin() + index * nb_quadrature_points; auto gradu_end = gradu_it + nb_quadrature_points; auto stress_it = stress_view.begin() + index * nb_quadrature_points; auto stress_end = stress_it + nb_quadrature_points; auto epot_quad = epot_on_quad_points.begin(); Matrix grad_u(dim, dim); for (auto data : zip(tuple::get<"grad_u"_h>() = range(gradu_it, gradu_end), tuple::get<"sigma"_h>() = range(stress_it, stress_end), tuple::get<"Epot"_h>() = epot_on_quad_points)) { this->computePotentialEnergyOnQuad(data, tuple::get<"Epot"_h>(data)); } } /* -------------------------------------------------------------------------- */ template Real MaterialElasticLinearAnisotropic::getCelerity( const Element & /*element*/) const { return std::sqrt(this->eigC(0) / rho); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(elastic_anisotropic, MaterialElasticLinearAnisotropic); } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc index a1feeff7a..eba7f2bba 100644 --- a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc +++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc @@ -1,147 +1,147 @@ /** * @file material_elastic_orthotropic.cc * * @author Till Junge * @author Enrico Milanese * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Apr 09 2021 * * @brief Orthotropic elastic material * * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "material_elastic_orthotropic.hh" #include "solid_mechanics_model.hh" #include namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialElasticOrthotropic::MaterialElasticOrthotropic( SolidMechanicsModel & model, const ID & id) : MaterialElasticLinearAnisotropic(model, id) { AKANTU_DEBUG_IN(); this->registerParam("E1", E1, Real(0.), _pat_parsmod, "Young's modulus (n1)"); this->registerParam("E2", E2, Real(0.), _pat_parsmod, "Young's modulus (n2)"); this->registerParam("nu12", nu12, Real(0.), _pat_parsmod, "Poisson's ratio (12)"); this->registerParam("G12", G12, Real(0.), _pat_parsmod, "Shear modulus (12)"); if (Dim > 2) { this->registerParam("E3", E3, Real(0.), _pat_parsmod, "Young's modulus (n3)"); this->registerParam("nu13", nu13, Real(0.), _pat_parsmod, "Poisson's ratio (13)"); this->registerParam("nu23", nu23, Real(0.), _pat_parsmod, "Poisson's ratio (23)"); this->registerParam("G13", G13, Real(0.), _pat_parsmod, "Shear modulus (13)"); this->registerParam("G23", G23, Real(0.), _pat_parsmod, "Shear modulus (23)"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::initMaterial() { AKANTU_DEBUG_IN(); MaterialElasticLinearAnisotropic::initMaterial(); AKANTU_DEBUG_ASSERT(not this->finite_deformation, "finite deformation not possible in material orthotropic " "(TO BE IMPLEMENTED)"); updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::updateInternalParameters() { this->C.zero(); this->Cprime.zero(); /* 1) construction of temporary material frame stiffness tensor------------ */ // http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 auto nu21 = nu12 * E2 / E1; auto nu31 = nu13 * E3 / E1; auto nu32 = nu23 * E3 / E2; // Full (i.e. dim^2 by dim^2) stiffness tensor in material frame if (Dim == 1) { AKANTU_ERROR("Dimensions 1 not implemented: makes no sense to have " "orthotropy for 1D"); } Real Gamma; if (Dim == 3) { Gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13); } if (Dim == 2) { Gamma = 1 / (1 - nu12 * nu21); } // Lamé's first parameters this->Cprime(0, 0) = E1 * (1 - nu23 * nu32) * Gamma; this->Cprime(1, 1) = E2 * (1 - nu13 * nu31) * Gamma; if (Dim == 3) { this->Cprime(2, 2) = E3 * (1 - nu12 * nu21) * Gamma; } // normalised poisson's ratio's this->Cprime(1, 0) = this->Cprime(0, 1) = E1 * (nu21 + nu31 * nu23) * Gamma; if (Dim == 3) { this->Cprime(2, 0) = this->Cprime(0, 2) = E1 * (nu31 + nu21 * nu32) * Gamma; this->Cprime(2, 1) = this->Cprime(1, 2) = E2 * (nu32 + nu12 * nu31) * Gamma; } // Lamé's second parameters (shear moduli) if (Dim == 3) { this->Cprime(3, 3) = G23; this->Cprime(4, 4) = G13; this->Cprime(5, 5) = G12; } else { this->Cprime(2, 2) = G12; } /* 1) rotation of C into the global frame */ this->rotateCprime(); - this->C.eigh(this->eigC); + this->C.eig(this->eigC); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(elastic_orthotropic, MaterialElasticOrthotropic); } // namespace akantu diff --git a/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function.cc b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function.cc new file mode 100644 index 000000000..44fd35754 --- /dev/null +++ b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function.cc @@ -0,0 +1,37 @@ +/** + * @file damaged_weight_function.cc + * + * @author Nicolas Richart + * + * @date creation sam jan 22 2022 + * + * @brief A Documented file. + * + * @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 . + * + */ +/* -------------------------------------------------------------------------- */ + +#include "damaged_weight_function.hh" + +namespace akantu { + +INSTANTIATE_NL_NEIGHBORHOOD(damage_wf, DamagedWeightFunction); + +} // namespace akantu diff --git a/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.hh b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.hh index 9f5825ef7..b101ee66c 100644 --- a/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.hh +++ b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.hh @@ -1,87 +1,85 @@ /** * @file damaged_weight_function_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Cyprien Wolff * * @date creation: Mon Aug 24 2015 * @date last modification: Fri Jan 15 2016 * * @brief Implementation of inline function of damaged weight function * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "damaged_weight_function.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline Real DamagedWeightFunction::operator()(Real r, const __attribute__((unused)) IntegrationPoint & q1, const IntegrationPoint & q2) { /// compute the weight UInt quad = q2.global_num; Array & dam_array = (*this->damage)(q2.type, q2.ghost_type); Real D = dam_array(quad); Real Radius_t = 0; Real Radius_init = this->R2; // if(D <= 0.5) // { // Radius_t = 2*D*Radius_init; // } // else // { // Radius_t = 2*Radius_init*(1-D); // } // Radius_t = Radius_init * (1 - D); Radius_init *= Radius_init; Radius_t *= Radius_t; if (Radius_t < Math::getTolerance()) { Radius_t = 0.001 * Radius_init; } Real expb = (2 * std::log(0.51)) / (std::log(1.0 - 0.49 * Radius_t / Radius_init)); Int expb_floor = std::floor(expb); Real b = expb_floor + expb_floor % 2; Real alpha = std::max(0., 1. - r * r / Radius_init); Real w = std::pow(alpha, b); return w; } /* -------------------------------------------------------------------------- */ inline void DamagedWeightFunction::init() { this->damage = &(this->manager.registerWeightFunctionInternal("damage")); } -INSTANTIATE_NL_NEIGHBORHOOD(damage_wf, DamagedWeightFunction); - } // namespace akantu diff --git a/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.cc b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.cc new file mode 100644 index 000000000..2a0501d9a --- /dev/null +++ b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.cc @@ -0,0 +1,37 @@ +/** + * @file remove_damaged_weight_function.cc + * + * @author Nicolas Richart + * + * @date creation sam jan 22 2022 + * + * @brief A Documented file. + * + * @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 . + * + */ +/* -------------------------------------------------------------------------- */ + +#include "remove_damaged_weight_function.hh" + +namespace akantu { + +INSTANTIATE_NL_NEIGHBORHOOD(remove_wf, RemoveDamagedWeightFunction); + +} // namespace akantu diff --git a/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.hh b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.hh index 57d98ae0d..054646c9f 100644 --- a/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.hh +++ b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.hh @@ -1,110 +1,108 @@ /** * @file remove_damaged_weight_function_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Cyprien Wolff * * @date creation: Mon Aug 24 2015 * @date last modification: Fri Apr 09 2021 * * @brief Implementation of inline function of remove damaged weight function * * * @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 . * */ /* -------------------------------------------------------------------------- */ //#include "remove_damaged_weight_function.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_REMOVE_DAMAGED_WEIGHT_FUNCTION_INLINE_IMPL_HH_ #define AKANTU_REMOVE_DAMAGED_WEIGHT_FUNCTION_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline Real RemoveDamagedWeightFunction::operator()( Real r, const __attribute__((unused)) IntegrationPoint & q1, const IntegrationPoint & q2) { /// compute the weight UInt quad = q2.global_num; if (q1 == q2) { return 1.; } Array & dam_array = (*this->damage)(q2.type, q2.ghost_type); Real D = dam_array(quad); Real w = 0.; if (D < damage_limit * (1 - Math::getTolerance())) { Real alpha = std::max(0., 1. - r * r / this->R2); w = alpha * alpha; } return w; } /* -------------------------------------------------------------------------- */ inline void RemoveDamagedWeightFunction::init() { this->damage = &(this->manager.registerWeightFunctionInternal("damage")); } /* -------------------------------------------------------------------------- */ inline Int RemoveDamagedWeightFunction::getNbData(const Array & elements, const SynchronizationTag & tag) const { if (tag == SynchronizationTag::_mnl_weight) { return this->manager.getModel().getNbIntegrationPoints(elements) * sizeof(Real); } return 0; } /* -------------------------------------------------------------------------- */ inline void RemoveDamagedWeightFunction::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag == SynchronizationTag::_mnl_weight) { DataAccessor::packElementalDataHelper( *damage, buffer, elements, true, this->manager.getModel().getFEEngine()); } } /* -------------------------------------------------------------------------- */ inline void RemoveDamagedWeightFunction::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag == SynchronizationTag::_mnl_weight) { DataAccessor::unpackElementalDataHelper( *damage, buffer, elements, true, this->manager.getModel().getFEEngine()); } } -INSTANTIATE_NL_NEIGHBORHOOD(remove_wf, RemoveDamagedWeightFunction); - } // namespace akantu #endif /* AKANTU_REMOVE_DAMAGED_WEIGHT_FUNCTION_INLINE_IMPL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc index 2ca89a228..18a2db33d 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc @@ -1,547 +1,547 @@ /** * @file fragment_manager.cc * * @author Aurelia Isabel Cuba Ramos * @author Marco Vocialta * * @date creation: Thu Jan 23 2014 * @date last modification: Mon Mar 29 2021 * * @brief Group manager to handle fragments * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "fragment_manager.hh" #include "aka_iterators.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "material_cohesive.hh" #include "mesh_iterators.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ FragmentManager::FragmentManager(SolidMechanicsModelCohesive & model, bool dump_data, const ID & id) : GroupManager(model.getMesh(), id), model(model), mass_center(0, model.getSpatialDimension(), "mass_center"), mass(0, model.getSpatialDimension(), "mass"), velocity(0, model.getSpatialDimension(), "velocity"), inertia_moments(0, model.getSpatialDimension(), "inertia_moments"), principal_directions( 0, model.getSpatialDimension() * model.getSpatialDimension(), "principal_directions"), quad_coordinates("quad_coordinates", id), mass_density("mass_density", id), nb_elements_per_fragment(0, 1, "nb_elements_per_fragment"), dump_data(dump_data) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); /// compute quadrature points' coordinates quad_coordinates.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(quad_coordinates, spatial_dimension, // spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getMesh().getNodes(), quad_coordinates); /// store mass density per quadrature point mass_density.initialize(mesh, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(mass_density, 1, spatial_dimension, // _not_ghost); storeMassDensityPerIntegrationPoint(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ class CohesiveElementFilter : public GroupManager::ClusteringFilter { public: CohesiveElementFilter(const SolidMechanicsModelCohesive & model, const Real max_damage = 1.) : model(model), is_unbroken(max_damage) {} bool operator()(const Element & el) const override { if (Mesh::getKind(el.type) == _ek_regular) { return true; } const auto & mat_indexes = model.getMaterialByElement(el.type, el.ghost_type); const auto & mat_loc_num = model.getMaterialLocalNumbering(el.type, el.ghost_type); const auto & mat = static_cast( model.getMaterial(mat_indexes(el.element))); auto el_index = mat_loc_num(el.element); auto nb_quad_per_element = model.getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(el.type, el.ghost_type); const auto & damage_array = mat.getDamage(el.type, el.ghost_type); AKANTU_DEBUG_ASSERT(Int(nb_quad_per_element * el_index) < damage_array.size(), "This quadrature point is out of range"); const auto * element_damage = damage_array.data() + nb_quad_per_element * el_index; auto unbroken_quads = std::count_if( element_damage, element_damage + nb_quad_per_element, is_unbroken); return (unbroken_quads > 0); } private: struct IsUnbrokenFunctor { IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {} bool operator()(const Real & x) const { return x < max_damage; } const Real max_damage; }; const SolidMechanicsModelCohesive & model; const IsUnbrokenFunctor is_unbroken; }; /* -------------------------------------------------------------------------- */ void FragmentManager::buildFragments(Real damage_limit) { AKANTU_DEBUG_IN(); if (mesh.isDistributed()) { auto & cohesive_synchronizer = model.getCohesiveSynchronizer(); cohesive_synchronizer.synchronize(model, SynchronizationTag::_smmc_damage); } auto & mesh_facets = mesh.getMeshFacets(); UInt spatial_dimension = model.getSpatialDimension(); std::string fragment_prefix("fragment"); /// generate fragments global_nb_fragment = createClusters(spatial_dimension, mesh_facets, fragment_prefix, CohesiveElementFilter(model, damage_limit)); nb_fragment = getNbElementGroups(spatial_dimension); fragment_index.resize(nb_fragment); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { auto name = std::get<0>(data).getName(); /// get fragment index std::string fragment_index_string = name.substr(fragment_prefix.size() + 1); std::get<1>(data) = std::stoul(fragment_index_string); } /// compute fragments' mass computeMass(); if (dump_data) { createDumpDataArray(fragment_index, "fragments", true); createDumpDataArray(mass, "fragments mass"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeMass() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// create a unit field per quadrature point, since to compute mass /// it's neccessary to integrate only density ElementTypeMapArray unit_field("unit_field", id); unit_field.initialize(model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); integrateFieldOnFragments(unit_field, mass); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeCenterOfMass() { AKANTU_DEBUG_IN(); /// integrate position multiplied by density integrateFieldOnFragments(quad_coordinates, mass_center); /// divide it by the fragments' mass Real * mass_storage = mass.data(); Real * mass_center_storage = mass_center.data(); UInt total_components = mass_center.size() * mass_center.getNbComponent(); for (UInt i = 0; i < total_components; ++i) { mass_center_storage[i] /= mass_storage[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeVelocity() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// compute velocity per quadrature point ElementTypeMapArray velocity_field("velocity_field", id); velocity_field.initialize( model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getVelocity(), velocity_field); /// integrate on fragments integrateFieldOnFragments(velocity_field, velocity); /// divide it by the fragments' mass Real * mass_storage = mass.data(); Real * velocity_storage = velocity.data(); UInt total_components = velocity.size() * velocity.getNbComponent(); for (UInt i = 0; i < total_components; ++i) { velocity_storage[i] /= mass_storage[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Given the distance @f$ \mathbf{r} @f$ between a quadrature point * and its center of mass, the moment of inertia is computed as \f[ * I_\mathrm{CM} = \mathrm{tr}(\mathbf{r}\mathbf{r}^\mathrm{T}) * \mathbf{I} - \mathbf{r}\mathbf{r}^\mathrm{T} \f] for more * information check Wikipedia * (http://en.wikipedia.org/wiki/Moment_of_inertia#Identities_for_a_skew-symmetric_matrix) * */ void FragmentManager::computeInertiaMoments() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); computeCenterOfMass(); /// compute local coordinates products with respect to the center of match ElementTypeMapArray moments_coords("moments_coords", id); moments_coords.initialize(model.getFEEngine(), _nb_component = spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); /// loop over fragments for (auto && data : zip(iterateElementGroups(), make_view(mass_center, spatial_dimension))) { const auto & el_list = std::get<0>(data).getElements(); auto & mass_center = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); auto & moments_coords_array = moments_coords(type); const auto & quad_coordinates_array = quad_coordinates(type); const auto & el_list_array = el_list(type); auto moments_begin = moments_coords_array.begin(spatial_dimension, spatial_dimension); auto quad_coordinates_begin = quad_coordinates_array.begin(spatial_dimension); Vector relative_coords(spatial_dimension); for (Int el = 0; el < el_list_array.size(); ++el) { auto global_el = el_list_array(el); /// loop over quadrature points for (Int q = 0; q < nb_quad_per_element; ++q) { auto global_q = global_el * nb_quad_per_element + q; auto && moments_matrix = moments_begin[global_q]; const auto & quad_coord_vector = quad_coordinates_begin[global_q]; /// to understand this read the documentation written just /// before this function relative_coords = quad_coord_vector; relative_coords -= mass_center; moments_matrix = relative_coords * relative_coords.transpose(); Real trace = moments_matrix.trace(); moments_matrix = -1. * moments_matrix + trace * Matrix::Identity(spatial_dimension, spatial_dimension); } } } } /// integrate moments Array integrated_moments(global_nb_fragment, spatial_dimension * spatial_dimension); integrateFieldOnFragments(moments_coords, integrated_moments); /// compute and store principal moments inertia_moments.resize(global_nb_fragment); principal_directions.resize(global_nb_fragment); auto integrated_moments_it = integrated_moments.begin(spatial_dimension, spatial_dimension); auto inertia_moments_it = inertia_moments.begin(spatial_dimension); auto principal_directions_it = principal_directions.begin(spatial_dimension, spatial_dimension); for (UInt frag = 0; frag < global_nb_fragment; ++frag, ++integrated_moments_it, ++inertia_moments_it, ++principal_directions_it) { - integrated_moments_it->eigh(*inertia_moments_it, *principal_directions_it); + integrated_moments_it->eig(*inertia_moments_it, *principal_directions_it); } if (dump_data) { createDumpDataArray(inertia_moments, "moments of inertia"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeAllData(Real damage_limit) { AKANTU_DEBUG_IN(); buildFragments(damage_limit); computeVelocity(); computeInertiaMoments(); computeNbElementsPerFragment(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::storeMassDensityPerIntegrationPoint() { AKANTU_DEBUG_IN(); auto spatial_dimension = model.getSpatialDimension(); for (auto type : mesh.elementTypes(_spatial_dimension = spatial_dimension, _element_kind = _ek_regular)) { auto & mass_density_array = mass_density(type); auto nb_element = mesh.getNbElement(type); auto nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); mass_density_array.resize(nb_element * nb_quad_per_element); const auto & mat_indexes = model.getMaterialByElement(type); auto * mass_density_it = mass_density_array.data(); /// store mass_density for each element and quadrature point for (Int el = 0; el < nb_element; ++el) { auto & mat = model.getMaterial(mat_indexes(el)); for (Int q = 0; q < nb_quad_per_element; ++q, ++mass_density_it) { *mass_density_it = mat.getRho(); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::integrateFieldOnFragments( ElementTypeMapArray & field, Array & output) { AKANTU_DEBUG_IN(); auto spatial_dimension = model.getSpatialDimension(); auto nb_component = output.getNbComponent(); /// integration part output.resize(global_nb_fragment); output.zero(); auto output_begin = output.begin(nb_component); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & el_list = std::get<0>(data).getElements(); auto fragment_index = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); const auto & density_array = mass_density(type); auto & field_array = field(type); const auto & elements = el_list(type); /// generate array to be integrated by filtering fragment's elements Array integration_array(elements.size() * nb_quad_per_element, nb_component); auto field_array_begin = make_view(field_array, nb_quad_per_element, nb_component).begin(); auto density_array_begin = make_view(density_array, nb_quad_per_element).begin(); for (auto && data : enumerate(make_view( integration_array, nb_quad_per_element, nb_component))) { auto global_el = elements(std::get<0>(data)); auto & int_array = std::get<1>(data); int_array = field_array_begin[global_el]; /// multiply field by density const auto & density_vector = density_array_begin[global_el]; for (Int i = 0; i < nb_quad_per_element; ++i) { for (Int j = 0; j < nb_component; ++j) { int_array(i, j) *= density_vector(i); } } } /// integrate the field over the fragment Array integrated_array(elements.size(), nb_component); model.getFEEngine().integrate(integration_array, integrated_array, nb_component, type, _not_ghost, elements); /// sum over all elements and store the result Vector output_tmp(output_begin[fragment_index]); output_tmp += std::accumulate(integrated_array.begin(nb_component), integrated_array.end(nb_component), Vector(nb_component)); } } /// sum output over all processors const auto & comm = mesh.getCommunicator(); comm.allReduce(output, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeNbElementsPerFragment() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); nb_elements_per_fragment.resize(global_nb_fragment); nb_elements_per_fragment.zero(); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & el_list = std::get<0>(data).getElements(); auto fragment_index = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_element = el_list(type).size(); nb_elements_per_fragment(fragment_index) += nb_element; } } /// sum values over all processors const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_elements_per_fragment, SynchronizerOperation::_sum); if (dump_data) { createDumpDataArray(nb_elements_per_fragment, "elements per fragment"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FragmentManager::createDumpDataArray(Array & data, std::string name, bool fragment_index_output) { AKANTU_DEBUG_IN(); if (data.empty()) { return; } auto & mesh_not_const = const_cast(mesh); auto && spatial_dimension = mesh.getSpatialDimension(); auto && nb_component = data.getNbComponent(); auto && data_begin = data.begin(nb_component); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & fragment = std::get<0>(data); auto fragment_idx = std::get<1>(data); /// loop over cluster types for (auto && type : fragment.elementTypes(spatial_dimension)) { /// init mesh data auto & mesh_data = mesh_not_const.getDataPointer( name, type, _not_ghost, nb_component); auto mesh_data_begin = mesh_data.begin(nb_component); /// fill mesh data for (const auto & elem : fragment.getElements(type)) { Vector md_tmp = mesh_data_begin[elem]; if (fragment_index_output) { md_tmp(0) = fragment_idx; } else { md_tmp = data_begin[fragment_idx]; } } } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/test/test_common/test_array.cc b/test/test_common/test_array.cc index 5dc7655c6..34c70e04e 100644 --- a/test/test_common/test_array.cc +++ b/test/test_common/test_array.cc @@ -1,292 +1,322 @@ /** * @file test_array.cc * * @author Nicolas Richart * * @date creation: Thu Nov 09 2017 * @date last modification: Wed Nov 18 2020 * * @brief Test the arry class * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { class NonTrivial { public: NonTrivial() = default; NonTrivial(int a) : a(a){}; bool operator==(const NonTrivial & rhs) { return a == rhs.a; } int a{0}; }; bool operator==(const int & a, const NonTrivial & rhs) { return a == rhs.a; } std::ostream & operator<<(std::ostream & stream, const NonTrivial & _this) { stream << _this.a; return stream; } /* -------------------------------------------------------------------------- */ using TestTypes = ::testing::Types; /* -------------------------------------------------------------------------- */ ::testing::AssertionResult AssertType(const char * /*a_expr*/, const char * /*b_expr*/, const std::type_info & a, const std::type_info & b) { if (std::type_index(a) == std::type_index(b)) return ::testing::AssertionSuccess(); return ::testing::AssertionFailure() << debug::demangle(a.name()) << " != " << debug::demangle(b.name()) << ") are different"; } /* -------------------------------------------------------------------------- */ template class ArrayConstructor : public ::testing::Test { protected: using type = T; void SetUp() override { type_str = debug::demangle(typeid(T).name()); } template decltype(auto) construct(P &&... params) { return std::make_unique>(std::forward

(params)...); } protected: std::string type_str; }; TYPED_TEST_SUITE(ArrayConstructor, TestTypes, ); TYPED_TEST(ArrayConstructor, ConstructDefault1) { auto array = this->construct(); EXPECT_EQ(0, array->size()); EXPECT_EQ(1, array->getNbComponent()); EXPECT_STREQ("", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault2) { auto array = this->construct(1000); EXPECT_EQ(1000, array->size()); EXPECT_EQ(1, array->getNbComponent()); EXPECT_STREQ("", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault3) { auto array = this->construct(1000, 10); EXPECT_EQ(1000, array->size()); EXPECT_EQ(10, array->getNbComponent()); EXPECT_STREQ("", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault4) { auto array = this->construct(1000, 10, "test"); EXPECT_EQ(1000, array->size()); EXPECT_EQ(10, array->getNbComponent()); EXPECT_STREQ("test", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault5) { auto array = this->construct(1000, 10, 1); EXPECT_EQ(1000, array->size()); EXPECT_EQ(10, array->getNbComponent()); EXPECT_EQ(1, array->operator()(10, 6)); EXPECT_STREQ("", array->getID().c_str()); } // TYPED_TEST(ArrayConstructor, ConstructDefault6) { // typename TestFixture::type defaultv[2] = {0, 1}; // auto array = this->construct(1000, 2, defaultv); // EXPECT_EQ(1000, array->size()); // EXPECT_EQ(2, array->getNbComponent()); // EXPECT_EQ(1, array->operator()(10, 1)); // EXPECT_EQ(0, array->operator()(603, 0)); // EXPECT_STREQ("", array->getID().c_str()); // } /* -------------------------------------------------------------------------- */ template class ArrayFixture : public ArrayConstructor { public: void SetUp() override { ArrayConstructor::SetUp(); array = this->construct(1000, 10); } void TearDown() override { array.reset(nullptr); } protected: std::unique_ptr> array; }; TYPED_TEST_SUITE(ArrayFixture, TestTypes, ); TYPED_TEST(ArrayFixture, Copy) { Array copy(*this->array); EXPECT_EQ(1000, copy.size()); EXPECT_EQ(10, copy.getNbComponent()); EXPECT_NE(this->array->data(), copy.data()); } TYPED_TEST(ArrayFixture, Set) { auto & arr = *(this->array); arr.set(12); EXPECT_EQ(12, arr(156, 5)); EXPECT_EQ(12, arr(520, 7)); EXPECT_EQ(12, arr(999, 9)); } TYPED_TEST(ArrayFixture, Resize) { auto & arr = *(this->array); auto * ptr = arr.data(); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_TRUE(arr.data() == nullptr or arr.data() == ptr); EXPECT_LE(0, arr.getAllocatedSize()); arr.resize(3000); EXPECT_EQ(3000, arr.size()); EXPECT_LE(3000, arr.getAllocatedSize()); ptr = arr.data(); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_TRUE(arr.data() == nullptr or arr.data() == ptr); EXPECT_LE(0, arr.getAllocatedSize()); } TYPED_TEST(ArrayFixture, PushBack) { auto & arr = *(this->array); auto * ptr = arr.data(); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_TRUE(arr.data() == nullptr or arr.data() == ptr); EXPECT_LE(0, arr.getAllocatedSize()); arr.resize(3000); EXPECT_EQ(3000, arr.size()); EXPECT_LE(3000, arr.getAllocatedSize()); ptr = arr.data(); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_TRUE(arr.data() == nullptr or arr.data() == ptr); EXPECT_LE(0, arr.getAllocatedSize()); } -TYPED_TEST(ArrayFixture, ViewVector) { +TYPED_TEST(ArrayFixture, ViewVectorDynamic) { auto && view = make_view(*this->array, 10); EXPECT_NO_THROW(view.begin()); { auto it = view.begin(); EXPECT_EQ(10, it->size()); EXPECT_PRED_FORMAT2(AssertType, typeid(*it), - typeid(Vector)); + typeid(VectorProxy)); EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), typeid(VectorProxy)); } } -TYPED_TEST(ArrayFixture, ViewMatrix) { +TYPED_TEST(ArrayFixture, ViewVectorStatic) { + auto && view = make_view<10>(*this->array); + EXPECT_NO_THROW(view.begin()); { - auto && view = make_view(*this->array, 2, 5); - - EXPECT_NO_THROW(view.begin()); - { - auto it = view.begin(); - EXPECT_EQ(10, it->size()); - EXPECT_EQ(2, it->size(0)); - EXPECT_EQ(5, it->size(1)); - - EXPECT_PRED_FORMAT2(AssertType, typeid(*it), - typeid(Matrix)); - EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), - typeid(MatrixProxy)); - } + auto it = view.begin(); + EXPECT_EQ(10, it->size()); + EXPECT_PRED_FORMAT2(AssertType, typeid(*it), + typeid(VectorProxy)); + EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), + typeid(VectorProxy)); + } +} + +TYPED_TEST(ArrayFixture, ViewMatrixStatic) { + auto && view = make_view(*this->array, 2, 5); + + EXPECT_NO_THROW(view.begin()); + { + auto it = view.begin(); + EXPECT_EQ(10, it->size()); + EXPECT_EQ(2, it->size(0)); + EXPECT_EQ(5, it->size(1)); + + EXPECT_PRED_FORMAT2(AssertType, typeid(*it), + typeid(MatrixProxy)); + EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), + typeid(MatrixProxy)); + } +} + +TYPED_TEST(ArrayFixture, ViewMatrixDynamic) { + auto && view = make_view<2, 5>(*this->array); + + EXPECT_NO_THROW(view.begin()); + { + auto it = view.begin(); + EXPECT_EQ(10, it->size()); + EXPECT_EQ(2, it->size(0)); + EXPECT_EQ(5, it->size(1)); + EXPECT_EQ(2, it->rows()); + EXPECT_EQ(5, it->cols()); + + EXPECT_PRED_FORMAT2(AssertType, typeid(*it), + typeid(MatrixProxy)); + EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), + typeid(MatrixProxy)); } } TYPED_TEST(ArrayFixture, ViewVectorWrong) { auto && view = make_view(*this->array, 11); EXPECT_THROW(view.begin(), debug::ArrayException); } TYPED_TEST(ArrayFixture, ViewMatrixWrong) { auto && view = make_view(*this->array, 3, 7); EXPECT_THROW(view.begin(), debug::ArrayException); } TYPED_TEST(ArrayFixture, ViewMatrixIter) { std::size_t count = 0; for (auto && mat : make_view(*this->array, 10, 10)) { EXPECT_EQ(100, mat.size()); EXPECT_EQ(10, mat.size(0)); EXPECT_EQ(10, mat.size(1)); EXPECT_PRED_FORMAT2(AssertType, typeid(mat), - typeid(Matrix)); + typeid(MatrixProxy)); ++count; } EXPECT_EQ(100, count); } TYPED_TEST(ArrayFixture, ConstViewVector) { const auto & carray = *this->array; auto && view = make_view(carray, 10); EXPECT_NO_THROW(view.begin()); { auto it = view.begin(); EXPECT_EQ(10, it->size()); EXPECT_PRED_FORMAT2(AssertType, typeid(*it), - typeid(Vector)); + typeid(VectorProxy)); EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), typeid(VectorProxy)); } } } // namespace diff --git a/test/test_common/test_tensors.cc b/test/test_common/test_tensors.cc index cabe39be8..7d0d7cf6f 100644 --- a/test/test_common/test_tensors.cc +++ b/test/test_common/test_tensors.cc @@ -1,555 +1,578 @@ /** * @file test_tensors.cc * * @author Nicolas Richart * * @date creation: Tue Nov 14 2017 * @date last modification: Tue Feb 05 2019 * * @brief test the tensors types * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_iterators.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { /* -------------------------------------------------------------------------- */ class TensorConstructorFixture : public ::testing::Test { public: void SetUp() override { for (auto & r : reference) { r = rand(); // google-test seeds srand() } } void TearDown() override {} template void compareToRef(const V & v) { for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i], v.data()[i]); } } protected: const int size_{24}; const std::array mat_size{{4, 6}}; // const std::array tens3_size{{4, 2, 3}}; std::array reference; }; /* -------------------------------------------------------------------------- */ class TensorFixture : public TensorConstructorFixture { public: TensorFixture() : vref(reference.data(), size_), mref(reference.data(), mat_size[0], mat_size[1]) {} protected: VectorProxy vref; MatrixProxy mref; }; /* -------------------------------------------------------------------------- */ // Vector ---------------------------------------------------------------------- TEST_F(TensorConstructorFixture, VectorDefaultConstruct) { Vector v; EXPECT_EQ(0, v.size()); EXPECT_EQ(nullptr, v.data()); } TEST_F(TensorConstructorFixture, VectorConstruct1) { Vector v(size_); EXPECT_EQ(size_, v.size()); } TEST_F(TensorConstructorFixture, VectorConstructWrapped) { VectorProxy v(reference.data(), size_); EXPECT_EQ(size_, v.size()); EXPECT_EQ(v.data(), reference.data()); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i], v(i)); EXPECT_DOUBLE_EQ(reference[i], v[i]); } } TEST_F(TensorConstructorFixture, VectorConstructInitializer) { Vector v{0., 1., 2., 3., 4., 5.}; EXPECT_EQ(6, v.size()); for (int i = 0; i < 6; ++i) { EXPECT_DOUBLE_EQ(i, v(i)); } } TEST_F(TensorConstructorFixture, VectorConstructCopy1) { VectorProxy vref(reference.data(), reference.size()); Vector v(vref); EXPECT_EQ(size_, v.size()); compareToRef(v); } TEST_F(TensorConstructorFixture, VectorConstructProxy1) { VectorProxy vref(reference.data(), reference.size()); EXPECT_EQ(size_, vref.size()); compareToRef(vref); Vector v(vref); EXPECT_EQ(size_, v.size()); compareToRef(v); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorEqual) { Vector v; v = vref; compareToRef(v); EXPECT_EQ(size_, v.size()); } TEST_F(TensorFixture, VectorEqualProxy) { VectorProxy vref_proxy(vref); Vector v; v = vref; compareToRef(v); EXPECT_EQ(size_, v.size()); } TEST_F(TensorFixture, VectorEqualProxy2) { Vector v_store(size_); v_store.zero(); VectorProxy v(v_store.data(), size_); v = vref; compareToRef(v); compareToRef(v_store); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorSet) { Vector v(vref); compareToRef(v); double r = rand(); v.set(r); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(r, v[i]); } TEST_F(TensorFixture, VectorClear) { Vector v(vref); compareToRef(v); v.zero(); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0, v[i]); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorDivide) { Vector v; double r = rand(); v = vref / r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] / r, v[i]); } } TEST_F(TensorFixture, VectorMultiply1) { Vector v; double r = rand(); v = vref * r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } } TEST_F(TensorFixture, VectorMultiply2) { Vector v; double r = rand(); v = r * vref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } } TEST_F(TensorFixture, VectorAddition) { Vector v; v = vref + vref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * 2., v[i]); } } TEST_F(TensorFixture, VectorSubstract) { Vector v; v = vref - vref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(0., v[i]); } } TEST_F(TensorFixture, VectorDivideEqual) { Vector v(vref); double r = rand(); v /= r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] / r, v[i]); } } TEST_F(TensorFixture, VectorMultiplyEqual1) { Vector v(vref); double r = rand(); v *= r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } } TEST_F(TensorFixture, VectorMultiplyEqual2) { Vector v(vref); - v *= v; + v.array() *= v.array(); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * reference[i], v[i]); } } TEST_F(TensorFixture, VectorAdditionEqual) { Vector v(vref); v += vref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * 2., v[i]); } } TEST_F(TensorFixture, VectorSubstractEqual) { Vector v(vref); v -= vref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(0., v[i]); } } /* -------------------------------------------------------------------------- */ // Matrix ---------------------------------------------------------------------- TEST_F(TensorConstructorFixture, MatrixDefaultConstruct) { Matrix m; EXPECT_EQ(0, m.size()); EXPECT_EQ(0, m.rows()); EXPECT_EQ(0, m.cols()); EXPECT_EQ(nullptr, m.data()); } TEST_F(TensorConstructorFixture, MatrixConstruct1) { double r = rand(); Matrix m(mat_size[0], mat_size[1]); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); } TEST_F(TensorConstructorFixture, MatrixConstructWrapped) { MatrixProxy m(reference.data(), mat_size[0], mat_size[1]); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); for (int i = 0; i < mat_size[0]; ++i) { for (int j = 0; j < mat_size[1]; ++j) { EXPECT_DOUBLE_EQ(reference[i + j * mat_size[0]], m(i, j)); } } compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructInitializer) { Matrix m{{0., 1., 2.}, {3., 4., 5.}}; EXPECT_EQ(6, m.size()); EXPECT_EQ(2, m.rows()); EXPECT_EQ(3, m.cols()); int c = 0; for (int i = 0; i < 2; ++i) { for (int j = 0; j < 3; ++j, ++c) { EXPECT_DOUBLE_EQ(c, m(i, j)); } } } TEST_F(TensorConstructorFixture, MatrixConstructCopy1) { MatrixProxy mref(reference.data(), mat_size[0], mat_size[1]); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructProxy1) { MatrixProxy mref(reference.data(), mat_size[0], mat_size[1]); EXPECT_EQ(size_, mref.size()); EXPECT_EQ(mat_size[0], mref.size(0)); EXPECT_EQ(mat_size[1], mref.size(1)); compareToRef(mref); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); compareToRef(m); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixEqual) { Matrix m; m = mref; compareToRef(m); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); } TEST_F(TensorFixture, MatrixEqualProxy1) { MatrixProxy mref_proxy(mref.data(), mref.rows(), mref.cols()); Matrix m; m = mref; compareToRef(m); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); } TEST_F(TensorFixture, MatrixEqualProxy2) { Matrix m_store(mat_size[0], mat_size[1]); m_store.zero(); MatrixProxy m(m_store.data(), mat_size[0], mat_size[1]); m = mref; compareToRef(m); compareToRef(m_store); } TEST_F(TensorFixture, MatrixEqualSlice) { Matrix m(mat_size[0], mat_size[1]); m.zero(); for (unsigned int i = 0; i < m.cols(); ++i) { m(i) = Vector(mref(i)); } compareToRef(m); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixSet) { Matrix m(mref); compareToRef(m); double r = rand(); m.set(r); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(r, m.array()(i)); } } TEST_F(TensorFixture, MatrixClear) { Matrix m(mref); compareToRef(m); m.zero(); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(0, m.array()(i)); } } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixDivide) { Matrix m; double r = rand(); m = mref / r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] / r, m.array()(i)); } } TEST_F(TensorFixture, MatrixMultiply1) { Matrix m; double r = rand(); m = mref * r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * r, m.array()(i)); } } TEST_F(TensorFixture, MatrixMultiply2) { Matrix m; double r = rand(); m = r * mref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * r, m.array()(i)); } } TEST_F(TensorFixture, MatrixAddition) { Matrix m; m = mref + mref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * 2., m.array()(i)); } } TEST_F(TensorFixture, MatrixSubstract) { Matrix m; m = mref - mref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(0., m.array()(i)); } } TEST_F(TensorFixture, MatrixDivideEqual) { Matrix m(mref); double r = rand(); m /= r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] / r, m.array()(i)); } } TEST_F(TensorFixture, MatrixMultiplyEqual1) { Matrix m(mref); double r = rand(); m *= r; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * r, m.array()(i)); } } TEST_F(TensorFixture, MatrixAdditionEqual) { Matrix m(mref); m += mref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i] * 2., m.array()(i)); } } TEST_F(TensorFixture, MatrixSubstractEqual) { Matrix m(mref); m -= mref; for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(0., m.array()(i)); } } TEST_F(TensorFixture, MatrixIterator) { Matrix m(mref); UInt col_count = 0; for (auto && col : m) { VectorProxy col_hand(m.data() + col_count * m.rows(), m.rows()); Vector col_wrap(col); auto comp = (col_wrap - col_hand).lpNorm(); EXPECT_DOUBLE_EQ(0., comp); ++col_count; } } TEST_F(TensorFixture, MatrixIteratorZip) { Matrix m1(mref); Matrix m2(mref); UInt col_count = 0; for (auto && col : zip(m1, m2)) { Vector col1(std::get<0>(col)); Vector col2(std::get<1>(col)); auto comp = (col1 - col2).lpNorm(); EXPECT_DOUBLE_EQ(0., comp); ++col_count; } } #if defined(AKANTU_USE_LAPACK) TEST_F(TensorFixture, MatrixEigs) { - Matrix m{{0, 1, 0, 0}, {1., 0, 0, 0}, {0, 1, 0, 1}, {0, 0, 4, 0}}; + Matrix A{ + {0, 1., 0, 0}, {1., 0, 0, 0}, {0, 1., 0, 1.}, {0, 0, 4., 0}}; - Matrix eig_vects(4, 4); - Vector eigs(4); - m.eigh(eigs, eig_vects); + Matrix v; + Vector lambda; + A.eig(lambda, v); Vector eigs_ref{2, 1., -1., -2}; - auto lambda_v = m * eig_vects; + + auto Av = (A * v).eval(); + + // Eigen::PermutationMatrix + // perm(lambda.size()); perm.setIdentity(); + + // std::sort(perm.indices().data(), + // perm.indices().data() + perm.indices().size(), + // [&lambda](const Eigen::Index & a, const Eigen::Index & b) { + // return (lambda(a) - lambda(b)) > 0; + // }); + + // std::cout << v << std::endl; + // std::cout << lambda << std::endl; + + // std::cout << v * perm << std::endl; + // std::cout << perm.transpose() * lambda << std::endl; + + // std::cout << (Av(0) - lambda(0) * v(0)).eval() << std::endl; for (int i = 0; i < 4; ++i) { - EXPECT_NEAR(eigs_ref(i), eigs(i), 1e-14); - for (int j = 0; j < 4; ++j) { - EXPECT_NEAR(lambda_v(i, j), eigs(i) * eig_vects(i, j), 1e-14); - } + EXPECT_NEAR(eigs_ref(i), lambda(i), 1e-14); + } + + for (int i = 0; i < 4; ++i) { + auto lambda_v_minus_a_v = + (lambda(i) * v(i) - Av(i)).template lpNorm(); + + EXPECT_NEAR(lambda_v_minus_a_v, 0., 1e-14); } } #endif /* -------------------------------------------------------------------------- */ } // namespace diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc index b2f4701c0..8c3700afc 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc @@ -1,912 +1,912 @@ /** * @file test_elastic_materials.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author Enrico Milanese * * @date creation: Fri Nov 17 2017 * @date last modification: Wed Nov 18 2020 * * @brief Tests the Elastic materials * * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "test_gtest_utils.hh" #include "test_material_fixtures.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; using mat_types = ::testing::Types, Traits, Traits, Traits, Traits, Traits, Traits>; /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::setParams() { Real E = 3.; Real rho = 2; setParam("E", E); setParam("rho", rho); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Matrix eps{{2}}; Matrix sigma(1, 1); Real sigma_th = 2; this->computeStressOnQuad(make_named_tuple( tuple::get<"grad_u"_h>() = eps, tuple::get<"sigma"_h>() = sigma, tuple::get<"sigma"_h>() = sigma_th)); auto solution = E * eps(0, 0) + sigma_th; EXPECT_NEAR(sigma(0, 0), solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Real eps = 2, sigma = 2; Real epot = 0; this->computePotentialEnergyOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = Matrix{{eps}}, tuple::get<"sigma"_h>() = Matrix{{sigma}}), epot); Real solution = 2; EXPECT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { Matrix tangent(1, 1); this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); EXPECT_NEAR(tangent(0, 0), E, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testCelerity() { auto wave_speed = this->getCelerity(Element()); auto solution = std::sqrt(E / rho); EXPECT_NEAR(wave_speed, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::setParams() { Real E = 1.; Real nu = .3; Real rho = 2; setParam("E", E); setParam("nu", nu); setParam("rho", rho); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Real bulk_modulus_K = E / (3 * (1 - 2 * nu)); Real shear_modulus_mu = E / (2 * (1 + nu)); Matrix rotation_matrix = getRandomRotation(); auto grad_u = this->getComposedStrain(1.).block(0, 0, 2, 2); auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix); Matrix sigma_rot(2, 2); this->computeStressOnQuad(make_named_tuple( tuple::get<"grad_u"_h>() = grad_u_rot, tuple::get<"sigma"_h>() = sigma_rot, tuple::get<"sigma"_h>() = sigma_th)); auto sigma = this->reverseRotation(sigma_rot, rotation_matrix); auto identity = Matrix::Identity(2, 2); auto strain = 0.5 * (grad_u + grad_u.transpose()); auto deviatoric_strain = strain - 1. / 3. * strain.trace() * identity; auto sigma_expected = 2 * shear_modulus_mu * deviatoric_strain + (sigma_th + 2. * bulk_modulus_K) * identity; auto diff = sigma - sigma_expected; Real stress_error = diff.lpNorm() / sigma_expected.lpNorm(); EXPECT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma{{1, 2}, {2, 4}}; Matrix eps{{1, 0}, {0, 1}}; Real epot = 0; Real solution = 2.5; this->computePotentialEnergyOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = eps, tuple::get<"sigma"_h>() = sigma), epot); EXPECT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { Matrix tangent; /* Plane Strain */ // clang-format off Matrix solution{ {1 - nu, nu, 0}, {nu, 1 - nu, 0}, {0, 0, (1 - 2 * nu) / 2}, }; // clang-format on solution *= E / ((1 + nu) * (1 - 2 * nu)); this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); Real tangent_error = (tangent - solution).norm(); EXPECT_NEAR(tangent_error, 0, 1e-14); /* Plane Stress */ this->plane_stress = true; this->updateInternalParameters(); // clang-format off solution = Matrix{ {1, nu, 0}, {nu, 1, 0}, {0, 0, (1 - nu) / 2}, }; // clang-format on solution *= E / (1 - nu * nu); this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); tangent_error = (tangent - solution).norm(); EXPECT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testCelerity() { auto push_wave_speed = this->getPushWaveSpeed(Element()); auto celerity = this->getCelerity(Element()); Real K = E / (3 * (1 - 2 * nu)); Real mu = E / (2 * (1 + nu)); Real sol = std::sqrt((K + 4. / 3 * mu) / rho); EXPECT_NEAR(push_wave_speed, sol, 1e-14); EXPECT_NEAR(celerity, sol, 1e-14); auto shear_wave_speed = this->getShearWaveSpeed(Element()); sol = std::sqrt(mu / rho); EXPECT_NEAR(shear_wave_speed, sol, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::setParams() { Real E = 1.; Real nu = .3; Real rho = 2; setParam("E", E); setParam("nu", nu); setParam("rho", rho); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Real bulk_modulus_K = E / 3. / (1 - 2. * nu); Real shear_modulus_mu = 0.5 * E / (1 + nu); Matrix rotation_matrix = getRandomRotation(); auto grad_u = this->getComposedStrain(1.); auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix); Matrix sigma_rot(3, 3); this->computeStressOnQuad(make_named_tuple( tuple::get<"grad_u"_h>() = grad_u_rot, tuple::get<"sigma"_h>() = sigma_rot, tuple::get<"sigma"_h>() = sigma_th)); auto sigma = this->reverseRotation(sigma_rot, rotation_matrix); Matrix identity; identity.eye(); Matrix strain = 0.5 * (grad_u + grad_u.transpose()); Matrix deviatoric_strain = strain - 1. / 3. * strain.trace() * identity; Matrix sigma_expected = 2 * shear_modulus_mu * deviatoric_strain + (sigma_th + 3. * bulk_modulus_K) * identity; auto diff = sigma - sigma_expected; Real stress_error = diff.lpNorm(); EXPECT_NEAR(stress_error, 0., 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma{{1, 2, 3}, {2, 4, 5}, {3, 5, 6}}; Matrix eps{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; Real epot = 0; Real solution = 5.5; this->computePotentialEnergyOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = eps, tuple::get<"sigma"_h>() = sigma), epot); EXPECT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { Matrix tangent(6, 6); // clang-format off Matrix solution{ {1 - nu, nu, nu, 0, 0, 0}, {nu, 1 - nu, nu, 0, 0, 0}, {nu, nu, 1 - nu, 0, 0, 0}, {0, 0, 0, (1 - 2 * nu) / 2, 0, 0}, {0, 0, 0, 0, (1 - 2 * nu) / 2, 0}, {0, 0, 0, 0, 0, (1 - 2 * nu) / 2}, }; // clang-format on solution *= E / ((1 + nu) * (1 - 2 * nu)); this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); Real tangent_error = (tangent - solution).norm(); EXPECT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testCelerity() { auto push_wave_speed = this->getPushWaveSpeed(Element()); auto celerity = this->getCelerity(Element()); Real K = E / (3 * (1 - 2 * nu)); Real mu = E / (2 * (1 + nu)); Real sol = std::sqrt((K + 4. / 3 * mu) / rho); EXPECT_NEAR(push_wave_speed, sol, 1e-14); EXPECT_NEAR(celerity, sol, 1e-14); auto shear_wave_speed = this->getShearWaveSpeed(Element()); sol = std::sqrt(mu / rho); EXPECT_NEAR(shear_wave_speed, sol, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::setParams() { // Note: for this test material and canonical basis coincide Vector n1{1, 0}; Vector n2{0, 1}; Real E1 = 1.; Real E2 = 2.; Real nu12 = 0.1; Real G12 = 2.; Real rho = 2.5; *this->dir_vecs[0] = n1; *this->dir_vecs[1] = n2; this->E1 = E1; this->E2 = E2; this->nu12 = nu12; this->G12 = G12; this->rho = rho; } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { UInt Dim = 2; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticOrthotropic class (it is simply given by the // columns of the rotation_matrix; the lines give the material basis expressed // in the canonical frame of reference) *this->dir_vecs[0] = rotation_matrix(0); *this->dir_vecs[1] = rotation_matrix(1); // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); // gradient in material frame of reference auto grad_u = this->getComposedStrain(2.).block<2, 2>(0, 0); // gradient in canonical basis (we need to rotate *back* to the canonical // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); // stress in the canonical basis Matrix sigma_rot(2, 2); this->computeStressOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = grad_u_rot, tuple::get<"sigma"_h>() = sigma_rot)); // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); // construction of Cijkl engineering tensor in the *material* frame of // reference // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 Real nu21 = nu12 * E2 / E1; Real gamma = 1 / (1 - nu12 * nu21); Matrix C_expected(2 * Dim, 2 * Dim); C_expected.zero(); C_expected(0, 0) = gamma * E1; C_expected(1, 1) = gamma * E2; C_expected(2, 2) = G12; C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21; // epsilon is computed directly in the *material* frame of reference Matrix epsilon = 0.5 * (grad_u + grad_u.transpose()); // sigma_expected is computed directly in the *material* frame of reference Matrix sigma_expected(Dim, Dim); for (UInt i = 0; i < Dim; ++i) { for (UInt j = 0; j < Dim; ++j) { sigma_expected(i, i) += C_expected(i, j) * epsilon(j, j); } } sigma_expected(0, 1) = sigma_expected(1, 0) = C_expected(2, 2) * 2 * epsilon(0, 1); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.lpNorm(); EXPECT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma{{1, 2}, {2, 4}}; Matrix eps{{1, 0}, {0, 1}}; Real epot = 0; Real solution = 2.5; this->computePotentialEnergyOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = eps, tuple::get<"sigma"_h>() = sigma), epot); EXPECT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { // construction of Cijkl engineering tensor in the *material* frame of // reference // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 Real nu21 = nu12 * E2 / E1; Real gamma = 1 / (1 - nu12 * nu21); Matrix C_expected(3, 3); C_expected(0, 0) = gamma * E1; C_expected(1, 1) = gamma * E2; C_expected(2, 2) = G12; C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21; Matrix tangent(3, 3); this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); Real tangent_error = (tangent - C_expected).norm(); EXPECT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testCelerity() { // construction of Cijkl engineering tensor in the *material* frame of // reference // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 Real nu21 = nu12 * E2 / E1; Real gamma = 1 / (1 - nu12 * nu21); Matrix C_expected(3, 3); C_expected(0, 0) = gamma * E1; C_expected(1, 1) = gamma * E2; C_expected(2, 2) = G12; C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21; Vector eig_expected(3); - C_expected.eigh(eig_expected); + C_expected.eig(eig_expected); auto celerity_expected = std::sqrt(eig_expected(0) / rho); auto celerity = this->getCelerity(Element()); EXPECT_NEAR(celerity_expected, celerity, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::setParams() { Vector n1{1, 0, 0}; Vector n2{0, 1, 0}; Vector n3{0, 0, 1}; Real E1 = 1.; Real E2 = 2.; Real E3 = 3.; Real nu12 = 0.1; Real nu13 = 0.2; Real nu23 = 0.3; Real G12 = 2.; Real G13 = 3.; Real G23 = 1.; Real rho = 2.3; *this->dir_vecs[0] = n1; *this->dir_vecs[1] = n2; *this->dir_vecs[2] = n3; this->E1 = E1; this->E2 = E2; this->E3 = E3; this->nu12 = nu12; this->nu13 = nu13; this->nu23 = nu23; this->G12 = G12; this->G13 = G13; this->G23 = G23; this->rho = rho; } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { UInt Dim = 3; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticOrthotropic class (it is simply given by the // columns of the rotation_matrix; the lines give the material basis expressed // in the canonical frame of reference) *this->dir_vecs[0] = rotation_matrix(0); *this->dir_vecs[1] = rotation_matrix(1); *this->dir_vecs[2] = rotation_matrix(2); // set internal Cijkl matrix expressed in the canonical frame of reference this->updateInternalParameters(); // gradient in material frame of reference auto grad_u = this->getComposedStrain(2.); // gradient in canonical basis (we need to rotate *back* to the canonical // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); // stress in the canonical basis Matrix sigma_rot(3, 3); this->computeStressOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = grad_u_rot, tuple::get<"sigma"_h>() = sigma_rot)); // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); // construction of Cijkl engineering tensor in the *material* frame of // reference // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 Real nu21 = nu12 * E2 / E1; Real nu31 = nu13 * E3 / E1; Real nu32 = nu23 * E3 / E2; Real gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13); Matrix C_expected(6, 6); C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32); C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31); C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21); C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23); C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32); C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31); C_expected(3, 3) = G23; C_expected(4, 4) = G13; C_expected(5, 5) = G12; // epsilon is computed directly in the *material* frame of reference Matrix epsilon = 0.5 * (grad_u + grad_u.transpose()); // sigma_expected is computed directly in the *material* frame of reference Matrix sigma_expected(Dim, Dim); for (UInt i = 0; i < Dim; ++i) { for (UInt j = 0; j < Dim; ++j) { sigma_expected(i, i) += C_expected(i, j) * epsilon(j, j); } } sigma_expected(0, 1) = C_expected(5, 5) * 2 * epsilon(0, 1); sigma_expected(0, 2) = C_expected(4, 4) * 2 * epsilon(0, 2); sigma_expected(1, 2) = C_expected(3, 3) * 2 * epsilon(1, 2); sigma_expected(1, 0) = sigma_expected(0, 1); sigma_expected(2, 0) = sigma_expected(0, 2); sigma_expected(2, 1) = sigma_expected(1, 2); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.lpNorm(); EXPECT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma{{1, 2, 3}, {2, 4, 5}, {3, 5, 6}}; Matrix eps{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; Real epot = 0; Real solution = 5.5; this->computePotentialEnergyOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = eps, tuple::get<"sigma"_h>() = sigma), epot); EXPECT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeTangentModuli() { // Note: for this test material and canonical basis coincide UInt Dim = 3; // construction of Cijkl engineering tensor in the *material* frame of // reference // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 auto nu21 = nu12 * E2 / E1; auto nu31 = nu13 * E3 / E1; auto nu32 = nu23 * E3 / E2; auto gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13); Matrix C_expected(2 * Dim, 2 * Dim); C_expected.zero(); C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32); C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31); C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21); C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23); C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32); C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31); C_expected(3, 3) = G23; C_expected(4, 4) = G13; C_expected(5, 5) = G12; Matrix tangent(6, 6); this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); Real tangent_error = (tangent - C_expected).norm(); EXPECT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testCelerity() { // Note: for this test material and canonical basis coincide UInt Dim = 3; // construction of Cijkl engineering tensor in the *material* frame of // reference // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 auto nu21 = nu12 * E2 / E1; auto nu31 = nu13 * E3 / E1; auto nu32 = nu23 * E3 / E2; auto gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13); Matrix C_expected(2 * Dim, 2 * Dim); C_expected.zero(); C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32); C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31); C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21); C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23); C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32); C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31); C_expected(3, 3) = G23; C_expected(4, 4) = G13; C_expected(5, 5) = G12; Vector eig_expected(6); - C_expected.eigh(eig_expected); + C_expected.eig(eig_expected); auto celerity_expected = std::sqrt(eig_expected(0) / rho); auto celerity = this->getCelerity(Element()); EXPECT_NEAR(celerity_expected, celerity, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::setParams() { Matrix C{{1.0, 0.3, 0.4}, {0.3, 2.0, 0.1}, {0.4, 0.1, 1.5}}; Cprime.block<3, 3>(0, 0) = C; this->rho = 2.7; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticLinearAnisotropic class (it is simply given by // the columns of the rotation_matrix; the lines give the material basis // expressed in the canonical frame of reference) *this->dir_vecs[0] = rotation_matrix(0); *this->dir_vecs[1] = rotation_matrix(1); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Matrix C{{1.0, 0.3, 0.4}, {0.3, 2.0, 0.1}, {0.4, 0.1, 1.5}}; Matrix rotation_matrix; rotation_matrix(0) = *this->dir_vecs[0]; rotation_matrix(1) = *this->dir_vecs[1]; // gradient in material frame of reference auto grad_u = this->getComposedStrain(1.).block<2, 2>(0, 0); // gradient in canonical basis (we need to rotate *back* to the canonical // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); // stress in the canonical basis Matrix sigma_rot; this->computeStressOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = grad_u_rot, tuple::get<"sigma"_h>() = sigma_rot)); // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); // epsilon is computed directly in the *material* frame of reference Matrix epsilon = (grad_u + grad_u.transpose()) / 2.; Vector epsilon_voigt = VoigtHelper<2>::matrixToVoigtWithFactors(epsilon); // sigma_expected is computed directly in the *material* frame of reference Vector sigma_voigt = C * epsilon_voigt; Matrix sigma_expected = VoigtHelper<2>::voigtToMatrix(sigma_voigt); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.lpNorm(); EXPECT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma{{1, 2}, {2, 4}}; Matrix eps{{1, 0}, {0, 1}}; Real epot = 0; Real solution = 2.5; this->computePotentialEnergyOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = eps, tuple::get<"sigma"_h>() = sigma), epot); EXPECT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial< MaterialElasticLinearAnisotropic<2>>::testComputeTangentModuli() { Matrix tangent(3, 3); this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); Real tangent_error = (tangent - C).norm(); EXPECT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testCelerity() { Vector eig_expected(3); - C.eigh(eig_expected); + C.eig(eig_expected); auto celerity_expected = std::sqrt(eig_expected(0) / this->rho); auto celerity = this->getCelerity(Element()); EXPECT_NEAR(celerity_expected, celerity, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::setParams() { // Note: for this test material and canonical basis coincide Matrix C{ {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, {0.3, 2.0, 0.1, 0.2, 0.3, 0.2}, {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, {0.3, 0.2, 0.1, 2.4, 0.1, 0.4}, {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, {0.1, 0.2, 0.3, 0.4, 0.1, 1.2}}; for (auto i = 0u; i < C.rows(); ++i) for (auto j = 0u; j < C.cols(); ++j) this->Cprime(i, j) = C(i, j); this->rho = 2.9; // material frame of reference is rotate by rotation_matrix starting from // canonical basis Matrix rotation_matrix = getRandomRotation(); // canonical basis as expressed in the material frame of reference, as // required by MaterialElasticLinearAnisotropic class (it is simply given by // the columns of the rotation_matrix; the lines give the material basis // expressed in the canonical frame of reference) *this->dir_vecs[0] = rotation_matrix(0); *this->dir_vecs[1] = rotation_matrix(1); *this->dir_vecs[2] = rotation_matrix(2); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Matrix C{ {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, {0.3, 2.0, 0.1, 0.2, 0.3, 0.2}, {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, {0.3, 0.2, 0.1, 2.4, 0.1, 0.4}, {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, {0.1, 0.2, 0.3, 0.4, 0.1, 1.2}}; Matrix rotation_matrix; rotation_matrix(0) = *this->dir_vecs[0]; rotation_matrix(1) = *this->dir_vecs[1]; rotation_matrix(2) = *this->dir_vecs[2]; // gradient in material frame of reference auto grad_u = this->getComposedStrain(2.); // gradient in canonical basis (we need to rotate *back* to the canonical // basis) auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix); // stress in the canonical basis Matrix sigma_rot; this->computeStressOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = grad_u_rot, tuple::get<"sigma"_h>() = sigma_rot)); // stress in the material reference (we need to apply the rotation) auto sigma = this->applyRotation(sigma_rot, rotation_matrix); // epsilon is computed directly in the *material* frame of reference Matrix epsilon = (grad_u + grad_u.transpose()) / 2.; Vector epsilon_voigt = VoigtHelper<3>::matrixToVoigtWithFactors(epsilon); // sigma_expected is computed directly in the *material* frame of reference Vector sigma_voigt = C * epsilon_voigt; Matrix sigma_expected = VoigtHelper<3>::voigtToMatrix(sigma_voigt); // sigmas are checked in the *material* frame of reference auto diff = sigma - sigma_expected; Real stress_error = diff.lpNorm(); EXPECT_NEAR(stress_error, 0., 1e-13); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testEnergyDensity() { Matrix sigma{{1, 2, 3}, {2, 4, 5}, {3, 5, 6}}; Matrix eps{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; Real epot = 0; Real solution = 5.5; this->computePotentialEnergyOnQuad( make_named_tuple(tuple::get<"grad_u"_h>() = eps, tuple::get<"sigma"_h>() = sigma), epot); EXPECT_NEAR(epot, solution, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial< MaterialElasticLinearAnisotropic<3>>::testComputeTangentModuli() { Matrix tangent; this->computeTangentModuliOnQuad( make_named_tuple(tuple::get<"tangent_moduli"_h>() = tangent)); Real tangent_error = (tangent - C).norm(); EXPECT_NEAR(tangent_error, 0, 1e-14); } /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testCelerity() { Vector eig_expected; - C.eigh(eig_expected); + C.eig(eig_expected); auto celerity_expected = std::sqrt(eig_expected(0) / this->rho); auto celerity = this->getCelerity(Element()); EXPECT_NEAR(celerity_expected, celerity, 1e-14); } /* -------------------------------------------------------------------------- */ namespace { template class TestElasticMaterialFixture : public ::TestMaterialFixture {}; TYPED_TEST_SUITE(TestElasticMaterialFixture, mat_types, ); TYPED_TEST(TestElasticMaterialFixture, ComputeStress) { this->material->testComputeStress(); } TYPED_TEST(TestElasticMaterialFixture, EnergyDensity) { this->material->testEnergyDensity(); } TYPED_TEST(TestElasticMaterialFixture, ComputeTangentModuli) { this->material->testComputeTangentModuli(); } TYPED_TEST(TestElasticMaterialFixture, ComputeCelerity) { this->material->testCelerity(); } } // namespace