diff --git a/src/materials/material_hyper_elasto_plastic1.cc b/src/materials/material_hyper_elasto_plastic1.cc index 7e71297..3f9dda0 100644 --- a/src/materials/material_hyper_elasto_plastic1.cc +++ b/src/materials/material_hyper_elasto_plastic1.cc @@ -1,64 +1,66 @@ /** * @file material_hyper_elasto_plastic1.cc * * @author Till Junge * * @date 21 Feb 2018 * * @brief implementation for MaterialHyperElastoPlastic1 * * 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_plastic1.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialHyperElastoPlastic1:: MaterialHyperElastoPlastic1(std::string name, Real young, Real poisson, Real tau_y0, Real H) : 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), young{young}, poisson{poisson}, lambda{Hooke::compute_lambda(young, poisson)}, mu{Hooke::compute_mu(young, poisson)}, K{Hooke::compute_K(young, poisson)}, tau_y0{tau_y0}, H{H}, - C{Hooke::compute_C_T4(lambda, mu)}, + // the factor .5 comes from equation (18) in Geers 2003 + // (https://doi.org/10.1016/j.cma.2003.07.014) + C{0.5*Hooke::compute_C_T4(lambda, mu)}, internal_variables{F_prev_field.get_map(), be_prev_field.get_map(), plast_flow_field.get_map()} {} /* ---------------------------------------------------------------------- */ template void MaterialHyperElastoPlastic1::save_history_variables() { this->plast_flow_field.cycle(); this->F_prev_field.cycle(); this->be_prev_field.cycle(); } template class MaterialHyperElastoPlastic1< twoD, twoD>; template class MaterialHyperElastoPlastic1< twoD, threeD>; template class MaterialHyperElastoPlastic1; } // muSpectre diff --git a/src/materials/material_hyper_elasto_plastic1.hh b/src/materials/material_hyper_elasto_plastic1.hh index b8ba661..2ceaef9 100644 --- a/src/materials/material_hyper_elasto_plastic1.hh +++ b/src/materials/material_hyper_elasto_plastic1.hh @@ -1,301 +1,301 @@ /** * @file material_hyper_elasto_plastic1.hh * * @author Till Junge * * @date 20 Feb 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) * * 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_PLASTIC1_H #define MATERIAL_HYPER_ELASTO_PLASTIC1_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 MaterialHyperElastoPlastic1; /** * 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>; /** * 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 MaterialHyperElastoPlastic1: 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 MaterialHyperElastoPlastic1() = delete; //! Constructor with name and material properties MaterialHyperElastoPlastic1(std::string name, Real young, Real poisson, Real tau_y0, Real H); //! Copy constructor MaterialHyperElastoPlastic1(const MaterialHyperElastoPlastic1 &other) = delete; //! Move constructor MaterialHyperElastoPlastic1(MaterialHyperElastoPlastic1 &&other) = delete; //! Destructor virtual ~MaterialHyperElastoPlastic1() = default; //! Copy assignment operator MaterialHyperElastoPlastic1& operator=(const MaterialHyperElastoPlastic1 &other) = delete; //! Move assignment operator MaterialHyperElastoPlastic1& operator=(MaterialHyperElastoPlastic1 &&other) = delete; /** * evaluates Kirchhoff stress given the current placement gradient * Fₜ, the previous Gradient Fₜ₋₁ and the cumulated plastic flow * εₚ */ template inline decltype(auto) evaluate_stress(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t plast_flow); /** * evaluates Kirchhoff stress and stiffness given the current placement gradient * Fₜ, the previous Gradient Fₜ₋₁ and the cumulated plastic flow * εₚ */ template inline decltype(auto) evaluate_stress_tangent(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t plast_flow); /** * The statefields need to be cycled at the end of each load increment */ virtual void save_history_variables() override; /** * return the internals tuple */ typename traits::InternalVariables & get_internals() { return this->internal_variables;}; 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); using LColl_t = LocalFieldCollection; //! storage for cumulated plastic flow εₚ StateField> plast_flow_field; //! storage for previous gradient Fᵗ StateField> F_prev_field; //! storage for elastic left Cauchy-Green deformation tensor bₑᵗ StateField> be_prev_field; // material properties const Real young; //!< Young's modulus const Real poisson; //!< Poisson's ratio const Real lambda; //!< first Lamé constant const Real mu; //!< second Lamé constant (shear modulus) const Real K; //!< Bulk modulus const Real tau_y0; //!< initial yield stress const Real H; //!< hardening modulus const T4Mat C; //!< stiffness tensor typename traits::InternalVariables internal_variables; private: }; //----------------------------------------------------------------------------// template template decltype(auto) MaterialHyperElastoPlastic1:: stress_n_internals_worker(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t eps_p) { // 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()}; - auto && ln_be_star{logm(std::move(be_star))}; - auto && tau_star{.5*Hooke::evaluate_stress(this->lambda, this->mu, ln_be_star)}; + Mat_t ln_be_star{logm(std::move(be_star))}; + Mat_t tau_star{.5*Hooke::evaluate_stress(this->lambda, this->mu, ln_be_star)}; // deviatoric part of Kirchhoff stress - auto && tau_d_star{tau_star - tau_star.trace()/DimM*tau_star.Identity()}; + 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())}; - auto && N_star{3*.5*tau_d_star/tau_eq_star}; + 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) Real phi_star{std::max(tau_eq_star - this->tau_y0 - this->H * eps_p.old(), 0.)}; // return mapping - auto && Del_gamma{phi_star/(this->H + 3 * this->mu)}; + Real Del_gamma{phi_star/(this->H + 3 * this->mu)}; auto && tau{tau_star - 2*Del_gamma*this->mu*N_star}; - //auto && tau_eq{tau_eq_star - 3*this->mu*Del_gamma}; + /////auto && tau_eq{tau_eq_star - 3*this->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}; 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) MaterialHyperElastoPlastic1:: evaluate_stress(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t eps_p) { auto retval(std::move(std::get<0>(this->stress_n_internals_worker (std::forward(F), F_prev, be_prev, eps_p)))); return retval; } //----------------------------------------------------------------------------// template template decltype(auto) MaterialHyperElastoPlastic1:: evaluate_stress_tangent(grad_t && F, StrainStRef_t F_prev, StrainStRef_t be_prev, FlowStRef_t eps_p) { //! 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)}; 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*this->mu/tau_eq_star; auto && a1 = this->mu/(this->H + 3*this->mu); return std::make_tuple(std::move(tau), T4Mat{ ((this->K/2. - this->mu/3 + a0*this->mu)*Matrices::Itrac() + (1 - 3*a0) * this->mu*Matrices::Isymm() + 2*this->mu * (a0 - a1)*Matrices::outer(N_star, N_star))}); } else { return std::make_tuple(std::move(tau), T4Mat{this->C}); } } } // muSpectre #endif /* MATERIAL_HYPER_ELASTO_PLASTIC1_H */ diff --git a/tests/elasto-plasticity.py b/tests/elasto-plasticity.py index dcb4f2d..f7ad35e 100644 --- a/tests/elasto-plasticity.py +++ b/tests/elasto-plasticity.py @@ -1,282 +1,287 @@ import numpy as np import scipy.sparse.linalg as sp import itertools import sys # turn of warning for zero division # (which occurs in the linearization of the logarithmic strain) np.seterr(divide='ignore', invalid='ignore') # ----------------------------------- GRID ------------------------------------ Nx = 1 # number of voxels in x-direction Ny = 1 # number of voxels in y-direction Nz = 1 # number of voxels in z-direction shape = [Nx,Ny,Nz] # number of voxels as list: [Nx,Ny,Nz] # ----------------------------- TENSOR OPERATIONS ----------------------------- # tensor operations / products: np.einsum enables index notation, avoiding loops # e.g. ddot42 performs $C_ij = A_ijkl B_lk$ for the entire grid trans2 = lambda A2 : np.einsum('ijxyz ->jixyz ',A2 ) ddot22 = lambda A2,B2: np.einsum('ijxyz ,jixyz ->xyz ',A2,B2) ddot42 = lambda A4,B2: np.einsum('ijklxyz,lkxyz ->ijxyz ',A4,B2) ddot44 = lambda A4,B4: np.einsum('ijklxyz,lkmnxyz->ijmnxyz',A4,B4) dot11 = lambda A1,B1: np.einsum('ixyz ,ixyz ->xyz ',A1,B1) dot22 = lambda A2,B2: np.einsum('ijxyz ,jkxyz ->ikxyz ',A2,B2) dot24 = lambda A2,B4: np.einsum('ijxyz ,jkmnxyz->ikmnxyz',A2,B4) dot42 = lambda A4,B2: np.einsum('ijklxyz,lmxyz ->ijkmxyz',A4,B2) dyad22 = lambda A2,B2: np.einsum('ijxyz ,klxyz ->ijklxyz',A2,B2) dyad11 = lambda A1,B1: np.einsum('ixyz ,jxyz ->ijxyz ',A1,B1) # eigenvalue decomposition of 2nd-order tensor: return in convention i,j,x,y,z # NB requires to swap default order of NumPy (in in/output) def eig2(A2): swap1i = lambda A1: np.einsum('xyzi ->ixyz ',A1) swap2 = lambda A2: np.einsum('ijxyz->xyzij',A2) swap2i = lambda A2: np.einsum('xyzij->ijxyz',A2) vals,vecs = np.linalg.eig(swap2(A2)) vals = swap1i(vals) vecs = swap2i(vecs) return vals,vecs # logarithm of grid of 2nd-order tensors def ln2(A2): vals,vecs = eig2(A2) return sum([np.log(vals[i])*dyad11(vecs[:,i],vecs[:,i]) for i in range(3)]) # exponent of grid of 2nd-order tensors def exp2(A2): vals,vecs = eig2(A2) return sum([np.exp(vals[i])*dyad11(vecs[:,i],vecs[:,i]) for i in range(3)]) # determinant of grid of 2nd-order tensors def det2(A2): return (A2[0,0]*A2[1,1]*A2[2,2]+A2[0,1]*A2[1,2]*A2[2,0]+A2[0,2]*A2[1,0]*A2[2,1])-\ (A2[0,2]*A2[1,1]*A2[2,0]+A2[0,1]*A2[1,0]*A2[2,2]+A2[0,0]*A2[1,2]*A2[2,1]) # inverse of grid of 2nd-order tensors def inv2(A2): A2det = det2(A2) A2inv = np.empty([3,3,Nx,Ny,Nz]) A2inv[0,0] = (A2[1,1]*A2[2,2]-A2[1,2]*A2[2,1])/A2det A2inv[0,1] = (A2[0,2]*A2[2,1]-A2[0,1]*A2[2,2])/A2det A2inv[0,2] = (A2[0,1]*A2[1,2]-A2[0,2]*A2[1,1])/A2det A2inv[1,0] = (A2[1,2]*A2[2,0]-A2[1,0]*A2[2,2])/A2det A2inv[1,1] = (A2[0,0]*A2[2,2]-A2[0,2]*A2[2,0])/A2det A2inv[1,2] = (A2[0,2]*A2[1,0]-A2[0,0]*A2[1,2])/A2det A2inv[2,0] = (A2[1,0]*A2[2,1]-A2[1,1]*A2[2,0])/A2det A2inv[2,1] = (A2[0,1]*A2[2,0]-A2[0,0]*A2[2,1])/A2det A2inv[2,2] = (A2[0,0]*A2[1,1]-A2[0,1]*A2[1,0])/A2det return A2inv # ------------------------ INITIATE (IDENTITY) TENSORS ------------------------ # identity tensor (single tensor) i = np.eye(3) # identity tensors (grid) I = np.einsum('ij,xyz' , i ,np.ones([Nx,Ny,Nz])) I4 = np.einsum('ijkl,xyz->ijklxyz',np.einsum('il,jk',i,i),np.ones([Nx,Ny,Nz])) I4rt = np.einsum('ijkl,xyz->ijklxyz',np.einsum('ik,jl',i,i),np.ones([Nx,Ny,Nz])) I4s = (I4+I4rt)/2. II = dyad22(I,I) # ------------------------------------ FFT ------------------------------------ # projection operator (only for non-zero frequency, associated with the mean) # NB: vectorized version of "hyper-elasticity.py" # - allocate / support function Ghat4 = np.zeros([3,3,3,3,Nx,Ny,Nz]) # projection operator x = np.zeros([3 ,Nx,Ny,Nz],dtype='int64') # position vectors q = np.zeros([3 ,Nx,Ny,Nz],dtype='int64') # frequency vectors delta = lambda i,j: np.float(i==j) # Dirac delta function # - set "x" as position vector of all grid-points [grid of vector-components] x[0],x[1],x[2] = np.mgrid[:Nx,:Ny,:Nz] # - convert positions "x" to frequencies "q" [grid of vector-components] for i in range(3): freq = np.arange(-(shape[i]-1)/2,+(shape[i]+1)/2,dtype='int64') q[i] = freq[x[i]] # - compute "Q = ||q||", and "norm = 1/Q" being zero for the mean (Q==0) # NB: avoid zero division q = q.astype(np.float) Q = dot11(q,q) Z = Q==0 Q[Z] = 1. norm = 1./Q norm[Z] = 0. # - set projection operator [grid of tensors] for i, j, l, m in itertools.product(range(3), repeat=4): Ghat4[i,j,l,m] = norm*delta(i,m)*q[j]*q[l] # (inverse) Fourier transform (for each tensor component in each direction) fft = lambda x: np.fft.fftshift(np.fft.fftn (np.fft.ifftshift(x),[Nx,Ny,Nz])) ifft = lambda x: np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(x),[Nx,Ny,Nz])) # functions for the projection 'G', and the product 'G : K^LT : (delta F)^T' G = lambda A2 : np.real( ifft( ddot42(Ghat4,fft(A2)) ) ).reshape(-1) K_dF = lambda dFm: trans2(ddot42(K4,trans2(dFm.reshape(3,3,Nx,Ny,Nz)))) G_K_dF = lambda dFm: G(K_dF(dFm)) # --------------------------- CONSTITUTIVE RESPONSE --------------------------- # constitutive response to a certain loading and history # NB: completely uncoupled from the FFT-solver, but implemented as a regular # grid of quadrature points, to have an efficient code; # each point is completely independent, just evaluated at the same time def constitutive(F,F_t,be_t,ep_t): # function to compute linearization of the logarithmic Finger tensor def dln2_d2(A2): vals,vecs = eig2(A2) K4 = np.zeros([3,3,3,3,Nx,Ny,Nz]) for m, n in itertools.product(range(3),repeat=2): gc = (np.log(vals[n])-np.log(vals[m]))/(vals[n]-vals[m]) gc[vals[n]==vals[m]] = (1.0/vals[m])[vals[n]==vals[m]] K4 += gc*dyad22(dyad11(vecs[:,m],vecs[:,n]),dyad11(vecs[:,m],vecs[:,n])) return K4 # elastic stiffness tensor C4e = K*II+2.*mu*(I4s-1./3.*II) # trial state Fdelta = dot22(F,inv2(F_t)) be_s = dot22(Fdelta,dot22(be_t,trans2(Fdelta))) lnbe_s = ln2(be_s) tau_s = ddot42(C4e,lnbe_s)/2. taum_s = ddot22(tau_s,I)/3. taud_s = tau_s-taum_s*I print("taud_s =\n{}".format(taud_s.reshape([3, 3]))) taueq_s = np.sqrt(3./2.*ddot22(taud_s,taud_s)) print("taueq_s = {}".format(taueq_s)) N_s = 3./2.*taud_s/taueq_s + print("N_s =\n{}".format(N_s.reshape([3, 3]))) + print("tauy0 = {}".format(tauy0)) phi_s = taueq_s-(tauy0+H*ep_t) - print("phi_s = {}".format(phi_s)) + print("phi_s(before) = {}".format(phi_s)) print("ep_t = {}".format(ep_t)) phi_s = 1./2.*(phi_s+np.abs(phi_s)) print("phi_s = {}".format(phi_s)) # return map dgamma = phi_s/(H+3.*mu) print("dgamma = {}".format(dgamma)) ep = ep_t + dgamma tau = tau_s -2.*dgamma*N_s*mu lnbe = lnbe_s-2.*dgamma*N_s be = exp2(lnbe) P = dot22(tau,trans2(inv2(F))) # consistent tangent operator a0 = dgamma*mu/taueq_s a1 = mu/(H+3.*mu) C4ep = ((K-2./3.*mu)/2.+a0*mu)*II+(1.-3.*a0)*mu*I4s+2.*mu*(a0-a1)*dyad22(N_s,N_s) dlnbe4_s = dln2_d2(be_s) dbe4_s = 2.*dot42(I4s,be_s) K4 = (C4e/2.)*(phi_s<=0.).astype(np.float)+C4ep*(phi_s>0.).astype(np.float) + C4 = K4.copy() K4 = ddot44(K4,ddot44(dlnbe4_s,dbe4_s)) K4 = dot42(-I4rt,tau)+K4 K4 = dot42(dot24(inv2(F),K4),trans2(inv2(F))) - return tau, P,K4,be,ep + return tau, C4, P,K4,be,ep # phase indicator: square inclusion of volume fraction (3*3*15)/(11*13*15) phase = np.zeros([Nx,Ny,Nz]); phase[:3,:3,:] = 1. # function to convert material parameters to grid of scalars param = lambda M0,M1: M0*np.ones([Nx,Ny,Nz])*(1.-phase)+\ M1*np.ones([Nx,Ny,Nz])* phase # material parameters K = param(0.833,0.833) # bulk modulus mu = param(0.386,0.386) # shear modulus -H = param(0.004,0.008) # hardening modulus -tauy0 = param(0.003,0.006) # initial yield stress +H = param(0.004,0.004) # hardening modulus +tauy0 = param(0.003,0.003) # initial yield stress # ---------------------------------- LOADING ---------------------------------- # stress, deformation gradient, plastic strain, elastic Finger tensor # NB "_t" signifies that it concerns the value at the previous increment ep_t = np.zeros([ Nx,Ny,Nz]) P = np.zeros([3,3,Nx,Ny,Nz]) F = np.array(I,copy=True) F_t = np.array(I,copy=True) be_t = np.array(I,copy=True) # initialize macroscopic incremental loading ninc = 50 lam = 0.0 barF = np.array(I,copy=True) barF_t = np.array(I,copy=True) # initial tangent operator: the elastic tangent K4 = K*II+2.*mu*(I4s-1./3.*II) F = np.array(I,copy=True) F[0, 1] = 1.e-5 -tau, P,K4,be,ep = constitutive(F, F_t, be_t, ep_t) +tau, C4, P,K4,be,ep = constitutive(F, F_t, be_t, ep_t) # end-of-increment: update history barF_t = np.array(barF,copy=True) F_t = np.array(F ,copy=True) be_t = np.array(be ,copy=True) ep_t = np.array(ep ,copy=True) print("tau1 =\n{}".format(tau.reshape([3, 3]))) print("F_t1 =\n{}".format(F_t.reshape([3, 3]))) print("be_t1 =\n{}".format(be_t.reshape([3, 3]))) print("ep_t1 =\n{}".format(ep_t.reshape([1]))) +print("C4 =\n{}".format(C4.reshape([9, 9]))) F = np.array(I,copy=True) F[0, 1] = .2 -tau, P,K4,be,ep = constitutive(F, F_t, be_t, ep_t) +tau, C4, P,K4,be,ep = constitutive(F, F_t, be_t, ep_t) print("tau2 =\n{}".format(tau.reshape([3, 3]))) print("tau2 =\n{}".format(tau.reshape([3, 3]))) -print("F_t2 =\n{}".format(F_t.reshape([3, 3]))) -print("be_t2 =\n{}".format(be_t.reshape([3, 3]))) -print("ep_t2 =\n{}".format(ep_t.reshape([1]))) +print("F_t2 =\n{}".format(F.reshape([3, 3]))) +print("be_t2 =\n{}".format(be.reshape([3, 3]))) +print("ep_t2 =\n{}".format(ep.reshape([1]))) +print("C4 =\n{}".format(C4.reshape([9, 9]))) sys.exit() # incremental deformation for inc in range(1,ninc): print('=============================') print('inc: {0:d}'.format(inc)) # set macroscopic deformation gradient (pure-shear) lam += 0.2/float(ninc) barF = np.array(I,copy=True) barF[0,0] = (1.+lam) barF[1,1] = 1./(1.+lam) # store normalization Fn = np.linalg.norm(F) # first iteration residual: distribute "barF" over grid using "K4" b = -G_K_dF(barF-barF_t) F += barF-barF_t # parameters for Newton iterations: normalization and iteration counter Fn = np.linalg.norm(F) iiter = 0 # iterate as long as the iterative update does not vanish while True: # solve linear system using the Conjugate Gradient iterative solver dFm,_ = sp.cg(tol=1.e-8, A = sp.LinearOperator(shape=(F.size,F.size),matvec=G_K_dF,dtype='float'), b = b, ) # add solution of linear system to DOFs F += dFm.reshape(3,3,Nx,Ny,Nz) # compute residual stress and tangent, convert to residual P,K4,be,ep = constitutive(F,F_t,be_t,ep_t) b = -G(P) # check for convergence, print convergence info to screen print('{0:10.2e}'.format(np.linalg.norm(dFm)/Fn)) if np.linalg.norm(dFm)/Fn<1.e-5 and iiter>0: break # update Newton iteration counter iiter += 1 # end-of-increment: update history barF_t = np.array(barF,copy=True) F_t = np.array(F ,copy=True) be_t = np.array(be ,copy=True) ep_t = np.array(ep ,copy=True) diff --git a/tests/test_material_hyper_elasto_plastic1.cc b/tests/test_material_hyper_elasto_plastic1.cc index 078a326..4cc0ef5 100644 --- a/tests/test_material_hyper_elasto_plastic1.cc +++ b/tests/test_material_hyper_elasto_plastic1.cc @@ -1,147 +1,343 @@ /** * @file test_material_hyper_elasto_plastic1.cc * * @author Till Junge * * @date 25 Feb 2018 * * @brief Tests for the large-strain Simo-type plastic law implemented * using MaterialMuSpectre * * 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 "boost/mpl/list.hpp" #include "materials/material_hyper_elasto_plastic1.hh" #include "materials/materials_toolbox.hh" #include "tests.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_hyper_elasto_plastic_1); template struct MaterialFixture { using Mat = Mat_t; constexpr static Real K{.833}; // bulk modulus constexpr static Real mu{.386}; // shear modulus constexpr static Real H{.004}; // hardening modulus constexpr static Real tau_y0{.003}; // initial yield stress constexpr static Real young{ MatTB::convert_elastic_modulus(K, mu)}; constexpr static Real poisson{ MatTB::convert_elastic_modulus(K, mu)}; MaterialFixture():mat("Name", young, poisson, tau_y0, H){}; constexpr static Dim_t sdim{Mat_t::sdim()}; constexpr static Dim_t mdim{Mat_t::mdim()}; Mat_t mat; }; using mats = boost::mpl::list>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mats, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_stress, Fix, mats, Fix) { + + // This test uses precomputed reference values (computed using + // elasto-plasticity.py) for the 3d case only + + // need higher tol because of printout precision of reference solutions + constexpr Real hi_tol{1e-8}; constexpr Dim_t mdim{Fix::mdim}, sdim{Fix::sdim}; - constexpr bool verbose{true}; + constexpr bool has_precomputed_values{(mdim == sdim) && (mdim == threeD)}; + constexpr bool verbose{false}; using Strain_t = Eigen::Matrix; using traits = MaterialMuSpectre_traits>; using LColl_t = typename traits::LFieldColl_t; using StrainStField_t = StateField< TensorField>; using FlowStField_t = StateField< ScalarField>; // using StrainStRef_t = typename traits::LStrainMap_t::reference; // using ScalarStRef_t = typename traits::LScalarMap_t::reference; // create statefields LColl_t coll{}; coll.add_pixel({0}); coll.initialise(); StrainStField_t F_("previous gradient", coll); StrainStField_t be_("previous elastic strain", coll); FlowStField_t eps_("plastic flow", coll); auto F_prev{F_.get_map()}; F_prev[0].current() = Strain_t::Identity(); auto be_prev{be_.get_map()}; be_prev[0].current() = Strain_t::Identity(); auto eps_prev{eps_.get_map()}; eps_prev[0].current() = 0; // elastic deformation Strain_t F{Strain_t::Identity()}; F(0, 1) = 1e-5; F_.cycle(); be_.cycle(); eps_.cycle(); Strain_t stress{Fix::mat.evaluate_stress(F, F_prev[0], be_prev[0], eps_prev[0])}; + if (has_precomputed_values) { + Strain_t tau_ref{}; + tau_ref << + 1.92999522e-11, 3.86000000e-06, 0.00000000e+00, + 3.86000000e-06, -1.93000510e-11, 0.00000000e+00, + 0.00000000e+00, 0.00000000e+00, -2.95741950e-17; + Real error{(tau_ref-stress).norm()}; + BOOST_CHECK_LT(error, hi_tol); + + Strain_t be_ref{}; + be_ref << + 1.00000000e+00, 1.00000000e-05, 0.00000000e+00, + 1.00000000e-05, 1.00000000e+00, 0.00000000e+00, + 0.00000000e+00, 0.00000000e+00, 1.00000000e+00; + error = (be_ref-be_prev[0].current()).norm(); + BOOST_CHECK_LT(error, hi_tol); + + Real ep_ref{0}; + error = ep_ref-eps_prev[0].current(); + BOOST_CHECK_LT(error, hi_tol); + + } if (verbose) { std::cout << "τ =" << std::endl << stress << std::endl; std::cout << "F =" << std::endl << F << std::endl; std::cout << "Fₜ =" << std::endl << F_prev[0].current() << std::endl; std::cout << "bₑ =" << std::endl << be_prev[0].current() << std::endl; std::cout << "εₚ =" << std::endl << eps_prev[0].current() << std::endl; } F_.cycle(); be_.cycle(); eps_.cycle(); - std::cout << "Post Cycle" << std::endl; // plastic deformation F(0, 1) = .2; stress = Fix::mat.evaluate_stress(F, F_prev[0], be_prev[0], eps_prev[0]); + if (has_precomputed_values) { + Strain_t tau_ref{}; + tau_ref << + 1.98151335e-04, 1.98151335e-03, 0.00000000e+00, + 1.98151335e-03, -1.98151335e-04, 0.00000000e+00, + 0.00000000e+00, 0.00000000e+00, 1.60615155e-16; + Real error{(tau_ref-stress).norm()}; + BOOST_CHECK_LT(error, hi_tol); + + Strain_t be_ref{}; + be_ref << + 1.00052666, 0.00513348, 0., + 0.00513348, 0.99949996, 0., + 0., 0., 1.; + error = (be_ref-be_prev[0].current()).norm(); + BOOST_CHECK_LT(error, hi_tol); + + Real ep_ref{0.11229988}; + error = ep_ref-eps_prev[0].current(); + BOOST_CHECK_LT(error, hi_tol); + + } if (verbose) { + std::cout << "Post Cycle" << std::endl; std::cout << "τ =" << std::endl << stress << std::endl << "F =" << std::endl << F << std::endl << "Fₜ =" << std::endl << F_prev[0].current() << std::endl << "bₑ =" << std::endl << be_prev[0].current() << std::endl << "εₚ =" << std::endl << eps_prev[0].current() << std::endl; } } + BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_stiffness, Fix, mats, Fix) { + + // This test uses precomputed reference values (computed using + // elasto-plasticity.py) for the 3d case only + + // need higher tol because of printout precision of reference solutions + constexpr Real hi_tol{1e-8}; + constexpr Dim_t mdim{Fix::mdim}, sdim{Fix::sdim}; + constexpr bool has_precomputed_values{(mdim == sdim) && (mdim == threeD)}; + constexpr bool verbose{has_precomputed_values && false}; + using Strain_t = Eigen::Matrix; + using Stiffness_t = T4Mat; + using traits = MaterialMuSpectre_traits>; + using LColl_t = typename traits::LFieldColl_t; + using StrainStField_t = StateField< + TensorField>; + using FlowStField_t = StateField< + ScalarField>; + + // using StrainStRef_t = typename traits::LStrainMap_t::reference; + // using ScalarStRef_t = typename traits::LScalarMap_t::reference; + + // create statefields + LColl_t coll{}; + coll.add_pixel({0}); + coll.initialise(); + + StrainStField_t F_("previous gradient", coll); + StrainStField_t be_("previous elastic strain", coll); + FlowStField_t eps_("plastic flow", coll); + + auto F_prev{F_.get_map()}; + F_prev[0].current() = Strain_t::Identity(); + auto be_prev{be_.get_map()}; + be_prev[0].current() = Strain_t::Identity(); + auto eps_prev{eps_.get_map()}; + eps_prev[0].current() = 0; + // elastic deformation + Strain_t F{Strain_t::Identity()}; + F(0, 1) = 1e-5; + + F_.cycle(); + be_.cycle(); + eps_.cycle(); + Strain_t stress{}; + Stiffness_t stiffness{}; + std::tie(stress, stiffness) = + Fix::mat.evaluate_stress_tangent(F, + F_prev[0], + be_prev[0], + eps_prev[0]); + + if (has_precomputed_values) { + Strain_t tau_ref{}; + tau_ref << + 1.92999522e-11, 3.86000000e-06, 0.00000000e+00, + 3.86000000e-06, -1.93000510e-11, 0.00000000e+00, + 0.00000000e+00, 0.00000000e+00, -2.95741950e-17; + Real error{(tau_ref-stress).norm()}; + BOOST_CHECK_LT(error, hi_tol); + + Strain_t be_ref{}; + be_ref << + 1.00000000e+00, 1.00000000e-05, 0.00000000e+00, + 1.00000000e-05, 1.00000000e+00, 0.00000000e+00, + 0.00000000e+00, 0.00000000e+00, 1.00000000e+00; + error = (be_ref-be_prev[0].current()).norm(); + BOOST_CHECK_LT(error, hi_tol); + + Real ep_ref{0}; + error = ep_ref-eps_prev[0].current(); + BOOST_CHECK_LT(error, hi_tol); + + Stiffness_t C4_ref{}; + C4_ref << + 0.67383333, 0., 0., 0., 0.28783333, 0., 0., 0., 0.28783333, + 0., 0.193, 0., 0.193, 0., 0., 0., 0., 0., + 0., 0., 0.193, 0., 0., 0., 0.193, 0., 0., + 0., 0.193, 0., 0.193, 0., 0., 0., 0., 0., + 0.28783333, 0., 0., 0., 0.67383333, 0., 0., 0., 0.28783333, + 0., 0., 0., 0., 0., 0.193, 0., 0.193, 0., + 0., 0., 0.193, 0., 0., 0., 0.193, 0., 0., + 0., 0., 0., 0., 0., 0.193, 0., 0.193, 0., + 0.28783333, 0., 0., 0., 0.28783333, 0., 0., 0., 0.67383333; + error = (C4_ref - stiffness).norm(); + BOOST_CHECK_LT(error, hi_tol); + + } + if (verbose) { + std::cout << "C₄ =" << std::endl << stiffness << std::endl; + } + F_.cycle(); + be_.cycle(); + eps_.cycle(); + + // plastic deformation + F(0, 1) = .2; + std::tie(stress, stiffness) = + Fix::mat.evaluate_stress_tangent(F, + F_prev[0], + be_prev[0], + eps_prev[0]); + + if (has_precomputed_values) { + Strain_t tau_ref{}; + tau_ref << + 1.98151335e-04, 1.98151335e-03, 0.00000000e+00, + 1.98151335e-03, -1.98151335e-04, 0.00000000e+00, + 0.00000000e+00, 0.00000000e+00, 1.60615155e-16; + Real error{(tau_ref-stress).norm()}; + BOOST_CHECK_LT(error, hi_tol); + + Strain_t be_ref{}; + be_ref << + 1.00052666, 0.00513348, 0., + 0.00513348, 0.99949996, 0., + 0., 0., 1.; + error = (be_ref-be_prev[0].current()).norm(); + BOOST_CHECK_LT(error, hi_tol); + + Real ep_ref{0.11229988}; + error = ep_ref-eps_prev[0].current(); + BOOST_CHECK_LT(error, hi_tol); + + Stiffness_t C4_ref{}; + C4_ref << + +4.23106224e-01, -4.27959704e-04, 0.00000000e+00, -4.27959704e-04, 4.13218286e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.13175490e-01, + -4.27959704e-04, 7.07167743e-04, 0.00000000e+00, 7.07167743e-04, 4.27959704e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.79121029e-18, + +0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, + -4.27959704e-04, 7.07167743e-04, 0.00000000e+00, 7.07167743e-04, 4.27959704e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.79121029e-18, + +4.13218286e-01, 4.27959704e-04, 0.00000000e+00, 4.27959704e-04, 4.23106224e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.13175490e-01, + +0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, + +0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 0.00000000e+00, + +0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, 4.98676478e-03, 0.00000000e+00, + +4.13175490e-01, 2.79121029e-18, 0.00000000e+00, 2.79121029e-18, 4.13175490e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.23149020e-01; + error = (C4_ref - stiffness).norm(); + BOOST_CHECK_LT(error, hi_tol); + + } + if (verbose) { + std::cout << "Post Cycle" << std::endl; + std::cout << "C₄ =" << std::endl << stiffness << std::endl; + } + } + BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/tests.hh b/tests/tests.hh index 4784f2b..1c83eb3 100644 --- a/tests/tests.hh +++ b/tests/tests.hh @@ -1,61 +1,61 @@ /** * @file tests.hh * * @author Till Junge * * @date 10 May 2017 * * @brief common defs for tests * * Copyright © 2017 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 "common/common.hh" # if defined(__INTEL_COMPILER) //# pragma warning ( disable : 383 ) # elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu # pragma clang diagnostic push # pragma clang diagnostic ignored "-Weffc++" # elif (defined(__GNUC__) || defined(__GNUG__)) # pragma GCC diagnostic push # pragma GCC diagnostic ignored "-Weffc++" # endif #include #include # if defined(__INTEL_COMPILER) //# pragma warning ( disable : 383 ) # elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu # pragma clang diagnostic pop # pragma clang diagnostic ignored "-Weffc++" # elif (defined(__GNUC__) || defined(__GNUG__)) # pragma GCC diagnostic pop # pragma GCC diagnostic ignored "-Weffc++" # endif #ifndef TESTS_H #define TESTS_H namespace muSpectre { - const Real tol = 1e-14*100; //it's in percent + constexpr Real tol = 1e-14*100; //it's in percent } // muSpectre #endif /* TESTS_H */