diff --git a/bin/profile_material_linear_elastic4.cc b/bin/profile_material_linear_elastic4.cc index 2064ea4..4068f99 100644 --- a/bin/profile_material_linear_elastic4.cc +++ b/bin/profile_material_linear_elastic4.cc @@ -1,101 +1,105 @@ /** * @file profile_material_linear_elastic4.cc * * @author Richard Leute * * @date 19 Apr 2018 * * @brief program to profile material_linear_elastic4 for possible bottlenecks * * Copyright © 2018 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 GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "cell/cell_factory.hh" #include "materials/material_linear_elastic4.hh" #include "solver/solvers.hh" #include "solver/solver_cg.hh" #include #include using namespace muSpectre; int main() { std::cout << "startig program\n"; const Ccoord_t resolution{25, 25, 25}; const Rcoord_t lengths{5.0, 5.0, 5.0}; constexpr Formulation formulation{Formulation::finite_strain}; auto cell{make_cell(resolution, lengths, formulation)}; //initialize my material using Material_t = MaterialLinearElastic4; auto mat{std::make_unique("material")}; //introduce random number generators for Young and Poisson //make random numbers as in std::uniform_real_distribution cppreference.com std::random_device rd; std::mt19937 gen(rd()); //set a seed for the random number generation to make reproducable results gen.seed(15); // for Young [5,10) std::uniform_real_distribution<> dis_Y(5, 10); //for Poisson [0.1, 0.4) std::uniform_real_distribution<> dis_P(0.1, 0.4); //fill the pixels for (const auto && pixel:cell) { mat->add_pixel(pixel, dis_Y(gen), dis_P(gen)); } cell.add_material(std::move(mat)); //add the material to the cell cell.initialise(FFT_PlanFlags::measure); //fft initialization, make faster fft //set the deformation //DelF = np.array([[0, 0.01, 0], // [0, 0 , 0], // [0, 0 , 0]]) - Grad_t DelF{Grad_t::Zero()}; + //Grad_t DelF{Grad_t::Zero()}; + //DelF(0, 1) = 0.01; + //GradIncrements grads{DelF}; + + + Eigen::MatrixXd DelF{Eigen::MatrixXd::Zero(threeD, threeD)}; DelF(0, 1) = 0.01; //set solver constants const Real newton_tol {1e-6}; //tolerance for newton algo const Real cg_tol {1e-6}; //tolerance for cg algo const Real equil_tol {1e-6}; //tolerance for equilibrium / tol for div(P) ? const Uint maxiter {100}; //maximum cg iterations const Uint verbose {0}; //verbosity of solver - GradIncrements grads{DelF}; - SolverCG cg{cell, cg_tol, maxiter, bool(verbose)}; - auto newton_result {newton_cg(cell, grads, cg, - newton_tol, equil_tol, verbose)[0]}; + SolverCG cg{cell, cg_tol, maxiter, bool(verbose)}; + auto newton_result {newton_cg(cell, DelF, cg, + newton_tol, equil_tol, verbose)}; //print some messages of the solver //std::cout << "gradient F:\n" << newton_result.grad << "\n"; //std::cout << "\n\nstress:\n" << newton_result.stress << "\n"; std::cout << "convergence_test: " << newton_result.success << "\n"; std::cout << "status: " << newton_result.status << "\n"; std::cout << "message: " << newton_result.message << "\n"; std::cout << "# iterations: " << newton_result.nb_it << "\n"; std::cout << "# cell evaluations: " << newton_result.nb_fev << "\n"; std::cout << "program finished correct\n"; } diff --git a/src/materials/material_hyper_elasto_plastic2.cc b/src/materials/material_hyper_elasto_plastic2.cc index 6c054d5..7a9483b 100644 --- a/src/materials/material_hyper_elasto_plastic2.cc +++ b/src/materials/material_hyper_elasto_plastic2.cc @@ -1,107 +1,112 @@ /** * @file material_hyper_elasto_plastic2.cc * * @author Richard Leute * * @date 08 May 2018 * * @brief implementation for MaterialHyperElastoPlastic2 * * Copyright © 2018 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 GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material_hyper_elasto_plastic2.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialHyperElastoPlastic2:: MaterialHyperElastoPlastic2(std::string name) : Parent{name}, - plast_flow_field("cumulated plastic flow εₚ", this->internal_fields), - F_prev_field("Previous placement gradient Fᵗ", this->internal_fields), - be_prev_field("Previous left Cauchy-Green deformation bₑᵗ", - this->internal_fields), + plast_flow_field{make_statefield>> + ("cumulated plastic flow εₚ", this->internal_fields)}, + F_prev_field{make_statefield>> + ("Previous placement gradient Fᵗ", this->internal_fields)}, + be_prev_field{make_statefield>> + ("Previous left Cauchy-Green deformation bₑᵗ", + this->internal_fields)}, lambda_field{make_field("local first Lame constant", this->internal_fields )}, mu_field{make_field("local second lame constant", this->internal_fields )}, tau_y0_field{make_field("local yield stress", this->internal_fields )}, H_field{make_field("local hardening modulus", this->internal_fields )}, internal_variables{F_prev_field.get_map(), be_prev_field.get_map(), plast_flow_field.get_map(), lambda_field.get_const_map(), mu_field.get_const_map(), tau_y0_field.get_const_map(), H_field.get_const_map()} {} /* ---------------------------------------------------------------------- */ template void MaterialHyperElastoPlastic2::save_history_variables() { this->plast_flow_field.cycle(); this->F_prev_field.cycle(); this->be_prev_field.cycle(); } /* ---------------------------------------------------------------------- */ template - void MaterialHyperElastoPlastic2::initialise(bool /*stiffness*/) { + void MaterialHyperElastoPlastic2::initialise() { Parent::initialise(); this->F_prev_field.get_map().current() = Eigen::Matrix::Identity(); this->be_prev_field.get_map().current() = Eigen::Matrix::Identity(); this->save_history_variables(); } /* ---------------------------------------------------------------------- */ template void MaterialHyperElastoPlastic2:: add_pixel(const Ccoord_t & /*pixel*/) { throw std::runtime_error ("This material needs pixels with Youngs modulus, Poisson ratio,\n" "Yield stress and Hardening modulus."); } /* ---------------------------------------------------------------------- */ template void MaterialHyperElastoPlastic2:: add_pixel(const Ccoord_t & pixel, const Real & Youngs_modulus, const Real & Poisson_ratio, const Real & tau_y0, const Real & H) { this->internal_fields.add_pixel(pixel); this->tau_y0_field.push_back(tau_y0); this->H_field.push_back(H); // compute lambda and mu auto lambda{Hooke::compute_lambda(Youngs_modulus, Poisson_ratio)}; auto mu{Hooke::compute_mu(Youngs_modulus, Poisson_ratio)}; this->lambda_field.push_back(lambda); this->mu_field.push_back(mu); } template class MaterialHyperElastoPlastic2< twoD, twoD>; template class MaterialHyperElastoPlastic2< twoD, threeD>; template class MaterialHyperElastoPlastic2; } // muSpectre diff --git a/src/materials/material_hyper_elasto_plastic2.hh b/src/materials/material_hyper_elasto_plastic2.hh index fa722ea..9126af1 100644 --- a/src/materials/material_hyper_elasto_plastic2.hh +++ b/src/materials/material_hyper_elasto_plastic2.hh @@ -1,374 +1,374 @@ /** * @file material_hyper_elasto_plastic2.hh * * @author Richard Leute * * @date 08 May 2018 * * @brief Material for logarithmic hyperelasto-plasticity, as defined in de * Geus 2017 (https://doi.org/10.1016/j.cma.2016.12.032) and further * explained in Geers 2003 (https://doi.org/10.1016/j.cma.2003.07.014). * In difference to material_hyper_elasto_plastic1.hh one can choose * arbitrary material constants (Youngs modulus, Poisson ratio, Yield * stress and Hardening modulus) for each pixel. * * Copyright © 2018 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 GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MATERIAL_HYPER_ELASTO_PLASTIC2_H #define MATERIAL_HYPER_ELASTO_PLASTIC2_H #include "materials/material_muSpectre_base.hh" #include "materials/materials_toolbox.hh" #include "common/eigen_tools.hh" #include "common/statefield.hh" #include namespace muSpectre { template class MaterialHyperElastoPlastic2; /** * traits for hyper-elastoplastic material */ 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::Gradient}; //! declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::Kirchhoff}; //! local field collection used for internals using LFieldColl_t = LocalFieldCollection; //! storage type for plastic flow measure (εₚ in the papers) using LScalarMap_t = StateFieldMap>; /** * storage type for for previous gradient Fᵗ and elastic left * Cauchy-Green deformation tensor bₑᵗ */ using LStrainMap_t = StateFieldMap< MatrixFieldMap>; /** * storage type for the Const Real Material Variables: * lambda, mu, tau_y0 and H. */ using LConstRealMaterialVariableMap_t = ScalarFieldMap; /** * format in which to receive internals (previous gradient Fᵗ, * previous elastic lef Cauchy-Green deformation tensor bₑᵗ, and * the plastic flow measure εₚ */ using InternalVariables = std::tuple; }; /** * Material implementation for hyper-elastoplastic constitutive law */ template class MaterialHyperElastoPlastic2: public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre , DimS, DimM>; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! shortcut to traits using traits = MaterialMuSpectre_traits; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! type in which the previous strain state is referenced using StrainStRef_t = typename traits::LStrainMap_t::reference; //! type in which the previous plastic flow is referenced using FlowStRef_t = typename traits::LScalarMap_t::reference; //! Default constructor MaterialHyperElastoPlastic2() = delete; //! Constructor with name and material properties MaterialHyperElastoPlastic2(std::string name); //! Copy constructor MaterialHyperElastoPlastic2(const MaterialHyperElastoPlastic2 &other) = delete; //! Move constructor MaterialHyperElastoPlastic2(MaterialHyperElastoPlastic2 &&other) = delete; //! Destructor virtual ~MaterialHyperElastoPlastic2() = default; //! Copy assignment operator MaterialHyperElastoPlastic2& operator=(const MaterialHyperElastoPlastic2 &other) = delete; //! Move assignment operator MaterialHyperElastoPlastic2& operator=(MaterialHyperElastoPlastic2 &&other) = delete; /** * evaluates Kirchhoff stress given the current placement gradient * Fₜ, the previous Gradient Fₜ₋₁ and the cumulated plastic flow * εₚ * Additionally the constant fields: lambda (Lames first constant), * mu (shear modulus/Lames second constant), tau_y0 (yield stress) * and H (hardening modulus). */ template inline decltype(auto) evaluate_stress(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t plast_flow, const Real lambda, const Real mu, const Real tau_y0, const Real H); /** * evaluates Kirchhoff stress and stiffness given the current placement gradient * Fₜ, the previous Gradient Fₜ₋₁ and the cumulated plastic flow * εₚ * Additionally the constant fields: lambda (Lames first constant), * mu (shear modulus/Lames second constant), K (bulk modulus), yield stress * and the hardening modulus. */ template inline decltype(auto) evaluate_stress_tangent(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t plast_flow, const Real lambda, const Real mu, const Real tau_y0, const Real H); /** * The statefields need to be cycled at the end of each load increment */ virtual void save_history_variables() override; /** * set the previous gradients to identity */ - virtual void initialise(bool stiffness=false) override final; + virtual void initialise() override final; /** * return the internals tuple */ typename traits::InternalVariables & get_internals() { return this->internal_variables;}; /** * overload add_pixel to write into loacal stiffness tensor */ void add_pixel(const Ccoord_t & pixel) override final; /** * overload add_pixel to write into local stiffness tensor */ void add_pixel(const Ccoord_t & pixel, const Real & Young_modulus, const Real & Poisson_ratio, const Real & tau_y0, const Real & H); protected: /** * worker function computing stresses and internal variables */ template inline decltype(auto) stress_n_internals_worker(grad_t && F, StrainStRef_t& F_prev, StrainStRef_t& be_prev, FlowStRef_t& plast_flow, const Real lambda, const Real mu, const Real tau_y0, const Real H); //! Local FieldCollection type for field storage using LColl_t = LocalFieldCollection; - //! storage for cumulated plastic flow εₚ - StateField> plast_flow_field; + //! storage for cumulated plastic flow εₚ //WHY do I need here References??? + StateField> & plast_flow_field; //! storage for previous gradient Fᵗ - StateField> F_prev_field; + StateField> & F_prev_field; //! storage for elastic left Cauchy-Green deformation tensor bₑᵗ - StateField> be_prev_field; + StateField> & be_prev_field; /** storage for the first Lame constant (lambda), second Lame constant (mu), * Yield_stress (tau_y0) and Hardening_modulus (H). */ using Field_t = MatrixField, Real, oneD, oneD>; Field_t & lambda_field; //!< first Lamé constant Field_t & mu_field; //!< second Lamé constant (shear modulus) Field_t & tau_y0_field; //!< initial yield stress Field_t & H_field; //!< hardening modulus //! Field maps and state field maps over internal fields typename traits::InternalVariables internal_variables; private: }; //----------------------------------------------------------------------------// template template decltype(auto) MaterialHyperElastoPlastic2:: stress_n_internals_worker(grad_t && F, StrainStRef_t& F_prev, StrainStRef_t& be_prev, FlowStRef_t& eps_p, const Real lambda, const Real mu, const Real tau_y0, const Real H) { // the notation in this function follows Geers 2003 // (https://doi.org/10.1016/j.cma.2003.07.014). // computation of trial state using Mat_t = Eigen::Matrix; auto && f{F*F_prev.old().inverse()}; Mat_t be_star{f*be_prev.old()*f.transpose()}; Mat_t ln_be_star{logm(std::move(be_star))}; Mat_t tau_star{.5*Hooke::evaluate_stress(lambda, mu, ln_be_star)}; // deviatoric part of Kirchhoff stress Mat_t tau_d_star{tau_star - tau_star.trace()/DimM*tau_star.Identity()}; auto && tau_eq_star{std::sqrt(3*.5*(tau_d_star.array()* tau_d_star.transpose().array()).sum())}; Mat_t N_star{3*.5*tau_d_star/tau_eq_star}; // this is eq (27), and the std::min enforces the Kuhn-Tucker relation (16) //// MIN or MAX ??? Real phi_star{std::max(tau_eq_star - tau_y0 - H * eps_p.old(), 0.)}; // return mapping Real Del_gamma{phi_star/(H + 3 * mu)}; auto && tau{tau_star - 2*Del_gamma*mu*N_star}; /////auto && tau_eq{tau_eq_star - 3*mu*Del_gamma}; // update the previous values to the new ones F_prev.current() = F; ln_be_star -= 2*Del_gamma*N_star; be_prev.current() = expm(std::move(ln_be_star)); eps_p.current() += Del_gamma; // transmit info whether this is a plastic step or not bool is_plastic{phi_star > 0}; // compute tau as for a linear elastic material for test purpose std::cout << "tau_star*:\n" << tau_star << std::endl; std::cout << "tau_eq_star:\n" << tau_eq_star << std::endl; std::cout << "Del_gamma:\n" << Del_gamma << std::endl; std::cout << "N_star:\n" << N_star << std::endl; return std::tuple (std::move(tau), std::move(tau_eq_star), std::move(Del_gamma), std::move(N_star), std::move(is_plastic)); } //----------------------------------------------------------------------------// template template decltype(auto) MaterialHyperElastoPlastic2:: evaluate_stress(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t eps_p, const Real lambda, const Real mu, const Real tau_y0, const Real H) { auto retval(std::move(std::get<0>(this->stress_n_internals_worker (std::forward(F), F_prev, be_prev, eps_p, lambda, mu, tau_y0, H)))); return retval; } //----------------------------------------------------------------------------// template template decltype(auto) MaterialHyperElastoPlastic2:: evaluate_stress_tangent(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t eps_p, const Real lambda, const Real mu, const Real tau_y0, const Real H) { //! after the stress computation, all internals are up to date auto && vals{this->stress_n_internals_worker (std::forward(F), F_prev, be_prev, eps_p, lambda, mu, tau_y0, H)}; auto && tau {std::get<0>(vals)}; auto && tau_eq_star{std::get<1>(vals)}; auto && Del_gamma {std::get<2>(vals)}; auto && N_star {std::get<3>(vals)}; auto && is_plastic {std::get<4>(vals)}; if (is_plastic) { auto && a0 = Del_gamma* mu/tau_eq_star; auto && a1 = mu/(H + 3*mu); // compute bulk modulus K from first(lambda) and second(mu) Lame constants auto K{Hooke::compute_K(lambda, mu)}; //does this work? const Real K_test{lambda + 2*mu/3}; //K computed by hard for test reasons //test of compute K std::cout << K << " = " << K_test << "?" << std::endl; return std::make_tuple(std::move(tau), T4Mat{ ((K/2. - mu/3 + a0*mu)*Matrices::Itrac() + (1 - 3*a0) * mu*Matrices::Isymm() + 2*mu * (a0 - a1)*Matrices::outer(N_star, N_star))}); } else { // compute stiffness tensor C // the factor .5 comes from equation (18) in Geers 2003 // (https://doi.org/10.1016/j.cma.2003.07.014) auto C{0.5*Hooke::compute_C_T4(lambda, mu)}; return std::make_tuple(std::move(tau), T4Mat{C}); std::cout << "elastic" << std::endl; } } } // muSpectre #endif /* MATERIAL_HYPER_ELASTO_PLASTIC2_H */ diff --git a/src/materials/material_linear_elastic4.hh b/src/materials/material_linear_elastic4.hh index 85c2a57..9e7edcd 100644 --- a/src/materials/material_linear_elastic4.hh +++ b/src/materials/material_linear_elastic4.hh @@ -1,206 +1,211 @@ /** * @file material_linear_elastic4.hh * * @author Richard Leute * * @date 15 March 2018 * * @brief linear elastic material with distribution of stiffness properties. * In difference to material_linear_elastic3 two Lame constants are * stored per pixel instead of the whole elastic matrix C. * Uses the MaterialMuSpectre facilities to keep it simple. * * Copyright © 2018 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 GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_2_H #define MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_2_H #include "materials/material_linear_elastic1.hh" #include "common/field.hh" #include "common/tensor_algebra.hh" #include namespace muSpectre { template class MaterialLinearElastic4; /** * traits for objective linear elasticity 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 type using LLameConstantMap_t = ScalarFieldMap; //! elasticity without internal variables using InternalVariables = std::tuple; }; /** * implements objective linear elasticity with an eigenstrain per pixel */ template class MaterialLinearElastic4: 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 using Stiffness_t = Eigen::TensorFixedSize >; //! 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 MaterialLinearElastic4() = delete; //! Construct by name MaterialLinearElastic4(std::string name); //! Copy constructor MaterialLinearElastic4(const MaterialLinearElastic4 &other) = delete; //! Move constructor MaterialLinearElastic4(MaterialLinearElastic4 &&other) = delete; //! Destructor virtual ~MaterialLinearElastic4() = default; //! Copy assignment operator MaterialLinearElastic4& operator=(const MaterialLinearElastic4 &other) = delete; //! Move assignment operator MaterialLinearElastic4& operator=(MaterialLinearElastic4 &&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); /** * 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); /** * 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) override 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); protected: //! storage for first Lame constant 'lambda' //! and second Lame constant(shear modulus) 'mu' using Field_t = MatrixField, Real, oneD, oneD>; Field_t & lambda_field; Field_t & mu_field; //! tuple for iterable eigen_field InternalVariables internal_variables; private: }; /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic4:: - evaluate_stress(s_t && E, const Real & lambda, const Real & mu) { - return Hooke::evaluate_stress(lambda, mu, std::forward(E)); + evaluate_stress(s_t && E, const Real & lambda, const Real & mu) + -> decltype(auto) { + auto C = Hooke::compute_C_T4(lambda, mu); + return Matrices::tensmult(C, E); } + /*{ + return Hooke::evaluate_stress(lambda, mu, std::forward(E)); + }*/ /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic4:: evaluate_stress_tangent(s_t && E, const Real & lambda, const Real & mu) -> decltype(auto) { T4Mat C = Hooke::compute_C_T4(lambda, mu); return std::make_tuple(this->evaluate_stress(std::forward(E), lambda, mu), C); } } // muSpectre #endif /* MATERIAL_LINEAR_ELASTIC_RANDOM_STIFFNESS_2_H */ diff --git a/src/solver/solver_cg.hh b/src/solver/solver_cg.hh index bc1cc18..fa03b33 100644 --- a/src/solver/solver_cg.hh +++ b/src/solver/solver_cg.hh @@ -1,103 +1,103 @@ /** * file solver_cg.hh * * @author Till Junge * * @date 24 Apr 2018 * * @brief class fo a simple implementation of a conjugate gradient solver. * This follows algorithm 5.2 in Nocedal's Numerical Optimization * (p 112) * * Copyright © 2018 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 GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef SOLVER_CG_H #define SOLVER_CG_H #include "solver/solver_base.hh" namespace muSpectre { /** * implements the `muSpectre::SolverBase` interface using a * conjugate gradient solver. This particular class is useful for * trouble shooting, as it can be made very verbose, but for * production runs, it is probably better to use * `muSpectre::SolverCGEigen`. */ class SolverCG: public SolverBase { public: using Parent = SolverBase; //!< standard short-hand for base class //! for storage of fields using Vector_t = Parent::Vector_t; //! Input vector for solvers using Vector_ref = Parent::Vector_ref; //! Input vector for solvers using ConstVector_ref = Parent::ConstVector_ref; //! Output vector for solvers using Vector_map = Parent::Vector_map; //! Default constructor SolverCG() = delete; //! Copy constructor SolverCG(const SolverCG &other) = delete; /** * Constructor takes a Cell, tolerance, max number of iterations * and verbosity flag as input */ SolverCG(Cell & cell, Real tol, Uint maxiter, bool verbose=false); //! Move constructor SolverCG(SolverCG &&other) = default; //! Destructor virtual ~SolverCG() = default; //! Copy assignment operator SolverCG& operator=(const SolverCG &other) = delete; //! Move assignment operator SolverCG& operator=(SolverCG &&other) = default; //! initialisation does not need to do anything in this case void initialise() override final {}; //! returns the solver's name std::string get_name() const override final {return "CG";} //! the actual solver Vector_map solve(const ConstVector_ref rhs) override final; - + protected: Vector_t r_k; //!< residual Vector_t p_k; //!< search direction Vector_t Ap_k; //!< directional stiffness Vector_t x_k; //!< current solution private: }; } // muSpectre #endif /* SOLVER_CG_H */ diff --git a/tests/python_muSpectre_tools_test.py b/tests/python_muSpectre_tools_test.py index cff8712..4c83e1f 100644 --- a/tests/python_muSpectre_tools_test.py +++ b/tests/python_muSpectre_tools_test.py @@ -1,750 +1,750 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file python_muSpectre_tools_test.py @author Richard Leute @date 24 May 2018 @brief test the behaviour of tools.py @section LICENSE Copyright © 2018 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 GNU Emacs; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. """ import unittest import numpy as np import scipy.misc as sm import itertools import tempfile import os import time from python_test_imports import µ def init_X_F_Chi(res, lens): """ Setup all the needed parameters for initialization of the deformation gradient F and the corresponding deformation map/field Chi_X. Parameters ---------- res: list, int List of grid resoultion in each direction [Nx, Ny, Nz] lens: list, float List of box lengths in each direction [Lx, Ly, Lz] Returns ------- ndim: int Dimension of the structure, derived from len(res). X: np.ndarray, float initial grid point positions. F: np.ndarray initial deformation gradient (right shape, all values set to zero) Chi_X: np.ndarray initial deformation field (right shape, all values set to zero) """ # initialization of parameters res = np.array(res) #grid resolution lens = np.array(lens) #box lengths ndim = len(res) #dimension of the problem X = µ.tools.compute_real_space_lattice(res, lens) F = np.zeros((ndim,) + X.shape) Chi_X = np.zeros(X.shape) return ndim, X, F, Chi_X def nth_order_central_diff_derivative(function_values, grid_spacing, order): """ Compute the first derivative of a function with values 'function values' with the central difference approximation to the order 'order'. The function values are on a rectangualr grid with constant grid spacing 'grid_spacing' for each direction. CAUTION The function is assumed to be periodic (pbc)! Parameters ---------- function_values: numpy.ndarray, float value of the function on an equally spaced grid, with grid spacing 'grid_spacing'. grid_spacing: float or np.array of floats gives the grid spacing distance in each grid direction. For a single value the grid spacing is assumed similar in each direction. order: int, order >= 1 gives the accuracy order of the central difference approximation. This means for first order (order=1), three grid points are used to compute the derivative; order=2, 5 grid points; etc. Returns ------- deriv: np.ndarray, float central difference derivative """ f = function_values d = grid_spacing n = order dim = function_values.shape[0] #dimension of the problem (1D-3D) # Get the stencils/ central difference weights Np = 2*n + 1 #number of gridpoints used weights = sm.central_diff_weights(Np) # Central difference approx. deriv = np.zeros(shape=(dim,)+f.shape) for ax in range(dim): #get derivatives in each direction (x,y,z) for i in range(Np): #get the fuction values on the left and right hand side deriv[:,ax] += weights[i] * np.roll(f, n-i, axis=ax+1) * 1/d[ax] return deriv class MuSpectre_tools_Check(unittest.TestCase): """ Check the implementation of all muSpectre.tools functions. TODO ---- * tests for: -> compute_real_space_lattice() -> compute_reciprocal_space_lattice() -> extend_grid_by_pbc() * better tests for: -> test_point_to_volume_grid() (off diagonal F and nonaff displacements) """ def setUp(self): self.lengths = np.array([2.4, 3.7, 4.1]) self.resolutions = np.array([4, 19, 13]) #tolerance for the norm(Chi_X - integral(F)) self.norm_tol = 1e-8 #tollerance for the central difference approximation test self.norm_tol_cent_diff = 1e-8 #tollerance for i_DFT test self.norm_tol_i_DFT = 1e-8 def test_nth_order_central_diff_derivative(self): """ Test of the nth order central difference approximation of the first derivative of a function on a grid. """ ### 1D res = np.array([50]) lens = np.array([7]) d = lens/res ndim, x, deriv, f = init_X_F_Chi(res, lens) f = np.sin(2*np.pi/lens * x) #f(x)=sin(a*x) approx_deriv = nth_order_central_diff_derivative(f, grid_spacing=d, order=4) deriv = 2*np.pi/lens * np.cos(2*np.pi/lens * x) #f'(x)=a*cos(a*x) self.assertLess(np.linalg.norm(deriv-approx_deriv), self.norm_tol_cent_diff) ### 3D res = np.array([46, 60, 41]) lens = np.array([4, 1.7, 5.3]) d = lens/res ndim, x, deriv, f = init_X_F_Chi(res, lens) for i in range(ndim): deriv[i,i,:,:,:] = 2*np.pi/lens[i] * np.cos(2*np.pi*x[i]/lens[i]) f[i,:,:,:] = np.sin(2*np.pi*x[i]/lens[i]) approx_deriv = nth_order_central_diff_derivative(f, grid_spacing=d, order=5) self.assertLess(np.linalg.norm(deriv-approx_deriv), self.norm_tol_cent_diff) def test_i_DFT(self): """ Test the inverse discrete Fourier transformation (i_DFT) by comparing it with numpy.fft.ifftn() and the given input function. """ for d in range(1,4): res = self.resolutions[:d] lens = self.lengths[:d] ndim = len(res) X = µ.tools.compute_real_space_lattice(res, lens) #real signal f_real = np.sin(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) * X) + 4 f_real_q = np.fft.fftn(f_real, s=res) f_real_numpy = np.fft.ifftn(f_real_q, s=res) f_real_reconstr = µ.tools.i_DFT(f_real_q, res, lens, X) self.assertLess(np.linalg.norm(f_real - f_real_reconstr), self.norm_tol_i_DFT) self.assertLess(np.linalg.norm(f_real_numpy - f_real_reconstr), self.norm_tol_i_DFT) #complex signal comp = np.cos(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) * X) f_comp = np.sin(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) * X) + comp + 2 f_comp_q = np.fft.fftn(f_comp, s=res) f_comp_numpy = np.fft.ifftn(f_comp_q, s=res) f_comp_reconstr = µ.tools.i_DFT(f_comp_q, res, lens, X) self.assertLess(np.linalg.norm(f_comp - f_comp_reconstr), self.norm_tol_i_DFT) self.assertLess(np.linalg.norm(f_comp_numpy - f_comp_reconstr), self.norm_tol_i_DFT) def test_compute_displacements_integration(self): """ Proof if compute_displacements turns a deforamtion gradient F into its corresponding deformation field Chi_X. This is tested for different grid resolutions and cell lengths. One has to keep in minde that the displacement fluctuation field w(x) has to fulfill periodicity conditions of the form: w^-=w^+ on the surfaces delB^- , delB^+ of the periodic region B (see: P. Eisenlohr et al., International Journal of Plasticity 46 (2013), chap. 2.1 and appendix B; The mentioned w here corresponds to w tilde (displacement fluctuation field) in the paper). Additionally the implementation of the correction order is checked. This is tested for 2D and 3D grids. test off diagonal and non-diagonal deformation gradients. """ ###### test for the deformation gradient integration scheme ###### ### F=1I, diagonal deformation gradient res = [8, 15, 13] lens = [4, 2.3, 2.1] ndim, X, F, Chi_X = init_X_F_Chi(res, lens) for i in range(ndim): F[i,i,:,:,:] = 1 Chi_X[i,:,:,:] = X[i] initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### cosine diagonal deformation gradient # The resolution 30 and legth 2.2 in x-direction test a special case of # grid, here X = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1], 0:lens[2]:d[2]] # fails because it makes 31 grid points in x direction. res = [4, 30, 2] lens = [4, 2.2, 2] ndim, X, F, Chi_X = init_X_F_Chi(res, lens) for i in range(ndim): F[i,i,:,:,:] = 2 +0.2*2*np.pi/lens[i]*np.cos(2*np.pi*X[i]/lens[i]) Chi_X[i,:,:,:] = 2*X[i] + 0.2 * np.sin(2*np.pi*X[i]/lens[i]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### polynomial diagonal deformation gradient # Choose the prefactors of the polynomial such that at least Chi_X and F # have respectively the same values at the boundaries (here X_i=0 and # X_i=4). res = [2000, 1, 1] lens = [4, 4, 4] ndim, X, F, Chi_X = init_X_F_Chi(res, lens) for i in range(ndim): F[i,i,:,:,:] = 32*X[i] -24*X[i]**2 +4*X[i]**3 Chi_X[i,:,:,:] = 16*X[i]**2 -8*X[i]**3 +X[i]**4 if i == 0: #You have to subtract the average displacement field. This is #necessary because the Chi_X = F_av*X + w(X) and the boundary #conditions require w^-(X) = w^+(X) Chi_X[i,:,:,:] -= 128/15 initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### non-diagonal deformation gradient # This one needs a lot of grid points, think more than 1000^3, therefore # I should first improve the implementation of compute_diplacements. # Before it takes to much time/crashes #res = [100, 500, 500] #lens = [8, 8, 8] #ndim, X, F, Chi_X = init_X_F_Chi(res, lens) #F[0,0,:,:,:] = 8 -6*X[0] + 3/4*X[0]**2 #F[1,1,:,:,:] = -6*X[1] + 3/4*X[1]**2 #F[2,2,:,:,:] = -6*X[2] + 3/4*X[2]**2 #F[1,0,:,:,:] = 8 #F[2,0,:,:,:] = 8 #for i in range(ndim): # Chi_X[i,:,:,:] = 8*X[0] -3*X[i]**2 +1/4*X[i]**3 # #initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) ### non-diagonal deformation gradient ### 1D ### res = [5] lens = [8] ndim, X, F, Chi_X = init_X_F_Chi(res, lens) F[0,0,:] = 5 + 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) Chi_X[0,:] = 5*X[0] + 2*np.sin(2*np.pi*X[0]/lens[0]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### 2D ### res = [5, 5] lens = [8, 8] ndim, X, F, Chi_X = init_X_F_Chi(res, lens) F[0,0,:,:] = 5 + 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) F[1,1,:,:] = 5 + 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1]) F[1,0,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) for i in range(ndim): Chi_X[i,:,:] = 5*X[i] + np.sin(2*np.pi*X[i]/lens[i]) \ + np.sin(2*np.pi*X[0]/lens[0]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### 3D ### res = [5, 5, 5] lens = [8, 8, 8] ndim, X, F, Chi_X = init_X_F_Chi(res, lens) F[0,0,:,:,:] = 5 + 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) F[1,1,:,:,:] = 5 + 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1]) F[2,2,:,:,:] = 5 + 2*np.pi/lens[2]*np.cos(2*np.pi/lens[2]*X[2]) F[1,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) F[2,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) for i in range(ndim): Chi_X[i,:,:,:] = 5*X[i] + np.sin(2*np.pi*X[i]/lens[i]) \ + np.sin(2*np.pi*X[0]/lens[0]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ###### tests for higher order corrections ###### order_o = [1, 2] ### cosine diagonal deformation gradient res = [6, 10] lens = [7, 1.4] d = np.array(lens)/np.array(res) ndim, X, F, Chi_X = init_X_F_Chi(res, lens) for i in range(ndim): F[i,i,:,:] = 0.2*2*np.pi/lens[i] * np.cos(2*np.pi*X[i]/lens[i]) Chi_X[i,:,:] = 0.2 * np.sin(2*np.pi*X[i]/lens[i]) # zeroth order correction initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) # higher order correction for n in order_o: F = nth_order_central_diff_derivative(Chi_X, d, order=n) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### polynomial diagonal deformation gradient # Choose the prefactors of the polynomial such that at least Chi_X and F # have respectively the same values at the boundaries (here X_i=0 and # X_i=4). res = [2000, 1, 1] lens = [4, 4, 4] d = np.array(lens)/np.array(res) ndim, X, F, Chi_X = init_X_F_Chi(res, lens) for i in range(ndim): F[i,i,:,:,:] = 32*X[i] -24*X[i]**2 +4*X[i]**3 Chi_X[i,:,:,:] = 16*X[i]**2 -8*X[i]**3 +X[i]**4 if i == 0: #You have to subtract the average displacement field. This is #necessary because the Chi_X = F_av*X + w(X) and the boundary #conditions require w^-(X) = w^+(X) Chi_X[i,:,:,:] -= 128/15 # zeroth order correction initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) # higher order correction for n in order_o: F = nth_order_central_diff_derivative(Chi_X, d, order=n) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### non-diagonal deformation gradient ### 2D ### res = [5, 5] lens = [8, 8] d = np.array(lens)/np.array(res) ndim, X, F, Chi_X = init_X_F_Chi(res, lens) F[0,0,:,:] = 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) F[1,1,:,:] = 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1]) F[1,0,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) for i in range(ndim): Chi_X[i,:,:] = np.sin(2*np.pi*X[i]/lens[i]) \ + np.sin(2*np.pi*X[0]/lens[0]) # zeroth order correction initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) # higher order correction for n in order_o: F = nth_order_central_diff_derivative(Chi_X, d, order=n) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### 3D ### res = [5, 5, 5] lens = [8, 8, 8] d = np.array(lens)/np.array(res) ndim, X, F, Chi_X = init_X_F_Chi(res, lens) F[0,0,:,:,:] = 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) F[1,1,:,:,:] = 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1]) F[2,2,:,:,:] = 2*np.pi/lens[2]*np.cos(2*np.pi/lens[2]*X[2]) F[1,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) F[2,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) for i in range(ndim): Chi_X[i,:,:,:] = np.sin(2*np.pi*X[i]/lens[i]) \ + np.sin(2*np.pi*X[0]/lens[0]) # zeroth order correction initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) # higher order correction for n in order_o: F = nth_order_central_diff_derivative(Chi_X, d, order=n) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### real test, shear of material with two different hardnesses/ Young moduli. ### Can I see Gibbs ringing? #initialize res = [101, 101, 1] lens = [10, 10, 1] ndim = len(res) #dimension of the problem Nx, Ny, Nz = res Lx, Ly, Lz = lens gs = np.array(lens) / np.array(res) #grid spacing formulation = µ.Formulation.small_strain Young = [10, 20] #Youngs modulus for each phase (soft,hard) Poisson = [0.3, 0.3] #Poissons ratio for each phase order = 0 #correction order of Fourier integration #solver newton_tol = 1e-6 #tolerance for newton algo cg_tol = 1e-6 #tolerance for cg algo equil_tol = 1e-6 #tolerance for equilibrium maxiter = 100 verbose = 0 ### geometry (two slabs above each other in y-direction) phase = np.zeros((Nx,Ny,Nz), dtype=int) #0 for surrounding matrix boundary = res[1]//2 phase[:, boundary:, :] = 1 cell = µ.Cell(res, lens, formulation) mat = µ.material.MaterialLinearElastic4_3d.make(cell, "material") for i, pixel in enumerate(cell): #add Young and Poisson depending on the material index m_i = phase.flatten()[i] #m_i = material index (phase index) mat.add_pixel(pixel, Young[m_i], Poisson[m_i]) cell.initialise() DelF = np.array([[0 , 0.1, 0], [0.1, 0 , 0], [0 , 0 , 0]]) solver_newton = µ.solvers.SolverCG(cell, cg_tol, maxiter, verbose) r_newton = µ.solvers.newton_cg(cell, DelF, solver_newton, newton_tol, equil_tol, verbose) ### check strains strain = r_newton.grad.T.reshape(*res, 3, 3).transpose((4,3,0,1,2)) one = np.eye(ndim).reshape((ndim,)*2+(1,)*ndim) F = strain + one initial_pos, displacements = \ µ.tools.compute_displacements(F, res, lens, order) fin_pos = initial_pos + displacements ### analytic solution # Why L - gs / 2 ?! # The analytic signal doesn't know about the periodic boundary conditions, # thus the analytic box is L_ana = L-d; and the Heaviside step is inbetween two # grid points. # # solution : fourier analytic # Grid points : + + + # # # o + + + # # # # Step signal : ---------\ /- ----------| # \ / | # \---------/ |---------- # Length : |<--------- L --------->| |<- L_ana = L-d ->| # soft/hard : |< L_soft >|< L_hard >| | L_soft | L_hard | # grid spacing: |gs| l_soft_ana = gs[1]*(boundary - 1/2) #(Ly - gs[1]) / 2 , analytic height l_hard_ana = gs[1]*(Ny - boundary - 1/2) #(Ly - gs[1]) / 2 , analytic height l_soft = gs[1] * boundary #height material soft l_hard = gs[1] * (Ny-boundary) #height material hard mu = np.array(Young) / (2 * (1+np.array(Poisson))) #shear modul (mu_soft, mu_hard) gamma = 2*DelF[0,1] #mean deformation g_soft = (Ly*gamma) / (l_soft + l_hard * mu[0]/mu[1]) g_hard = (Ly*gamma) / (l_soft * mu[1]/mu[0] + l_hard) ###compute d (final positions) # initial grid X = µ.tools.compute_real_space_lattice(res, lens) d = np.zeros(X.shape)#X.copy() #initialize new positions #soft d[0, :, :boundary, :] += g_soft/2 * X[1, :, :boundary, :] #d_soft x coordinate d[1, :, :boundary, :] += gamma/2 * X[0, :, :boundary, :] #d_soft y coordinate #hard d[0, :, boundary:, :] += g_hard/2 * (X[1, :, boundary:, :]-l_soft_ana) \ + g_soft/2 * l_soft_ana d[1, :, boundary:, :] += gamma/2 * (X[0, :, boundary:, :]) #correct by integration constant # shift the analytic solution such that: # -> the average nonaffine deformation is zero (integral of the nonaffine # deformation gradient + integration const). # ==> choose the right integration const dim_axis = tuple(-np.arange(1,ndim+1)) #the last ndim entries F_homo = 1./(np.prod(res)) * F.sum(axis=dim_axis) F_homo = F_homo.reshape((ndim,)*2+(1,)*ndim) # integration constant = - (integral / number of gridpoints) int_const = - ((d[0] - F_homo[0,1] * X[1]).sum(axis=1))[0,0] / Ny #final positions = displacements + old positions + integration const. d += X + np.array([0, int_const, 0]).reshape((ndim,)+(1,)*ndim) #self.assertLess(np.linalg.norm(fin_pos - d), self.norm_tol) ### real test, integrate a delta function deformation gradient to a step function. ### Can I see Gibbs ringing? # does this test make sense? -> discuss with Till! def test_point_to_volume_grid(self): """ Test the implementation of the function point_to_volume_grid() TODO ---- *check affine displacements/and first implement them right... --> Do I need the nonaffine_displacement correction any more???!!! """ st = time.time() lens = self.lengths res = self.resolutions ndim = len(res) d = lens/res point_grid = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1], 0:lens[2]:d[2]] # for an orthogonal same spaced point_grid one can easily compute the # volume_grid by just increase (pbc=False) the mesh by one in each # direction and choose the points inbetween (d/2) the point_grid # positions. vg = np.mgrid[0-d[0]/2:lens[0]:d[0], 0-d[1]/2:lens[1]:d[1], 0-d[2]/2:lens[2]:d[2]] vg_pbc = np.mgrid[0-d[0]/2:lens[0]-d[0]/2:d[0], 0-d[1]/2:lens[1]-d[1]/2:d[1], 0-d[2]/2:lens[2]-d[2]/2:d[2]] # loop over 1D to 3D, for pbc=False/True for dim in range(1,4): s = list(np.s_[:dim,:,:,:]) try: s[dim+1] = 0 s[dim+2] = 0 except IndexError: pass s = tuple(s) F_eye = np.einsum('ij,k...->ijk...',np.eye(dim),np.ones(res[:dim])) calc_vg = µ.tools.point_to_volume_grid( point_grid = point_grid[s], resolutions = res[:dim], lengths = lens[:dim], F = F_eye, pbc = False, nonaffine_displacements = np.array(None) ) self.assertLess(np.linalg.norm(calc_vg - vg[s]), self.norm_tol) calc_vg_pbc = µ.tools.point_to_volume_grid( point_grid = point_grid[s], resolutions = res[:dim], lengths = lens[:dim], F = F_eye, pbc = True, nonaffine_displacements = np.array(None) ) self.assertLess(np.linalg.norm(calc_vg_pbc - vg_pbc[s]), self.norm_tol ) et = time.time() - print('By hand interpol {}: {} s'.format(res[:dim], et-st)) + print('In tools.test_point_to_volume_grid()\n' + ,'{}s , for a {} grid.\n'.format(et-st, res[:dim])) def test_point_to_volume_grid_fourier(self): """ Test the implementation of the function point_to_volume_grid_fourier(). Should do the same as point_to_volume_grid(). """ st = time.time() lens = self.lengths res = self.resolutions ndim = len(res) d = lens/res point_grid = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1], 0:lens[2]:d[2]] # for an orthogonal same spaced point_grid one can easily compute the # volume_grid by just increase (pbc=False) the mesh by one in each # direction and choose the points inbetween (d/2) the point_grid # positions. vg = np.mgrid[0-d[0]/2:lens[0]:d[0], 0-d[1]/2:lens[1]:d[1], 0-d[2]/2:lens[2]:d[2]] vg_pbc = np.mgrid[0-d[0]/2:lens[0]-d[0]/2:d[0], 0-d[1]/2:lens[1]-d[1]/2:d[1], 0-d[2]/2:lens[2]-d[2]/2:d[2]] # loop over 1D to 3D, for pbc=False/True for dim in range(1,4): s = list(np.s_[:dim,:,:,:]) try: s[dim+1] = 0 s[dim+2] = 0 except IndexError: pass s = tuple(s) F_eye = np.einsum('ij,k...->ijk...',np.eye(dim),np.ones(res[:dim])) - #print('dimension: ', dim) - #print('res: ', res[:dim]) calc_vg = µ.tools.point_to_volume_grid_fourier( point_grid = point_grid[s], resolutions = res[:dim], lengths = lens[:dim], F = F_eye, pbc = False ) self.assertLess(np.linalg.norm(calc_vg - vg[s]), self.norm_tol) calc_vg_pbc = µ.tools.point_to_volume_grid_fourier( point_grid = point_grid[s], resolutions = res[:dim], lengths = lens[:dim], F = F_eye, pbc = True ) self.assertLess(np.linalg.norm(calc_vg_pbc - vg_pbc[s]), self.norm_tol ) et = time.time() - print('Fourier interpol {}: {} s'.format(res[:dim], et-st)) + print('In tools.test_point_to_volume_grid_fourier()\n' + ,'{}s , for a {} grid.\n'.format(et-st, res[:dim])) def test_write_structure_as_vtk(self): """ Test the implementation of the function write_structure_as_vtk(). It will write 1D, 2D and 3D structures in a folder which will be deleted after the files were selfuccessfully written. TODO ---- Check not only if the files were written but also their content! """ #initialize some structure lens = self.lengths res = self.resolutions pos = {} #positions dictionary cd = {} #cell data dictionary pd = {} #point data dictionary for i in range(1,4): name = "{}d".format(i) pos[name] = init_X_F_Chi(res[:i], lens[:i])[1] cd_shape = tuple(np.array(pos[name].shape[1:])-1) + (1,)*(3-i) pd_shape = pos[name].shape[1:] + (1,)*(3-i) cd[name] = {'random values' : np.random.random_sample(cd_shape)} pd[name] = {'random values' : np.random.random_sample(pd_shape)} #create a temporary directory and run the tests in it dir_name = 'write_structure_as_vtk_tests' cwd_start = os.getcwd() #The temporary directory is atomatically cleand up after one is exiting #the with block. with tempfile.TemporaryDirectory(dir=cwd_start) as dir_name: #move to tmp dir os.chdir(dir_name) #start with the test (write file and check if they were written) for i in range(1,4): name = "{}d".format(i) µ.tools.write_structure_as_vtk('{}'.format(name), pos[name], None, None) µ.tools.write_structure_as_vtk('{}_cdpd'.format(name), pos[name], cd[name], pd[name]) #test if files were written error_message = '\nThere is a problem in the function ' + \ 'muSpectre.tools.write_structure_as_vtk().\n' +\ 'The file {} was not written.' assert os.path.exists('{}.vts'.format(name)) == 1, \ error_message.format(name+'_grid.vts') assert os.path.exists('{}_cdpd.vts'.format(name)) == 1, \ error_message.format(name+'_cdpd_grid.vts') os.chdir(cwd_start) if __name__ == '__main__': unittest.main()