diff --git a/language_bindings/python/CMakeLists.txt b/language_bindings/python/CMakeLists.txt index b074634..4820888 100644 --- a/language_bindings/python/CMakeLists.txt +++ b/language_bindings/python/CMakeLists.txt @@ -1,93 +1,94 @@ #============================================================================== # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief configuration for python binding using pybind11 # # @section LICENSE # # Copyright © 2018 Till Junge # # µSpectre 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, or (at # your option) any later version. # # µSpectre 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 # General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with µSpectre; see the file COPYING. If not, write to the # Free Software Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. # # Additional permission under GNU GPL version 3 section 7 # # If you modify this Program, or any covered work, by linking or combining it # with proprietary FFT implementations or numerical libraries, containing parts # covered by the terms of those libraries' licenses, the licensors of this # Program grant you additional permission to convey the resulting work. # ============================================================================= # FIXME! The user should have a choice to configure this path. execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-m" "site" "--user-site" RESULT_VARIABLE _PYTHON_SUCCESS OUTPUT_VARIABLE PYTHON_USER_SITE ERROR_VARIABLE _PYTHON_ERROR_VALUE) if(NOT _PYTHON_SUCCESS MATCHES 0) message(FATAL_ERROR "Python config failure:\n${_PYTHON_ERROR_VALUE}") endif() string(REGEX REPLACE "\n" "" PYTHON_USER_SITE ${PYTHON_USER_SITE}) set (PY_BINDING_SRCS ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_module.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_common.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_cell.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic1.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic2.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic3.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic4.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_hyper_elasto_plastic1.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic_generic.cc + ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_stochastic_plasticity.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_solvers.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_fftengine.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_projections.cc ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_field_collection.cc ) if (${USE_FFTWMPI}) add_definitions(-DWITH_FFTWMPI) endif(${USE_FFTWMPI}) if (${USE_PFFT}) add_definitions(-DWITH_PFFT) endif(${USE_PFFT}) find_package(PythonLibsNew ${MUSPECTRE_PYTHON_MAJOR_VERSION} MODULE REQUIRED) # On OS X, add -Wno-deprecated-declarations IF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") add_compile_options(-Wno-deprecated-declarations) ENDIF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") pybind11_add_module(pyMuSpectreLib ${PY_BINDING_SRCS}) target_link_libraries(pyMuSpectreLib PRIVATE muSpectre) # Want to rename the output, so that the python module is called muSpectre set_target_properties(pyMuSpectreLib PROPERTIES OUTPUT_NAME _muSpectre) target_include_directories(pyMuSpectreLib PUBLIC ${PYTHON_INCLUDE_DIRS}) add_custom_target(pyMuSpectre ALL SOURCES muSpectre/__init__.py muSpectre/fft.py muSpectre/tools.py) add_custom_command(TARGET pyMuSpectre POST_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/language_bindings/python/muSpectre $/muSpectre) install(TARGETS pyMuSpectreLib LIBRARY DESTINATION ${PYTHON_USER_SITE}) install(FILES muSpectre/__init__.py muSpectre/fft.py DESTINATION ${PYTHON_USER_SITE}/muSpectre) diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc index bfeef87..3f7f5c6 100644 --- a/language_bindings/python/bind_py_material.cc +++ b/language_bindings/python/bind_py_material.cc @@ -1,237 +1,240 @@ /** * @file bind_py_material.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for µSpectre's materials * * Copyright © 2018 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/common.hh" #include "materials/material_base.hh" #include "materials/material_evaluator.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic1_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic2_helper(py::module & mod); /** * python binding for the optionally objective form of Hooke's law */ template void add_material_linear_elastic1_helper(py::module & mod); template void add_material_linear_elastic2_helper(py::module & mod); template void add_material_linear_elastic3_helper(py::module & mod); template void add_material_linear_elastic4_helper(py::module & mod); template void add_material_hyper_elasto_plastic1_helper(py::module & mod); +template +void add_material_stochastic_plasticity_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template class PyMaterialBase : public muSpectre::MaterialBase { public: /* Inherit the constructors */ using Parent = muSpectre::MaterialBase; using Parent::Parent; /* Trampoline (need one for each virtual function) */ void save_history_variables() override { PYBIND11_OVERLOAD_PURE(void, // Return type Parent, // Parent class save_history_variables); // Name of function in C++ // (must match Python name) } /* Trampoline (need one for each virtual function) */ void initialise() override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class initialise); // Name of function in C++ (must match Python name) } void compute_stresses(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, muSpectre::Formulation form) override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class compute_stresses, // Name of function in C++ (must match Python name) F, P, form); } void compute_stresses_tangent(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, typename Parent::TangentField_t & K, muSpectre::Formulation form) override { PYBIND11_OVERLOAD_PURE( void, /* Return type */ Parent, /* Parent class */ compute_stresses, /* Name of function in C++ (must match Python name) */ F, P, K, form); } }; template void add_material_evaluator(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialEvaluator_" << Dim << "d"; std::string name{name_stream.str()}; using MatEval_t = muSpectre::MaterialEvaluator; py::class_(mod, name.c_str()) .def(py::init>>()) .def("save_history_variables", &MatEval_t::save_history_variables, "for materials with state variables") .def("evaluate_stress", [](MatEval_t & mateval, py::EigenDRef & grad, muSpectre::Formulation form) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return mateval.evaluate_stress(grad, form); }, "strain"_a, "formulation"_a, "Evaluates stress for a given strain and formulation " "(Piola-Kirchhoff 1 stress as a function of the placement gradient " "P = P(F) for formulation=Formulation.finite_strain and Cauchy " "stress as a function of the infinitesimal strain tensor σ = σ(ε) " "for formulation=Formulation.small_strain)") .def("evaluate_stress_tangent", [](MatEval_t & mateval, py::EigenDRef & grad, muSpectre::Formulation form) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return mateval.evaluate_stress_tangent(grad, form); }, "strain"_a, "formulation"_a, "Evaluates stress and tangent moduli for a given strain and " "formulation (Piola-Kirchhoff 1 stress as a function of the " "placement gradient P = P(F) for " "formulation=Formulation.finite_strain and Cauchy stress as a " "function of the infinitesimal strain tensor σ = σ(ε) for " "formulation=Formulation.small_strain). The tangent moduli are K = " "∂P/∂F for formulation=Formulation.finite_strain and C = ∂σ/∂ε for " "formulation=Formulation.small_strain.") .def("estimate_tangent", [](MatEval_t & evaluator, py::EigenDRef & grad, muSpectre::Formulation form, const Real step, const muSpectre::FiniteDiff diff_type) { if ((grad.cols() != Dim) or (grad.rows() != Dim)) { std::stringstream err{}; err << "need matrix of shape (" << Dim << "×" << Dim << ") but got (" << grad.rows() << "×" << grad.cols() << ")."; throw std::runtime_error(err.str()); } return evaluator.estimate_tangent(grad, form, step, diff_type); }, "strain"_a, "formulation"_a, "Delta_x"_a, "difference_type"_a = muSpectre::FiniteDiff::centred, "Numerical estimate of the tangent modulus using finite " "differences. The finite difference scheme as well as the finite " "step size can be chosen. If there are no special circumstances, " "the default scheme of centred finite differences yields the most " "accurate results at an increased computational cost."); } template void add_material_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialBase_" << Dim << "d"; std::string name{name_stream.str()}; using Material = muSpectre::MaterialBase; using MaterialTrampoline = PyMaterialBase; using FC_t = muSpectre::LocalFieldCollection; using FCBase_t = muSpectre::FieldCollectionBase; py::class_>(mod, name.c_str()) .def(py::init()) .def("save_history_variables", &Material::save_history_variables) .def("list_fields", &Material::list_fields) .def("get_real_field", &Material::get_real_field, "field_name"_a, py::return_value_policy::reference_internal) .def("size", &Material::size) .def("add_pixel", [](Material & mat, muSpectre::Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a) .def_property_readonly("collection", [](Material & material) -> FCBase_t & { return material.get_collection(); }, "returns the field collection containing internal " "fields of this material"); add_material_linear_elastic1_helper(mod); add_material_linear_elastic2_helper(mod); add_material_linear_elastic3_helper(mod); add_material_linear_elastic4_helper(mod); add_material_hyper_elasto_plastic1_helper(mod); add_material_linear_elastic_generic1_helper(mod); add_material_linear_elastic_generic2_helper(mod); + add_material_stochastic_plasticity_helper(mod); add_material_evaluator(mod); } void add_material(py::module & mod) { auto material{mod.def_submodule("material")}; material.doc() = "bindings for constitutive laws"; add_material_helper(material); add_material_helper(material); } diff --git a/language_bindings/python/bind_py_material_stochastic_plasticity.cc b/language_bindings/python/bind_py_material_stochastic_plasticity.cc new file mode 100644 index 0000000..d262d74 --- /dev/null +++ b/language_bindings/python/bind_py_material_stochastic_plasticity.cc @@ -0,0 +1,104 @@ +/** + * @file bind_py_material_stochastic_plasticity.cc + * + * @author Richard Leute + * + * @date 25 Jan 2019 + * + * @brief python binding for MaterialStochasticPlasticity + * + * Copyright © 2019 Till Junge, Richard Leute + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre 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 + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with µSpectre; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + * + * Additional permission under GNU GPL version 3 section 7 + * + * If you modify this Program, or any covered work, by linking or combining it + * with proprietary FFT implementations or numerical libraries, containing parts + * covered by the terms of those libraries' licenses, the licensors of this + * Program grant you additional permission to convey the resulting work. + */ + + +#include "common/common.hh" +#include "materials/material_stochastic_plasticity.hh" +#include "cell/cell_base.hh" +#include "common/field.hh" + +#include +#include +#include + +#include +#include + +using muSpectre::Dim_t; +using muSpectre::Real; +using pybind11::literals::operator""_a; +namespace py = pybind11; + +/** + * python binding for the optionally objective form of Hooke's law + * with per-pixel elastic properties + */ +template +void add_material_stochastic_plasticity_helper(py::module & mod) { + std::stringstream name_stream{}; + name_stream << "MaterialStochasticPlasticity_" << dim << 'd'; + const auto name{name_stream.str()}; + + using Mat_t = muSpectre::MaterialStochasticPlasticity; + using Sys_t = muSpectre::CellBase; + using StressField_t = + muSpectre::TensorField, Real, + muSpectre::secondOrder, dim>; + + py::class_, std::shared_ptr>( + mod, name.c_str()) + .def(py::init(), "name"_a) + .def_static("make", + [](Sys_t & sys, std::string n) -> Mat_t & { + return Mat_t::make(sys, n); + }, + "cell"_a, "name"_a, py::return_value_policy::reference, + py::keep_alive<1, 0>()) + .def("add_pixel", + [](Mat_t & mat, muSpectre::Ccoord_t pix, Real Young, + Real Poisson, Real plastic_increment, Real stress_threshold, + Eigen::Ref< + const Eigen::Matrix> + eigen_strain) { + mat.add_pixel(pix, Young, Poisson, plastic_increment, + stress_threshold, eigen_strain); + }, + "pixel"_a, "Young"_a, "Poisson"_a, "increment"_a, "threshold"_a, + "strain"_a) + .def_static("make_evaluator", []() { return Mat_t::make_evaluator(); }) + .def("identify_overloaded_pixels", + [](Mat_t & mat, StressField_t & stress) { + return mat.identify_overloaded_pixels( + stress); // TODO(RLeute): Here I + // need the right + // type conversion + // of stress + }, + "stress"_a); +} + +template void +add_material_stochastic_plasticity_helper(py::module &); +template void +add_material_stochastic_plasticity_helper(py::module &); diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 0cc3384..1470e00 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,721 +1,721 @@ /** * @file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * Copyright © 2017 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_ #define SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_ #include "common/common.hh" #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "materials/material_evaluator.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common//utilities.hh" #include #include #include #include namespace muSpectre { // Forward declaration for factory function template class CellBase; /** * material traits are used by `muSpectre::MaterialMuSpectre` to * break the circular dependence created by the curiously recurring * template parameter. These traits must define * - these `muSpectre::FieldMap`s: * - `StrainMap_t`: typically a `muSpectre::MatrixFieldMap` for a * constant second-order `muSpectre::TensorField` * - `StressMap_t`: typically a `muSpectre::MatrixFieldMap` for a * writable secord-order `muSpectre::TensorField` * - `TangentMap_t`: typically a `muSpectre::T4MatrixFieldMap` for a * writable fourth-order `muSpectre::TensorField` * - `strain_measure`: the expected strain type (will be replaced by the * small-strain tensor ε * `muspectre::StrainMeasure::Infinitesimal` in small * strain computations) * - `stress_measure`: the measure of the returned stress. Is used by * `muspectre::MaterialMuSpectre` to transform it into * Cauchy stress (`muspectre::StressMeasure::Cauchy`) in * small-strain computations and into first * Piola-Kirchhoff stress `muspectre::StressMeasure::PK1` * in finite-strain computations * - `InternalVariables`: a tuple of `muSpectre::FieldMap`s containing * internal variables */ template struct MaterialMuSpectre_traits {}; template class MaterialMuSpectre; /** * Base class for most convenient implementation of materials */ template class MaterialMuSpectre : public MaterialBase { public: /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = MatTB::NeedTangent; using Parent = MaterialBase; //!< base class //! global field collection using GFieldCollection_t = typename Parent::GFieldCollection_t; //! expected type for stress fields using StressField_t = typename Parent::StressField_t; //! expected type for strain fields using StrainField_t = typename Parent::StrainField_t; //! expected type for tangent stiffness fields using TangentField_t = typename Parent::TangentField_t; //! traits for the CRTP subclass using traits = MaterialMuSpectre_traits; //! Default constructor MaterialMuSpectre() = delete; //! Construct by name explicit MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre & other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre && other) = delete; //! Destructor virtual ~MaterialMuSpectre() = default; //! Factory template static Material & make(CellBase & cell, ConstructorArgs &&... args); /** Factory * takes all arguments after the name of the underlying * Material's constructor. E.g., if the underlying material is a * `muSpectre::MaterialLinearElastic1`, these would be Young's * modulus and Poisson's ratio. */ template static std::tuple, MaterialEvaluator> make_evaluator(ConstructorArgs &&... args); //! Copy assignment operator MaterialMuSpectre & operator=(const MaterialMuSpectre & other) = delete; //! Move assignment operator MaterialMuSpectre & operator=(MaterialMuSpectre && other) = delete; //! allocate memory, etc void initialise() override; using Parent::compute_stresses; using Parent::compute_stresses_tangent; //! computes stress void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form) final; //! computes stress and tangent modulus void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) final; protected: //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P) __attribute__((visibility("default"))); //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K) __attribute__((visibility("default"))); //! this iterable class is a default for simple laws that just take a strain //! the iterable is just a templated wrapper to provide a range to iterate //! over that does or does not include tangent moduli template class iterable_proxy; /** * inheriting classes with internal variables need to overload this function */ typename traits::InternalVariables & get_internals() { return static_cast(*this).get_internals(); } bool is_initialised{false}; //!< to handle double initialisation right private: }; /* ---------------------------------------------------------------------- */ template MaterialMuSpectre::MaterialMuSpectre(std::string name) : Parent(name) { using stress_compatible = typename traits::StressMap_t::template is_compatible; using strain_compatible = typename traits::StrainMap_t::template is_compatible; using tangent_compatible = typename traits::TangentMap_t::template is_compatible; static_assert((stress_compatible::value && stress_compatible::explain()), "The material's declared stress map is not compatible " "with the stress field. More info in previously shown " "assert."); static_assert((strain_compatible::value && strain_compatible::explain()), "The material's declared strain map is not compatible " "with the strain field. More info in previously shown " "assert."); static_assert((tangent_compatible::value && tangent_compatible::explain()), "The material's declared tangent map is not compatible " "with the tangent field. More info in previously shown " "assert."); } /* ---------------------------------------------------------------------- */ template template Material & MaterialMuSpectre::make(CellBase & cell, ConstructorArgs &&... args) { auto mat = std::make_unique(args...); auto & mat_ref = *mat; cell.add_material(std::move(mat)); return mat_ref; } /* ---------------------------------------------------------------------- */ template template std::tuple, MaterialEvaluator> MaterialMuSpectre::make_evaluator( ConstructorArgs &&... args) { auto mat = std::make_shared("name", args...); using Ret_t = std::tuple, MaterialEvaluator>; return Ret_t(mat, MaterialEvaluator{mat}); } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre::initialise() { if (!this->is_initialised) { this->internal_fields.initialise(); this->is_initialised = true; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre::compute_stresses( const StrainField_t & F, StressField_t & P, Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre::compute_stresses_tangent( const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P, K); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P, K); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre::compute_stresses_worker( const StrainField_t & F, StressField_t & P, TangentField_t & K) { /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this](Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent // moduli Stresses = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...); }, internal_variables); }; auto constitutive_law_finite_strain = [this](Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // TODO(junge): Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. // return value contains a tuple of rvalue_refs to both stress // and tangent moduli auto stress_tgt = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...); }, internal_variables); auto & stress = std::get<0>(stress_tgt); auto & tangent = std::get<1>(stress_tgt); Stresses = MatTB::PK1_stress( std::move(grad), std::move(stress), std::move(tangent)); }; iterable_proxy fields{*this, F, P, K}; for (auto && arglist : fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ // auto && stress_tgt = std::get<0>(tuples); // auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same( std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same( std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); static_assert(std::is_same( std::get<1>(arglist)))>>::value, "Type mismatch for tangent reference, check iterator" "value_type"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre::compute_stresses_worker( const StrainField_t & F, StressField_t & P) { /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this](Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent // moduli auto & sigma = std::get<0>(Stresses); sigma = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress(std::move(strain), internals...); }, internal_variables); }; auto constitutive_law_finite_strain = [this](Strains_t Strains, Stresses_t && Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // TODO(junge): Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. // return value contains a tuple of rvalue_refs to both stress // and tangent moduli auto && stress = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress(std::move(strain), internals...); }, internal_variables); auto & P = get<0>(Stresses); P = MatTB::PK1_stress( F, stress); }; iterable_proxy fields{*this, F, P}; for (auto && arglist : fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ // auto && stress_tgt = std::get<0>(tuples); // auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same( std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same( std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template template class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename MaterialMuSpectre::NeedTangent; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, StressField_t & P, std::enable_if_t & K) : material{mat}, strain_field{F}, stress_tup{P, K}, internals{material.get_internals()} {}; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P) : material{mat}, strain_field{F}, stress_tup{P}, internals{material.get_internals()} {}; //! Expected type for strain fields using StrainMap_t = typename traits::StrainMap_t; //! Expected type for stress fields using StressMap_t = typename traits::StressMap_t; //! Expected type for tangent stiffness fields using TangentMap_t = typename traits::TangentMap_t; //! expected type for strain values using Strain_t = typename traits::StrainMap_t::reference; //! expected type for stress values using Stress_t = typename traits::StressMap_t::reference; //! expected type for tangent stiffness values using Tangent_t = typename traits::TangentMap_t::reference; //! tuple of intervnal variables, depends on the material using InternalVariables = typename traits::InternalVariables; //! tuple containing a stress and possibly a tangent stiffness field using StressFieldTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness field map using StressMapTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness value ref using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! Copy constructor iterable_proxy(const iterable_proxy & other) = default; //! Move constructor iterable_proxy(iterable_proxy && other) = default; //! Destructor virtual ~iterable_proxy() = default; //! Copy assignment operator iterable_proxy & operator=(const iterable_proxy & other) = default; //! Move assignment operator iterable_proxy & operator=(iterable_proxy && other) = default; /** * dereferences into a tuple containing strains, and internal * variables, as well as maps to the stress and potentially * stiffness maps where to write the response of a pixel */ class iterator { public: //! type to refer to internal variables owned by a CRTP material using InternalReferences = MatTB::ReferenceTuple_t; //! return type to be unpacked per pixel my the constitutive law using value_type = std::tuple, Stress_tTup, InternalReferences>; using iterator_category = std::forward_iterator_tag; //!< stl conformance //! Default constructor iterator() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ explicit iterator(const iterable_proxy & it, bool begin = true) : it{it}, strain_map{it.strain_field}, stress_map(it.stress_tup), index{begin ? 0 : it.material.internal_fields.size()} {} //! Copy constructor iterator(const iterator & other) = default; //! Move constructor iterator(iterator && other) = default; //! Destructor virtual ~iterator() = default; //! Copy assignment operator iterator & operator=(const iterator & other) = default; //! Move assignment operator iterator & operator=(iterator && other) = default; //! pre-increment inline iterator & operator++(); //! dereference inline value_type operator*(); //! inequality inline bool operator!=(const iterator & other) const; protected: const iterable_proxy & it; //!< ref to the proxy StrainMap_t strain_map; //!< map onto the global strain field //! map onto the global stress field and possibly tangent stiffness StressMapTup stress_map; size_t index; //!< index or pixel currently referred to }; - //! returns iterator to first pixel if this material + //! returns iterator to first pixel in this material iterator begin() { return std::move(iterator(*this)); } //! returns iterator past the last pixel in this material iterator end() { return std::move(iterator(*this, false)); } protected: MaterialMuSpectre & material; //!< reference to the proxied material const StrainField_t & strain_field; //!< cell's global strain field //! references to the global stress field and perhaps tangent stiffness StressFieldTup stress_tup; //! references to the internal variables InternalVariables & internals; private: }; /* ---------------------------------------------------------------------- */ template template bool MaterialMuSpectre::iterable_proxy::iterator:: operator!=(const iterator & other) const { return (this->index != other.index); } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre::template iterable_proxy::iterator & MaterialMuSpectre::iterable_proxy::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre::template iterable_proxy< NeedTgT>::iterator::value_type MaterialMuSpectre::iterable_proxy::iterator::operator*() { const Ccoord_t pixel{ this->it.material.internal_fields.get_ccoord(this->index)}; auto && strain = std::make_tuple(this->strain_map[pixel]); auto && stresses = apply( [&pixel](auto &&... stress_tgt) { return std::make_tuple(stress_tgt[pixel]...); }, this->stress_map); auto && internal = this->it.material.get_internals(); const auto id{this->index}; auto && internals = apply( [id](auto &&... internals_) { return InternalReferences{internals_[id]...}; }, internal); return std::make_tuple(std::move(strain), std::move(stresses), std::move(internals)); } } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_ diff --git a/src/materials/material_stochastic_plasticity.cc b/src/materials/material_stochastic_plasticity.cc index 379d431..77c748a 100644 --- a/src/materials/material_stochastic_plasticity.cc +++ b/src/materials/material_stochastic_plasticity.cc @@ -1,107 +1,107 @@ /** * @file material_stochastic_plasticity.cc * * @author Richard Leute * * @date 24 Jan 2019 * * @brief material for stochastic plasticity as described in Z. Budrikis et al. * Nature Comm. 8:15928, 2017. It only works together with "python * -script", which performes the avalanche loop. This makes the material * slower but more easy to modify and test. * (copied from material_linear_elastic4.cc) * * Copyright © 2019 Till Junge, Richard Leute * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "materials/material_stochastic_plasticity.hh" #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template MaterialStochasticPlasticity::MaterialStochasticPlasticity( std::string name) : Parent{name}, lambda_field{make_field( "local first Lame constant", this->internal_fields)}, mu_field{ make_field("local second Lame constant(shear modulus)", this->internal_fields)}, plastic_increment_field{ make_field("plastic increment", this->internal_fields)}, stress_threshold_field{ make_field("threshold", this->internal_fields)}, eigen_strain_field{ make_field("eigen strain", this->internal_fields)}, - internal_variables{lambda_field.get_const_map(), - mu_field.get_const_map(), eigen_strain_field.get_map()} - { - } + overloaded_pixels{std::vector>()}, + internal_variables{ + lambda_field.get_const_map(), mu_field.get_const_map(), + eigen_strain_field.get_map()} {} /* ---------------------------------------------------------------------- */ template void MaterialStochasticPlasticity:: add_pixel(const Ccoord_t & /*pixel*/) { throw std::runtime_error ("This material needs pixels with Youngs modulus and Poisson ratio."); } /* ---------------------------------------------------------------------- */ template void MaterialStochasticPlasticity::add_pixel( const Ccoord_t & pixel, const Real & Young_modulus, const Real & Poisson_ratio, const Real & plastic_increment, const Real & stress_threshold, const Eigen::Ref> & eigen_strain) { - //check if the users input eigen strain has the right dimension + // check if the users input eigen strain has the right dimension if (eigen_strain.cols() != DimM || eigen_strain.rows() != DimM) { std::stringstream error{}; error << "Got a wrong shape " << std::to_string(eigen_strain.rows()) << "×" << std::to_string(eigen_strain.cols()) << " for the eigen strain matrix.\nI expected the shape: " << std::to_string(DimM) << "×" << std::to_string(DimM); throw std::runtime_error(error.str()); } this->internal_fields.add_pixel(pixel); // store the first(lambda) and second(mu) Lame constant in the field Real lambda = Hooke::compute_lambda(Young_modulus, Poisson_ratio); Real mu = Hooke::compute_mu(Young_modulus, Poisson_ratio); this->lambda_field.push_back(lambda); this->mu_field.push_back(mu); this->plastic_increment_field.push_back(plastic_increment); this->stress_threshold_field.push_back(stress_threshold); const Eigen::Map> strain_map( eigen_strain.data()); this->eigen_strain_field.push_back(strain_map); } template class MaterialStochasticPlasticity; template class MaterialStochasticPlasticity; template class MaterialStochasticPlasticity; } // namespace muSpectre diff --git a/src/materials/material_stochastic_plasticity.hh b/src/materials/material_stochastic_plasticity.hh index 25061ce..b353cf4 100644 --- a/src/materials/material_stochastic_plasticity.hh +++ b/src/materials/material_stochastic_plasticity.hh @@ -1,235 +1,272 @@ /** * @file material_stochastic_plasticity.hh * * @author Richard Leute * * @date 24 Jan 2019 * * @brief material for stochastic plasticity as described in Z. Budrikis et al. * Nature Comm. 8:15928, 2017. It only works together with "python * -script", which performes the avalanche loop. This makes the material * slower but more easy to modify and test. * (copied from material_linear_elastic4.hh) * * Copyright © 2019 Till Junge, Richard Leute * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_STOCHASTIC_PLASTICITY_HH_ #define SRC_MATERIALS_MATERIAL_STOCHASTIC_PLASTICITY_HH_ #include "materials/material_linear_elastic1.hh" #include "common/field.hh" #include "common/tensor_algebra.hh" #include namespace muSpectre { template class MaterialStochasticPlasticity; /** * traits for stochastic plasticity with eigenstrain */ template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! expected map type for strain fields using StrainMap_t = MatrixFieldMap; //! expected map type for stress fields using StressMap_t = MatrixFieldMap; //! expected map type for tangent stiffness fields using TangentMap_t = T4MatrixFieldMap; //! declare what type of strain measure your law takes as input constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; //! declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::PK2}; //! local field_collections used for internals using LFieldColl_t = LocalFieldCollection; //! local Lame constant, plastic increment, stress threshold type using ScalarMap_t = ScalarFieldMap; //! storage type for eigen strain (is updated from outside) using EigenStrainMap_t = MatrixFieldMap; - //! stochastic plasticity internal variables + //! storage type for an overloaded pixel vector + using Pixel_Vector_t = std::vector>; + + //! stochastic plasticity internal variables (Lame 1, Lame 2, eigen strain, + //! overloaded_pixels) using InternalVariables = std::tuple; }; /** * implements stochastic plasticity with an eigenstrain, Lameconstants and * plastic flow per pixel. */ template class MaterialStochasticPlasticity : public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! global field collection + //! Full type for stress fields + using StressField_t = + TensorField, Real, secondOrder, DimM>; + using Stiffness_t = Eigen::TensorFixedSize>; using EigenStrainArg_t = Eigen::Map>; //! traits of this material using traits = MaterialMuSpectre_traits; //! Type of container used for storing eigenstrain using InternalVariables = typename traits::InternalVariables; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! Default constructor MaterialStochasticPlasticity() = delete; //! Construct by name explicit MaterialStochasticPlasticity(std::string name); //! Copy constructor MaterialStochasticPlasticity(const MaterialStochasticPlasticity & other) = delete; //! Move constructor MaterialStochasticPlasticity(MaterialStochasticPlasticity && other) = delete; //! Destructor virtual ~MaterialStochasticPlasticity() = default; //! Copy assignment operator MaterialStochasticPlasticity & operator=(const MaterialStochasticPlasticity & other) = delete; //! Move assignment operator MaterialStochasticPlasticity & operator=(MaterialStochasticPlasticity && other) = delete; /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor), the first * Lame constant (lambda) and the second Lame constant (shear modulus/mu). */ template inline decltype(auto) evaluate_stress(s_t && E, const Real & lambda, const Real & mu, const EigenStrainArg_t & eigen_strain); /** * evaluates both second Piola-Kirchhoff stress and stiffness given * the Green-Lagrange strain (or Cauchy stress and stiffness if * called with a small strain tensor), the first Lame constant (lambda) and * the second Lame constant (shear modulus/mu). */ template inline decltype(auto) evaluate_stress_tangent(s_t && E, const Real & lambda, const Real & mu, const EigenStrainArg_t & eigen_strain); /** * return the empty internals tuple */ InternalVariables & get_internals() { return this->internal_variables; } /** * overload add_pixel to write into loacal stiffness tensor */ void add_pixel(const Ccoord_t & pixel) final; /** * overload add_pixel to write into local stiffness tensor */ void add_pixel(const Ccoord_t & pixel, const Real & Youngs_modulus, const Real & Poisson_ratio, const Real & plastic_increment, const Real & stress_threshold, const Eigen::Ref> & eigen_strain); + /** + * evaluate how many pixels have a higher stress than their stress threshold + */ + inline std::vector> & + identify_overloaded_pixels(StressField_t & stress_field); + protected: - //! storage for first Lame constant 'lambda' - //! second Lame constant(shear modulus) 'mu' - //! and the plastic flow epsilon_p + //! storage for first Lame constant 'lambda', + //! second Lame constant(shear modulus) 'mu', + //! plastic strain epsilon_p, + //! and a vector of overloaded (stress>stress_threshold) pixel coordinates using Field_t = ScalarField, Real>; using Tensor_Field_t = TensorField, Real, secondOrder, DimM>; Field_t & lambda_field; Field_t & mu_field; Field_t & plastic_increment_field; Field_t & stress_threshold_field; Tensor_Field_t & eigen_strain_field; + std::vector> overloaded_pixels; //! tuple for iterable eigen_field InternalVariables internal_variables; private: }; /* ---------------------------------------------------------------------- */ template template auto MaterialStochasticPlasticity::evaluate_stress( s_t && E, const Real & lambda, const Real & mu, const EigenStrainArg_t & eigen_strain) -> decltype(auto) { - - return Hooke::evaluate_stress(lambda, mu, E-eigen_strain); - + return Hooke::evaluate_stress(lambda, mu, E - eigen_strain); } /* ---------------------------------------------------------------------- */ template template auto MaterialStochasticPlasticity::evaluate_stress_tangent( s_t && E, const Real & lambda, const Real & mu, const EigenStrainArg_t & eigen_strain) -> decltype(auto) { T4Mat C = Hooke::compute_C_T4(lambda, mu); return std::make_tuple( this->evaluate_stress(std::forward(E), lambda, mu, eigen_strain), C); } + /* ---------------------------------------------------------------------- */ + template + std::vector> & + MaterialStochasticPlasticity::identify_overloaded_pixels( + StressField_t & stress_field) { + auto threshold_map{this->stress_threshold_field.get_map()}; + auto stress_map{stress_field.get_map()}; + std::vector> & overloaded_pixels_ref{ + this->overloaded_pixels}; + + // loop over all pixels and check if stress overcomes the threshold or not + for (const auto && pixel_threshold : threshold_map.enumerate()) { + const auto & pixel{std::get<0>(pixel_threshold)}; + const Real & threshold{std::get<1>(pixel_threshold)}; + const auto & stress{stress_map[pixel]}; + // check if stress is larger than threshold + std::cout << stress << threshold << std::endl; + // return & to Ccoord vector with pixels which overcome the threshold, + } + return overloaded_pixels_ref; + } + // relax_overloaded_pixels() //vector leeren + } // namespace muSpectre -#endif // SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC4_HH_ +#endif // SRC_MATERIALS_MATERIAL_STOCHASTIC_PLASTICITY_HH_ diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh index 383dcc6..01a9333 100644 --- a/src/materials/materials_toolbox.hh +++ b/src/materials/materials_toolbox.hh @@ -1,591 +1,603 @@ /** * @file materials_toolbox.hh * * @author Till Junge * * @date 02 Nov 2017 * * @brief collection of common continuum mechanics tools * * Copyright © 2017 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIALS_TOOLBOX_HH_ #define SRC_MATERIALS_MATERIALS_TOOLBOX_HH_ #include "common/common.hh" #include "common/tensor_algebra.hh" #include "common/eigen_tools.hh" #include "common/T4_map_proxy.hh" #include #include #include #include #include #include #include namespace muSpectre { namespace MatTB { /** * thrown when generic materials-related runtime errors occur * (mostly continuum mechanics problems) */ class MaterialsToolboxError : public std::runtime_error { public: //! constructor explicit MaterialsToolboxError(const std::string & what) : std::runtime_error(what) {} //! constructor explicit MaterialsToolboxError(const char * what) : std::runtime_error(what) {} }; /* ---------------------------------------------------------------------- */ /** * Flag used to designate whether the material should compute both stress * and tangent moduli or only stress */ enum class NeedTangent { yes, //!< compute both stress and tangent moduli no //!< compute only stress }; /** * struct used to determine the exact type of a tuple of references obtained * when a bunch of iterators over fiel_maps are dereferenced and their * results are concatenated into a tuple */ template struct ReferenceTuple { //! use this type using type = std::tuple; }; /** * specialisation for tuples */ // template <> template struct ReferenceTuple> { //! use this type using type = typename ReferenceTuple::type; }; /** * helper type for ReferenceTuple */ template using ReferenceTuple_t = typename ReferenceTuple::type; /* ---------------------------------------------------------------------- */ namespace internal { /** Structure for functions returning one strain measure as a * function of another **/ template struct ConvertStrain { static_assert((In == StrainMeasure::Gradient) || (In == StrainMeasure::Infinitesimal), "This situation makes me suspect that you are not using " "MatTb as intended. Disable this assert only if you are " "sure about what you are doing."); //! returns the converted strain template inline static decltype(auto) compute(Strain_t && input) { // transparent case, in which no conversion is required: // just a perfect forwarding static_assert((In == Out), "This particular strain conversion is not implemented"); return std::forward(input); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Green-Lagrange strain from the transformation gradient E = ¹/₂ (C - I) = ¹/₂ (Fᵀ·F - I) **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t && F) { return .5 * (F.transpose() * F - Strain_t::PlainObject::Identity()); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Left Cauchy-Green strain from the transformation gradient B = F·Fᵀ = V² **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t && F) { return F * F.transpose(); } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting Right Cauchy-Green strain from the transformation gradient C = Fᵀ·F = U² **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t && F) { return F.transpose() * F; } }; /* ---------------------------------------------------------------------- */ /** Specialisation for getting logarithmic (Hencky) strain from the transformation gradient E₀ = ¹/₂ ln C = ¹/₂ ln (Fᵀ·F) **/ template <> struct ConvertStrain { //! returns the converted strain template inline static decltype(auto) compute(Strain_t && F) { constexpr Dim_t dim{EigenCheck::tensor_dim::value}; return (.5 * logm(Eigen::Matrix{F.transpose() * F})) .eval(); } }; } // namespace internal /* ---------------------------------------------------------------------- */ //! set of functions returning one strain measure as a function of //! another template decltype(auto) convert_strain(Strain_t && strain) { return internal::ConvertStrain::compute(std::move(strain)); } namespace internal { //! Base template for elastic modulus conversion template struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & /*in1*/, const Real & /*in2*/) { // if no return has happened until now, the conversion is not // implemented yet static_assert( (In1 == In2), "This conversion has not been implemented yet, please add " "it here below as a specialisation of this function " "template. Check " "https://en.wikipedia.org/wiki/Lam%C3%A9_parameters for " "the formula."); return 0; } }; /** * Spectialisation for when the output is the first input */ template struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & A, const Real & /*B*/) { return A; } }; /** * Spectialisation for when the output is the second input */ template struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & /*A*/, const Real & B) { return B; } }; /** * Specialisation μ(E, ν) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & E, const Real & nu) { return E / (2 * (1 + nu)); } }; /** * Specialisation λ(E, ν) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & E, const Real & nu) { return E * nu / ((1 + nu) * (1 - 2 * nu)); } }; /** * Specialisation K(E, ν) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & E, const Real & nu) { return E / (3 * (1 - 2 * nu)); } }; /** * Specialisation E(K, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & K, const Real & G) { return 9 * K * G / (3 * K + G); } }; /** * Specialisation ν(K, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & K, const Real & G) { return (3 * K - 2 * G) / (2 * (3 * K + G)); } }; /** * Specialisation E(λ, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & lambda, const Real & G) { return G * (3 * lambda + 2 * G) / (lambda + G); } }; /** * Specialisation λ(K, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & K, const Real & mu) { return K - 2. * mu / 3.; } }; /** * Specialisation K(λ, µ) */ template <> struct Converter { //! wrapped function (raison d'être) inline constexpr static Real compute(const Real & lambda, const Real & G) { return lambda + (2 * G) / 3; } }; } // namespace internal /** * allows the conversion from any two distinct input moduli to a * chosen output modulus */ template inline constexpr Real convert_elastic_modulus(const Real & in1, const Real & in2) { // enforcing sanity static_assert((In1 != In2), "The input modulus types cannot be identical"); // enforcing independence from order in which moduli are supplied constexpr bool inverted{In1 > In2}; using Converter = std::conditional_t, internal::Converter>; if (inverted) { return Converter::compute(std::move(in2), std::move(in1)); } else { return Converter::compute(std::move(in1), std::move(in2)); } } //! static inline implementation of Hooke's law template struct Hooke { /** * compute Lamé's first constant * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_lambda(const Real & young, const Real & poisson) { return convert_elastic_modulus(young, poisson); } /** * compute Lamé's second constant (i.e., shear modulus) * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_mu(const Real & young, const Real & poisson) { return convert_elastic_modulus(young, poisson); } /** * compute the bulk modulus * @param young: Young's modulus * @param poisson: Poisson's ratio */ inline static constexpr Real compute_K(const Real & young, const Real & poisson) { return convert_elastic_modulus(young, poisson); } /** * compute the stiffness tensor * @param lambda: Lamé's first constant * @param mu: Lamé's second constant (i.e., shear modulus) */ inline static Eigen::TensorFixedSize> compute_C(const Real & lambda, const Real & mu) { return lambda * Tensors::outer(Tensors::I2(), Tensors::I2()) + 2 * mu * Tensors::I4S(); } /** * compute the stiffness tensor * @param lambda: Lamé's first constant * @param mu: Lamé's second constant (i.e., shear modulus) */ inline static T4Mat compute_C_T4(const Real & lambda, const Real & mu) { return lambda * Matrices::Itrac() + 2 * mu * Matrices::Isymm(); } /** * return stress * @param lambda: First Lamé's constant * @param mu: Second Lamé's constant (i.e. shear modulus) * @param E: Green-Lagrange or small strain tensor */ template inline static decltype(auto) evaluate_stress(const Real & lambda, const Real & mu, s_t && E) { return E.trace() * lambda * Strain_t::Identity() + 2 * mu * E; } /** * return stress and tangent stiffness * @param lambda: First Lamé's constant * @param mu: Second Lamé's constant (i.e. shear modulus) * @param E: Green-Lagrange or small strain tensor * @param C: stiffness tensor (Piola-Kirchhoff 2 (or σ) w.r.t to `E`) */ template inline static decltype(auto) evaluate_stress(const Real & lambda, const Real & mu, Tangent_t && C, s_t && E) { return std::make_tuple( std::move(evaluate_stress(lambda, mu, std::move(E))), std::move(C)); } }; namespace internal { /* ---------------------------------------------------------------------- */ template struct NumericalTangentHelper { using T4_t = T4Mat; using T2_t = Eigen::Matrix; using T2_vec = Eigen::Map>; template static inline T4_t compute(FunType && fun, const Eigen::MatrixBase & strain, Real delta); }; /* ---------------------------------------------------------------------- */ template template auto NumericalTangentHelper::compute( FunType && fun, const Eigen::MatrixBase & strain, Real delta) -> T4_t { static_assert((FinDif == FiniteDiff::forward) or (FinDif == FiniteDiff::backward), "Not implemented"); T4_t tangent{T4_t::Zero()}; const T2_t fun_val{fun(strain)}; for (Dim_t i{}; i < Dim * Dim; ++i) { T2_t strain2{strain}; T2_vec strain_vec{strain2.data()}; switch (FinDif) { case FiniteDiff::forward: { strain_vec(i) += delta; T2_t del_f_del{(fun(strain2) - fun_val) / delta}; tangent.col(i) = T2_vec(del_f_del.data()); break; } case FiniteDiff::backward: { strain_vec(i) -= delta; T2_t del_f_del{(fun_val - fun(strain2)) / delta}; tangent.col(i) = T2_vec(del_f_del.data()); break; } } static_assert(Int(decltype(tangent.col(i))::SizeAtCompileTime) == Int(T2_t::SizeAtCompileTime), "wrong column size"); } return tangent; } /* ---------------------------------------------------------------------- */ template struct NumericalTangentHelper { using T4_t = T4Mat; using T2_t = Eigen::Matrix; using T2_vec = Eigen::Map>; template static inline T4_t compute(FunType && fun, const Eigen::MatrixBase & strain, Real delta) { T4_t tangent{T4_t::Zero()}; for (Dim_t i{}; i < Dim * Dim; ++i) { T2_t strain1{strain}; T2_t strain2{strain}; T2_vec strain1_vec{strain1.data()}; T2_vec strain2_vec{strain2.data()}; strain1_vec(i) += delta; strain2_vec(i) -= delta; T2_t del_f_del{(fun(strain1).eval() - fun(strain2).eval()) / (2 * delta)}; tangent.col(i) = T2_vec(del_f_del.data()); static_assert(Int(decltype(tangent.col(i))::SizeAtCompileTime) == Int(T2_t::SizeAtCompileTime), "wrong column size"); } return tangent; } }; } // namespace internal /** * Helper function to numerically determine tangent, intended for * testing, rather than as a replacement for analytical tangents */ template inline T4Mat compute_numerical_tangent( FunType && fun, const Eigen::MatrixBase & strain, Real delta) { static_assert(Derived::RowsAtCompileTime == Dim, "can't handle dynamic matrix"); static_assert(Derived::ColsAtCompileTime == Dim, "can't handle dynamic matrix"); using T2_t = Eigen::Matrix; using T2_vec = Eigen::Map>; static_assert( std::is_convertible>::value, "Function argument 'fun' needs to be a function taking " "one second-rank tensor as input and returning a " "second-rank tensor"); static_assert(Dim_t(T2_t::SizeAtCompileTime) == Dim_t(T2_vec::SizeAtCompileTime), "wrong map size"); return internal::NumericalTangentHelper::compute( std::forward(fun), strain, delta); } + /** + * Computes the equivalent von Mises stress σ_{eq} on each pixel from a given + * Stress. + */ + /* inline compute_equivalent_von_Mises_stress(stress) { + //! compute Couchy stress σ + + //! compute deviatoric stress tensor σ^{dev}=σ-\frac{1}{3} tr(σ) I + + //! compute σ_{eq} = \sqrt{\frac{3}{2} σ^{dev} : σ^{dev}} + return equivalent_stress; + }*/ } // namespace MatTB } // namespace muSpectre #endif // SRC_MATERIALS_MATERIALS_TOOLBOX_HH_ diff --git a/tests/python_material_stochastic_plasticity_test.py b/tests/python_material_stochastic_plasticity_test.py index 4c497fe..837ffc9 100644 --- a/tests/python_material_stochastic_plasticity_test.py +++ b/tests/python_material_stochastic_plasticity_test.py @@ -1,84 +1,107 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file python_material_stochastic_plasticity.py @author Richard Leute @date 25 Jan 2019 @brief tests material_stochastic_plasticity for a single pixel. @section LICENSE Copyright © 2019 Till Junge µSpectre is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3, or (at your option) any later version. µSpectre 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 General Public License for more details. You should have received a copy of the GNU General Public License along with µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Additional permission under GNU GPL version 3 section 7 If you modify this Program, or any covered work, by linking or combining it with proprietary FFT implementations or numerical libraries, containing parts covered by the terms of those libraries' licenses, the licensors of this Program grant you additional permission to convey the resulting work. """ import unittest import numpy as np from muSpectre import material from muSpectre import Formulation from python_test_imports import µ class MaterialStochasticPlasticity_Check(unittest.TestCase): def setUp(self): self.young, self.poisson = 210e9, .33 - self.yield_stress = 1.6e9 + self.plastic_increment = 0.1 + self.yield_stress = 1.6e9 + self.eigen_strain = np.array([[.01, .002], [.002, 0.]]) + + #compute first (lam) and second (mu) Lame constants + #self.lam = (self.young * self.poisson) / \ + # ((1+self.poisson)*(1-2*self.poisson)) + #self.mu = self.young / (2*(1+self.poisson)) self.dim = 2 #spatial dimension self.rank_2_unit = np.eye(self.dim) def test_plastic_deformation(self): #Set up a material evaluator (first test for material_linear_elastic2, #later replace by material_stochastic_plasticity) LinMat2 = material.MaterialLinearElastic2_2d - StochPlastMat = material.MaterialStochasticPlasticity2_2d + StochPlastMat = material.MaterialStochasticPlasticity_2d #the factory returns a material and it's corresponding evaluator material2, evaluator2 = LinMat2.make_evaluator(self.young, self.poisson) + materialSP, evaluatorSP = StochPlastMat.make_evaluator() + #the material is empty(i.e., does not have any pixel / voxel), so a pixel #needs to be added.The coordinates are irrelevant, there just needs to #be one pixel. # when adding the pixel, we now need to specify also the per-pixel data: - eigenstrain = np.array([[.01, .002], [.002, 0.]]) - material2.add_pixel([0,0], eigenstrain) + material2.add_pixel([0,0], self.eigen_strain) + materialSP.add_pixel([0,0], self.young, self.poisson, + self.plastic_increment, self.yield_stress, + self.eigen_strain) # the stress and tangent can be evaluated for finite strain F = np.array([[1., .01],[0, 1.0]]) - P, K = evaluator2.evaluate_stress_tangent(F, Formulation.finite_strain) + P_L2, K_L2 = \ + evaluator2.evaluate_stress_tangent(F, Formulation.finite_strain) + print("P_L2:\n", P_L2) + P_SP, K_SP = \ + evaluatorSP.evaluate_stress_tangent(F, Formulation.finite_strain) + print("P_SP:\n", P_SP) #Compute the criteria for a plastic deformation, #compute the Couchy stress (sigma=1/J*P.F^T) ,J=det(F) - sigma = 1/np.linalg.det(F) * np.dot(P, F.T) + sigma = 1/np.linalg.det(F) * np.dot(P_SP, F.T) #deviatoric stress (sigma_dev=sigma-1/D*Tr(sigma)*I) sigma_dev = sigma - 1/self.dim * np.trace(sigma) * self.rank_2_unit #the equivalent stress (sigma_eq=sqrt(3/2*sigma_dev:sigma_dev)) sigma_eq = np.sqrt(3/2 * np.tensordot(sigma_dev, sigma_dev)) print("equivalent stress: ", sigma_eq, "\nyield stress: ", self.yield_stress) + + + # Check member function "identify_overloaded_pixels" + # I first have to handle the conversion of numpy.ndarray -> muSpectre:: + # TensorField, double, 2, 2> + #overloaded_pixels = materialSP.identify_overloaded_pixels(P_SP) + #print("overloaded_pixels: ", overloaded_pixels)