diff --git a/language_bindings/python/bind_py_material_stochastic_plasticity.cc b/language_bindings/python/bind_py_material_stochastic_plasticity.cc index d262d74..32da50e 100644 --- a/language_bindings/python/bind_py_material_stochastic_plasticity.cc +++ b/language_bindings/python/bind_py_material_stochastic_plasticity.cc @@ -1,104 +1,107 @@ /** * @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>; + + //! dynamic vector type for interactions with numpy/scipy/solvers etc. + using Vector_t = Eigen::Matrix; + using StressField_t = Eigen::Ref; + // 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/common/field_typed.hh b/src/common/field_typed.hh index 77f845d..40a8ef6 100644 --- a/src/common/field_typed.hh +++ b/src/common/field_typed.hh @@ -1,531 +1,531 @@ /** * file field_typed.hh * * @author Till Junge * * @date 10 Apr 2018 * * @brief Typed Field for dynamically sized fields and base class for fields * of tensors, matrices, etc * * 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. */ #ifndef SRC_COMMON_FIELD_TYPED_HH_ #define SRC_COMMON_FIELD_TYPED_HH_ #include "common/field_base.hh" #include "common/field_helpers.hh" #include namespace muSpectre { /** * forward-declaration */ template class TypedFieldMap; namespace internal { /* ---------------------------------------------------------------------- */ //! declaraton for friending template class FieldMap; } // namespace internal /** * Dummy intermediate class to provide a run-time polymorphic * typed field. Mainly for binding Python. TypedField specifies methods * that return typed Eigen maps and vectors in addition to pointers to the * raw data. */ template class TypedField : public internal::FieldBase { friend class internal::FieldMap; friend class internal::FieldMap; static constexpr bool Global{FieldCollection::is_global()}; public: using Parent = internal::FieldBase; //!< base class //! for type checks when mapping this field using collection_t = typename Parent::collection_t; //! for filling global fields from local fields and vice-versa using LocalField_t = std::conditional_t< Global, TypedField, TypedField>; //! for filling global fields from local fields and vice-versa using GlobalField_t = std::conditional_t< Global, TypedField, TypedField>; using Scalar = T; //!< for type checks using Base = Parent; //!< for uniformity of interface //! Plain Eigen type to map using EigenRep_t = Eigen::Array; using EigenVecRep_t = Eigen::Matrix; //! map returned when accessing entire field using EigenMap_t = Eigen::Map; //! map returned when accessing entire const field using EigenMapConst_t = Eigen::Map; //! Plain eigen vector to map using EigenVec_t = Eigen::Map; //! vector map returned when accessing entire field using EigenVecConst_t = Eigen::Map; //! associated non-const field map using FieldMap_t = TypedFieldMap; //! associated const field map using ConstFieldMap_t = TypedFieldMap; /** * type stored (unfortunately, we can't statically size the second * dimension due to an Eigen bug, i.e., creating a row vector * reference to a column vector does not raise an error :( */ using Stored_t = Eigen::Array; //! storage container using Storage_t = std::vector>; //! Default constructor TypedField() = delete; //! constructor TypedField(std::string unique_name, FieldCollection & collection, size_t nb_components); /** * constructor for field proxies which piggy-back on existing * memory. These cannot be registered in field collections and * should only be used for transient temporaries */ TypedField(std::string unique_name, FieldCollection & collection, Eigen::Ref> vec, size_t nb_components); //! Copy constructor TypedField(const TypedField & other) = delete; //! Move constructor TypedField(TypedField && other) = default; //! Destructor virtual ~TypedField() = default; //! Copy assignment operator TypedField & operator=(const TypedField & other) = delete; //! Move assignment operator TypedField & operator=(TypedField && other) = delete; //! return type_id of stored type const std::type_info & get_stored_typeid() const final; //! safe reference cast static TypedField & check_ref(Base & other); //! safe reference cast static const TypedField & check_ref(const Base & other); size_t size() const final; //! add a pad region to the end of the field buffer; required for //! using this as e.g. an FFT workspace void set_pad_size(size_t pad_size_) final; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) void set_zero() final; //! add a new value at the end of the field template inline void push_back(const Eigen::DenseBase & value); //! raw pointer to content (e.g., for Eigen maps) virtual T * data() { return this->get_ptr_to_entry(0); } //! raw pointer to content (e.g., for Eigen maps) virtual const T * data() const { return this->get_ptr_to_entry(0); } //! return a map representing the entire field as a single `Eigen::Array` EigenMap_t eigen(); //! return a map representing the entire field as a single `Eigen::Array` EigenMapConst_t eigen() const; //! return a map representing the entire field as a single Eigen vector EigenVec_t eigenvec(); //! return a map representing the entire field as a single Eigen vector EigenVecConst_t eigenvec() const; //! return a map representing the entire field as a single Eigen vector EigenVecConst_t const_eigenvec() const; /** * Convenience function to return a map onto this field. A map * allows iteration over all pixels. The map's iterator returns a * dynamically sized `Eigen::Map` the data associated with a * pixel. */ inline FieldMap_t get_map(); /** * Convenience function to return a map onto this field. A map * allows iteration over all pixels. The map's iterator returns a * dynamically sized `Eigen::Map` the data associated with a * pixel. */ inline ConstFieldMap_t get_map() const; /** * Convenience function to return a map onto this field. A map * allows iteration over all pixels. The map's iterator returns a * dynamically sized `Eigen::Map` the data associated with a * pixel. */ inline ConstFieldMap_t get_const_map() const; /** * creates a `TypedField` same size and type as this, but all * entries are zero. Convenience function */ inline TypedField & get_zeros_like(std::string unique_name) const; /** * Fill the content of the local field into the global field * (obviously only for pixels that actually are present in the * local field) */ template inline std::enable_if_t fill_from_local(const LocalField_t & local); /** * For pixels that are present in the local field, fill them with * the content of the global field at those pixels */ template inline std::enable_if_t fill_from_global(const GlobalField_t & global); protected: //! returns a raw pointer to the entry, for `Eigen::Map` inline T * get_ptr_to_entry(const size_t && index); //! returns a raw pointer to the entry, for `Eigen::Map` inline const T * get_ptr_to_entry(const size_t && index) const; //! set the storage size of this field inline void resize(size_t size) final; //! The actual storage container Storage_t values{}; /** * an unregistered typed field can be mapped onto an array of * existing values */ optional>> alt_values{}; /** * maintains a tally of the current size, as it cannot be reliably * determined from either `values` or `alt_values` alone. */ size_t current_size; /** * in order to accomodate both registered fields (who own and * manage their data) and unregistered temporary field proxies * (piggy-backing on a chunk of existing memory as e.g., a numpy * array) *efficiently*, the `get_ptr_to_entry` methods need to be * branchless. this means that we cannot decide on the fly whether * to return pointers pointing into values or into alt_values, we * need to maintain an (shudder) raw data pointer that is set * either at construction (for unregistered fields) or at any * resize event (which may invalidate existing pointers). For the * coder, this means that they need to be absolutely vigilant that * *any* operation on the values vector that invalidates iterators * needs to be followed by an update of data_ptr, or we will get * super annoying memory bugs. */ T * data_ptr{}; private: }; } // namespace muSpectre #include "common/field_map_dynamic.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ /* Implementations */ /* ---------------------------------------------------------------------- */ template TypedField::TypedField(std::string unique_name, FieldCollection & collection, size_t nb_components) : Parent(unique_name, nb_components, collection), current_size{0} {} /* ---------------------------------------------------------------------- */ template TypedField::TypedField( std::string unique_name, FieldCollection & collection, Eigen::Ref> vec, size_t nb_components) : Parent(unique_name, nb_components, collection), alt_values{vec}, current_size{vec.size() / nb_components}, data_ptr{vec.data()} { if (vec.size() % nb_components) { std::stringstream err{}; err << "The vector you supplied has a size of " << vec.size() << ", which is not a multiple of the number of components (" << nb_components << ")"; throw FieldError(err.str()); } if (current_size != collection.size()) { std::stringstream err{}; err << "The vector you supplied has the size for " << current_size - << " pixels with " << nb_components << "components each, but the " + << " pixels with " << nb_components << " components each, but the " << "field collection has " << collection.size() << " pixels."; throw FieldError(err.str()); } } /* ---------------------------------------------------------------------- */ //! return type_id of stored type template const std::type_info & TypedField::get_stored_typeid() const { return typeid(T); } /* ---------------------------------------------------------------------- */ template auto TypedField::eigen() -> EigenMap_t { return EigenMap_t(this->data(), this->get_nb_components(), this->size()); } /* ---------------------------------------------------------------------- */ template auto TypedField::eigen() const -> EigenMapConst_t { return EigenMapConst_t(this->data(), this->get_nb_components(), this->size()); } /* ---------------------------------------------------------------------- */ template auto TypedField::eigenvec() -> EigenVec_t { return EigenVec_t(this->data(), this->get_nb_components() * this->size(), 1); } /* ---------------------------------------------------------------------- */ template auto TypedField::eigenvec() const -> EigenVecConst_t { return EigenVecConst_t(this->data(), this->get_nb_components() * this->size(), 1); } /* ---------------------------------------------------------------------- */ template auto TypedField::const_eigenvec() const -> EigenVecConst_t { return EigenVecConst_t(this->data(), this->get_nb_components() * this->size()); } /* ---------------------------------------------------------------------- */ template auto TypedField::get_map() -> FieldMap_t { return FieldMap_t(*this); } /* ---------------------------------------------------------------------- */ template auto TypedField::get_map() const -> ConstFieldMap_t { return ConstFieldMap_t(*this); } /* ---------------------------------------------------------------------- */ template auto TypedField::get_const_map() const -> ConstFieldMap_t { return ConstFieldMap_t(*this); } /* ---------------------------------------------------------------------- */ template auto TypedField::get_zeros_like(std::string unique_name) const -> TypedField & { return make_field(unique_name, this->collection, this->nb_components); } /* ---------------------------------------------------------------------- */ template template std::enable_if_t TypedField::fill_from_local(const LocalField_t & local) { static_assert(IsGlobal == Global, "SFINAE parameter, do not touch"); if (not(local.get_nb_components() == this->get_nb_components())) { std::stringstream err_str{}; err_str << "Fields not compatible: You are trying to write a local " << local.get_nb_components() << "-component field into a global " << this->get_nb_components() << "-component field."; throw std::runtime_error(err_str.str()); } auto this_map{this->get_map()}; auto local_map{local.get_map()}; for (const auto && key_val : local_map.enumerate()) { const auto & key{std::get<0>(key_val)}; const auto & value{std::get<1>(key_val)}; this_map[key] = value; } } /* ---------------------------------------------------------------------- */ template template std::enable_if_t TypedField::fill_from_global( const GlobalField_t & global) { static_assert(IsLocal == not Global, "SFINAE parameter, do not touch"); if (not(global.get_nb_components() == this->get_nb_components())) { std::stringstream err_str{}; err_str << "Fields not compatible: You are trying to write a global " << global.get_nb_components() << "-component field into a local " << this->get_nb_components() << "-component field."; throw std::runtime_error(err_str.str()); } auto global_map{global.get_map()}; auto this_map{this->get_map()}; for (auto && key_val : this_map.enumerate()) { const auto & key{std::get<0>(key_val)}; auto & value{std::get<1>(key_val)}; value = global_map[key]; } } /* ---------------------------------------------------------------------- */ template void TypedField::resize(size_t size) { if (this->alt_values) { throw FieldError("Field proxies can't resize."); } this->current_size = size; this->values.resize(size * this->get_nb_components() + this->pad_size); this->data_ptr = &this->values.front(); } /* ---------------------------------------------------------------------- */ template void TypedField::set_zero() { std::fill(this->values.begin(), this->values.end(), T{}); } /* ---------------------------------------------------------------------- */ template auto TypedField::check_ref(Base & other) -> TypedField & { if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) { std::string err = "Cannot create a Reference of requested type " + ("for field '" + other.get_name() + "' of type '" + other.get_stored_typeid().name() + "'"); throw std::runtime_error(err); } return static_cast(other); } /* ---------------------------------------------------------------------- */ template auto TypedField::check_ref(const Base & other) -> const TypedField & { if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) { std::string err = "Cannot create a Reference of requested type " + ("for field '" + other.get_name() + "' of type '" + other.get_stored_typeid().name() + "'"); throw std::runtime_error(err); } return static_cast(other); } /* ---------------------------------------------------------------------- */ template size_t TypedField::size() const { return this->current_size; } /* ---------------------------------------------------------------------- */ template void TypedField::set_pad_size(size_t pad_size) { if (this->alt_values) { throw FieldError("You can't set the pad size of a field proxy."); } this->pad_size = pad_size; this->resize(this->size()); this->data_ptr = &this->values.front(); } /* ---------------------------------------------------------------------- */ template T * TypedField::get_ptr_to_entry(const size_t && index) { return this->data_ptr + this->get_nb_components() * std::move(index); } /* ---------------------------------------------------------------------- */ template const T * TypedField::get_ptr_to_entry( const size_t && index) const { return this->data_ptr + this->get_nb_components() * std::move(index); } /* ---------------------------------------------------------------------- */ template template void TypedField::push_back( const Eigen::DenseBase & value) { static_assert(not FieldCollection::Global, "You can only push_back data into local field collections"); if (value.cols() != 1) { std::stringstream err{}; err << "Expected a column vector, but received and array with " << value.cols() << " colums."; throw FieldError(err.str()); } if (value.rows() != static_cast(this->get_nb_components())) { std::stringstream err{}; err << "Expected a column vector of length " << this->get_nb_components() << ", but received one of length " << value.rows() << "."; throw FieldError(err.str()); } for (size_t i = 0; i < this->get_nb_components(); ++i) { this->values.push_back(value(i)); } ++this->current_size; this->data_ptr = &this->values.front(); } } // namespace muSpectre #endif // SRC_COMMON_FIELD_TYPED_HH_ diff --git a/src/materials/material_stochastic_plasticity.hh b/src/materials/material_stochastic_plasticity.hh index b353cf4..af102d8 100644 --- a/src/materials/material_stochastic_plasticity.hh +++ b/src/materials/material_stochastic_plasticity.hh @@ -1,272 +1,296 @@ /** * @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; //! 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>; + //! Full type for stress TypedField + //using StressTypedField_t = + // TypedField, Real, secondOrder, DimM>; + + //! dynamic vector type for interactions with numpy/scipy/solvers etc. + using Vector_t = Eigen::Matrix; + 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); + identify_overloaded_pixels(Eigen::Ref & stress_field); protected: //! 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); } /* ---------------------------------------------------------------------- */ 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) { + Eigen::Ref & stress_numpy_array) { auto threshold_map{this->stress_threshold_field.get_map()}; + + //std::cout << "traits: " << traits::GFieldCollection_t << std::endl; + + // convert Eigen::Ref to TypedField_t, then you can use get_map() + GlobalFieldCollection internal_stress_field; // here I should give + // the actual field + // collection of the + // material and not some + // empty one.... + //traits::GFieldCollection_t internal_stress_field; + TypedField, Real> stress_field( + "stress_field", internal_stress_field, stress_numpy_array, DimS * DimS); + 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; + std::cout << "\nYes! I'm here!!!" << std::endl; + std::cout << "stress numpy array:\n" << stress_numpy_array << std::endl; + std::cout << "stress converted:\n" << stress << std::endl; + std::cout << "threshold: " << threshold << std::endl; + std::cout << "pixel: " << pixel << std::endl; // return & to Ccoord vector with pixels which overcome the threshold, } return overloaded_pixels_ref; } - // relax_overloaded_pixels() //vector leeren + // relax_overloaded_pixels() //Funktion um vector zu leeren } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_STOCHASTIC_PLASTICITY_HH_ diff --git a/tests/python_material_stochastic_plasticity_test.py b/tests/python_material_stochastic_plasticity_test.py index 837ffc9..6d607cf 100644 --- a/tests/python_material_stochastic_plasticity_test.py +++ b/tests/python_material_stochastic_plasticity_test.py @@ -1,107 +1,109 @@ #!/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.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.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: 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_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_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) + #convert to accepted shape (m,1) + overloaded_pixels = \ + materialSP.identify_overloaded_pixels(P_SP.reshape(-1,1)) + print("overloaded_pixels: ", overloaded_pixels)