diff --git a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc index c15bc57d5..fa50b6cac 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc @@ -1,211 +1,212 @@ /** * @file material_plastic.cc * * @author Mohit Pundir * * @date creation: Wed Sep 09 2020 * @date last modification: Wed Sep 09 2020 * * @brief Implementation of the akantu::MaterialDruckerPrager class * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_drucker_prager.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template MaterialDruckerPrager::MaterialDruckerPrager( SolidMechanicsModel & model, const ID & id) : MaterialPlastic(model, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialDruckerPrager:: MaterialDruckerPrager(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : MaterialPlastic(model, dim, mesh, fe_engine, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::initialize() { this->registerParam("phi", phi, Real(0.), _pat_parsable | _pat_modifiable, "Internal friction angle in degrees"); this->registerParam("fc", fc, Real(1.), _pat_parsable | _pat_modifiable, "Compressive strength"); - this->registerParam("radial_return", radial_return_mapping, bool(true), - _pat_parsable | _pat_modifiable, - "Radial return mapping"); + this->registerParam("tolerance", tolerance, Real(1e-8), _pat_parsable | _pat_modifiable, + "projection tolerance"); + this->registerParam("max_iterations", max_iterations, UInt(100), _pat_parsable | _pat_modifiable, + "maximum iterations for projection"); this->registerParam("is_capped", is_capped, bool(false), _pat_parsable | _pat_modifiable, "Capping the yield surface"); this->registerParam("hydrostatic_min", hydrostatic_min, Real(1.), _pat_parsable | _pat_modifiable, "Minimum hydrostatic stress"); this->registerParam("deviatoric_max", deviatoric_max, Real(1.), _pat_parsable | _pat_modifiable, "Maximum deviatoric stress"); this->updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::updateInternalParameters() { MaterialElastic::updateInternalParameters(); // compute alpha and k parameters for Drucker-Prager Real phi_radian = this->phi * M_PI / 180.; this->alpha = (6.*sin(phi_radian))/(3.-sin(phi_radian)); Real cohesion = this->fc * (1. - sin(phi_radian))/(2.*cos(phi_radian)); this->k = (6. * cohesion * cos(phi_radian))/(3. - sin(phi_radian)); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialThermal::computeStress(el_type, ghost_type); // infinitesimal and finite deformation auto sigma_th_it = this->sigma_th(el_type, ghost_type).begin(); auto previous_sigma_th_it = this->sigma_th.previous(el_type, ghost_type).begin(); auto previous_gradu_it = this->gradu.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto previous_stress_it = this->stress.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto inelastic_strain_it = this->inelastic_strain(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto previous_inelastic_strain_it = this->inelastic_strain.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); // // Finite Deformations // if (this->finite_deformation) { auto previous_piola_kirchhoff_2_it = this->piola_kirchhoff_2.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto green_strain_it = this->green_strain(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); auto & inelastic_strain_tensor = *inelastic_strain_it; auto & previous_inelastic_strain_tensor = *previous_inelastic_strain_it; auto & previous_grad_u = *previous_gradu_it; auto & previous_sigma = *previous_piola_kirchhoff_2_it; auto & green_strain = *green_strain_it; this->template gradUToE(grad_u, green_strain); Matrix previous_green_strain(spatial_dimension, spatial_dimension); this->template gradUToE(previous_grad_u, previous_green_strain); Matrix F_tensor(spatial_dimension, spatial_dimension); this->template gradUToF(grad_u, F_tensor); computeStressOnQuad(green_strain, previous_green_strain, sigma, previous_sigma, inelastic_strain_tensor, previous_inelastic_strain_tensor, *sigma_th_it, *previous_sigma_th_it, F_tensor); ++sigma_th_it; ++inelastic_strain_it; ++previous_sigma_th_it; //++previous_stress_it; ++previous_gradu_it; ++green_strain_it; ++previous_inelastic_strain_it; ++previous_piola_kirchhoff_2_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } // Infinitesimal deformations else { MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); auto & inelastic_strain_tensor = *inelastic_strain_it; auto & previous_inelastic_strain_tensor = *previous_inelastic_strain_it; auto & previous_grad_u = *previous_gradu_it; auto & previous_sigma = *previous_stress_it; computeStressOnQuad( grad_u, previous_grad_u, sigma, previous_sigma, inelastic_strain_tensor, previous_inelastic_strain_tensor, *sigma_th_it, *previous_sigma_th_it); ++sigma_th_it; ++inelastic_strain_it; ++previous_sigma_th_it; ++previous_stress_it; ++previous_gradu_it; ++previous_inelastic_strain_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } AKANTU_DEBUG_OUT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::computeTangentModuli( __attribute__((unused)) ElementType el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(plastic_drucker_prager, MaterialDruckerPrager); } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh index 5897011b0..68e98012f 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh @@ -1,167 +1,190 @@ /** * @file material_drucker_pruger.hh * * @author Mohit Pundir * * @date creation: Wed Sep 09 2020 * @date last modification: Wed Sep 09 2020 * * @brief Specialization of the material class for isotropic * plasticity with Drucker-Pruger yield criterion * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_voigthelper.hh" #include "material_plastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DRUCKER_PRAGER_HH__ #define __AKANTU_MATERIAL_DRUCKER_PRAGER_HH__ namespace akantu { /** * Material plastic with a Drucker-Prager yield criterion */ template class MaterialDruckerPrager : public MaterialPlastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDruckerPrager(SolidMechanicsModel & model, const ID & id = ""); MaterialDruckerPrager(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); protected: using voigt_h = VoigtHelper; void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// compute the tangent stiffness matrix for an element type void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override; protected: /// Infinitesimal deformations inline void computeStressOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelas_strain, const Matrix & previous_inelas_strain, const Real & sigma_th, const Real & previous_sigma_th); /// Finite deformations inline void computeStressOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelas_strain, const Matrix & previous_inelas_strain, const Real & sigma_th, const Real & previous_sigma_th, const Matrix & F_tensor); inline void computeTangentModuliOnQuad( Matrix & tangent, const Matrix & grad_u, const Matrix & previous_grad_u, const Matrix & sigma_tensor, const Matrix & previous_sigma_tensor) const; inline Real computeYieldFunction(const Matrix & sigma); inline Real computeYieldStress(const Matrix & sigma); inline void computeDeviatoricStress(const Matrix & sigma, Matrix & sigma_dev); inline bool aboveThresholdStress(const Matrix &, Real & , Real & ); - inline void projectStressOnThreshold(Matrix &, Vector &, - Real &, Real &); + inline void projectionOnThreshold(Matrix &, Vector &, + Real &, Real &); - /// rcompute the alpha and k parameters + /// compute the alpha and k parameters void updateInternalParameters() override; + + /* ------------------------------------------------------------------------ */ + /* Methods for projection to the yield surface */ + /* ------------------------------------------------------------------------ */ + /// compute the objective function to minimize + virtual inline Real computeObjectiveFunction(const Matrix & sigma_guess, + const Matrix & sigma_trial, + Real & plastic_multiplier_guess, + Vector & gradient_f, + Vector & delta_inelastic_strain, + Real & yield_function); + + /// compute the jacobian of the objective function + virtual inline void computeJacobian(const Matrix &, Vector &); + + /// compute the hessian of the objective function + virtual inline void computeHessian(const Matrix &, Matrix &); + + + public: // closet point projection method to compute stress state on the // yield surface inline void computeGradientAndPlasticMultplier( const Matrix & sigma_tr, Real & plastic_multiplier_guess, Vector & gradient_f, Vector & delta_inelastic_strain, UInt max_iterations = 100, Real tolerance = 1e-8); public: /// get modified compressive strength AKANTU_GET_MACRO(K, k, Real); /// get modified friction angle AKANTU_GET_MACRO(Alpha, alpha, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: // Internal friction angle of the material Real phi; // Compressive strength of the material Real fc; // modified friction angle for Drucker-Prager Real alpha; // modified compressive strength for Drucker-Prager Real k; + // maximum number of iterations for projection + UInt max_iterations; + + // tolerance error for projection + Real tolerance; + // maximum deviatoric for capped Drucker-Prager Real deviatoric_max; // minimum hydrostatic for capped Drucker-Prager Real hydrostatic_min; - // radial return mapping - bool radial_return_mapping; - // capped bool is_capped; }; } // namespace akantu #include "material_drucker_prager_inline_impl.hh" #endif /*__AKANTU_MATERIAL_DRUCKER_PRAGER_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh index 194bb4c96..851952ff7 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh @@ -1,719 +1,727 @@ /** * @file material_linear_isotropic_hardening_inline_impl.hh * * @author Mohit Pundir * * @date creation: Wed Sep 09 2020 * @date last modification: Wed Sep 09 2020 * * @brief Implementation of the inline functions of the material * Drucker-Prager * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #include "material_drucker_prager.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /// Deviatoric Stress template inline void MaterialDruckerPrager::computeDeviatoricStress(const Matrix & sigma, Matrix & sigma_dev){ for (UInt i = 0; i < dim; ++i){ for (UInt j = 0; j < dim; ++j){ sigma_dev(i, j) = sigma(i, j); } } sigma_dev -= Matrix::eye(dim, sigma.trace() / dim); } /* -------------------------------------------------------------------------- */ /// Yield Stress template inline Real MaterialDruckerPrager::computeYieldStress(const Matrix & sigma) { return this->alpha * sigma.trace() - this->k; } /* -------------------------------------------------------------------------- */ /// Yield function template inline Real MaterialDruckerPrager::computeYieldFunction(const Matrix & sigma) { Matrix sigma_dev(dim, dim, 0); this->computeDeviatoricStress(sigma, sigma_dev); // compute deviatoric invariant J2 Real j2 = (1./ 2.) * sigma_dev.doubleDot(sigma_dev); Real sigma_dev_eff = std::sqrt(3. * j2); Real modified_yield_stress = computeYieldStress(sigma); return sigma_dev_eff + modified_yield_stress; } /* -------------------------------------------------------------------------- */ /// Threshold Stress template inline bool MaterialDruckerPrager::aboveThresholdStress(const Matrix & sigma, Real & k, Real & alpha) { Matrix sigma_dev(dim, dim, 0); this->computeDeviatoricStress(sigma, sigma_dev); Real I1 = sigma.trace(); Real J2 = (1./2) *sigma_dev.doubleDot(sigma_dev); if(I1 > k/alpha) { Real beta = sqrt(3.*J2)/(I1-k/alpha); Real gamma = tan(M_PI/2. - atan(alpha)); if(beta < gamma) return true; } return false; } + +/* -------------------------------------------------------------------------- */ +template +inline Real MaterialDruckerPrager::computeObjectiveFunction(const Matrix & sigma_guess, + const Matrix & sigma_trial, + Real & plastic_multiplier_guess, + Vector & gradient_f, + Vector & delta_inelastic_strain, + Real & yield_function){ + + UInt size = voigt_h::size; + + // elastic stifnness tensor + Matrix De(size, size, 0.); + MaterialElastic::computeTangentModuliOnQuad(De); + + // elastic compliance tensor + Matrix Ce(size, size, 0.); + Ce.inverse(De); + + // objective function to be computed + Vector f(size, 0.); + + // compute gradient + computeJacobian(sigma_guess, gradient_f); + + // compute yield function + yield_function = computeYieldFunction(sigma_guess); + + // compute increment strain + auto sigma_trial_voigt = voigt_h::matrixToVoigt(sigma_trial); + auto sigma_guess_voigt = voigt_h::matrixToVoigt(sigma_guess); + auto tmp = sigma_trial_voigt - sigma_guess_voigt; + delta_inelastic_strain.mul(Ce, tmp); + + // compute objective function + f.mul(De, gradient_f, plastic_multiplier_guess); + f = tmp - f; + + // compute error + Real error = std::max(f.norm(), std::abs(yield_function)); + return error; +} + +/* -------------------------------------------------------------------------- */ +template +inline void MaterialDruckerPrager::computeJacobian(const Matrix & sigma_guess, + Vector & gradient_f) { + + UInt size = voigt_h::size; + + // scaling matrix for computing gradient and hessian from voigt notation + Matrix scaling_matrix(size, size, 0.); + scaling_matrix.eye(1.); + for(auto i : arange(dim, size)) { + for(auto j : arange(dim, size)){ + scaling_matrix(i, j) *= 2.; + } + } + + // krocker delta vector in voigt notation + Vector kronecker_delta(size, 0.); + for(auto i : arange(dim)) { + kronecker_delta[i] = 1.; + } + + // compute deviatoric invariant + Matrix sigma_dev(dim, dim, 0); + + this->computeDeviatoricStress(sigma_guess, sigma_dev); + + Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); + + Vector sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); + + gradient_f.mul(scaling_matrix, sigma_dev_voigt, 3./ (2. * std::sqrt(3. * j2)) ); + gradient_f += alpha * kronecker_delta; +} + +/* -------------------------------------------------------------------------- */ +template +inline void MaterialDruckerPrager::computeHessian(const Matrix & sigma_guess, + Matrix & hessian_f) { + + UInt size = voigt_h::size; + + // scaling matrix for computing gradient and hessian from voigt notation + Matrix scaling_matrix(size, size, 0.); + scaling_matrix.eye(1.); + for(auto i : arange(dim, size)) { + for(auto j : arange(dim, size)){ + scaling_matrix(i, j) *= 2.; + } + } + + // krocker delta vector in voigt notation + Vector kronecker_delta(size, 0.); + for(auto i : arange(dim)) { + kronecker_delta[i] = 1.; + } + + // compute deviatoric invariant J2 + Matrix sigma_dev(dim, dim, 0); + + this->computeDeviatoricStress(sigma_guess, sigma_dev); + + Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); + auto sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); + + Vector temp(sigma_dev_voigt.size()); + temp.mul(scaling_matrix, sigma_dev_voigt); + + Matrix id(kronecker_delta.size(), kronecker_delta.size()); + id.outerProduct(kronecker_delta, kronecker_delta); + id *= -1./3.; + id += Matrix::eye(kronecker_delta.size(), 1.); + + Matrix tmp3(kronecker_delta.size(), kronecker_delta.size()); + tmp3.mul(scaling_matrix, id); + hessian_f.outerProduct(temp, temp); + hessian_f *= -9./(4.* pow(3.*j2, 3./2.)); + hessian_f += (3./(2.* pow(3.*j2, 1./2.)))*tmp3; +} + /* -------------------------------------------------------------------------- */ template inline void MaterialDruckerPrager::computeGradientAndPlasticMultplier( const Matrix & sigma_trial, Real & plastic_multiplier_guess, Vector & gradient_f, Vector & delta_inelastic_strain, UInt max_iterations, Real tolerance) { UInt size = voigt_h::size; - + // guess stress state at each iteration, initial guess is the trial state Matrix sigma_guess(sigma_trial); + // hessian matrix of yield surface + Matrix hessian_f(size, size, 0.); + + // variation in sigma at each iteration + Vector delta_sigma(size, 0.); + // plastic multiplier guess at each iteration, initial guess is zero plastic_multiplier_guess = 0.; // gradient of yield surface in voigt notation gradient_f.zero(); // plastic strain increment at each iteration delta_inelastic_strain.zero(); - - // variation in sigma at each iteration - Vector delta_sigma(size, 0.); - - // krocker delta vector in voigt notation - Vector kronecker_delta(size, 0.); - for(auto i : arange(dim)) { - kronecker_delta[i] = 1.; - } - - // hessian matrix of yield surface - Matrix hessian_f(size, size, 0.); - // scaling matrix for computing gradient and hessian from voigt notation - Matrix scaling_matrix(size, size, 0.); - scaling_matrix.eye(1.); - for(auto i : arange(dim, size)) { - for(auto j : arange(dim, size)){ - scaling_matrix(i, j) *= 2.; - } - } - // elastic stifnness tensor Matrix De(size, size, 0.); MaterialElastic::computeTangentModuliOnQuad(De); // elastic compliance tensor Matrix Ce(size, size, 0.); Ce.inverse(De); - - // objective function to be computed - Vector f(size, 0.); - + // yield function value at each iteration Real yield_function; /* ------------------------------------------------------- */ - /* projecting stress state to origin of yield function */ + /* Projecting stress state to origin of yield function */ /* if first invariant is greater than the threshold */ /* ------------------------------------------------------- */ if(this->aboveThresholdStress(sigma_guess, k, alpha) and this->alpha > 0) { - this->projectStressOnThreshold(sigma_guess, gradient_f, - plastic_multiplier_guess, tolerance); + this->projectionOnThreshold(sigma_guess, gradient_f, + plastic_multiplier_guess, tolerance); auto delta_sigma_final = sigma_trial - sigma_guess; auto delta_sigma_voigt = voigt_h::matrixToVoigt(delta_sigma_final); delta_inelastic_strain.mul(Ce, delta_sigma_voigt); - return; } /* -------------------------------------------------- */ /* If Drucker-prager is Capped */ /* -------------------------------------------------- */ /*Real I1 = sigma_guess.trace(); if(is_capped) { Matrix sigma_dev(dim, dim, 0); this->computeDeviatoricStress(sigma_guess, sigma_dev); Real J2 = (1./2) *sigma_dev.doubleDot(sigma_dev); if(I1 < hydrostatic_min and sqrt(3.*J2) < deviatoric_max) { auto projection_error = I1 - hydrostatic_min; while(tolerance < std::abs(projection_error)) { Matrix delta_sigma(dim, dim); Matrix jacobian(dim, dim); Matrix jacobian_inv(dim, dim); jacobian = Matrix::eye(dim, dim); jacobian_inv.inverse(jacobian); delta_sigma += -projection_error * jacobian_inv; sigma_guess += delta_sigma; projection_error = sigma_guess.trace() - hydrostatic_min; } auto delta_sigma_final = sigma_trial - sigma_guess; auto delta_sigma_voigt = voigt_h::matrixToVoigt(delta_sigma_final); delta_inelastic_strain.mul(Ce, delta_sigma_voigt); return; } Real beta = (sqrt(3.*J2) - deviatoric_max )/(I1- hydrostatic_min); Real gamma = tan(M_PI/2. - atan(alpha)); if((I1 > hydrostatic_min and beta > gamma) or (I1 < hydrostatic_min and sqrt(3.*J2) > deviatoric_max )) { // lambda function to compute gradient of yield surface in voigt notation auto compute_gradient_f3 = [&sigma_guess, &scaling_matrix, &kronecker_delta, &gradient_f](){ Matrix sigma_dev(dim, dim, 0); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma_guess(i, j); sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); Vector sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); // compute deviatoric invariant Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); gradient_f.mul(scaling_matrix, sigma_dev_voigt, 3./ (2. * std::sqrt(3. * j2)) ); }; auto projection_error = sqrt(3.*J2) - deviatoric_max; while( tolerance < projection_error ) { compute_gradient_f3(); Vector tmp1(size); tmp1.mul(Ce, gradient_f); auto denominator = gradient_f.dot(tmp1); auto delta_plastic_multiplier = projection_error; delta_plastic_multiplier /= denominator; plastic_multiplier_guess += delta_plastic_multiplier; delta_sigma.mul(Ce, gradient_f, delta_plastic_multiplier); // compute sigma_guess from voigt notation Matrix delta_sigma_mat(dim, dim); voigt_h::voigtToMatrix(delta_sigma, delta_sigma_mat); sigma_guess -= delta_sigma_mat; this->computeDeviatoricStress(sigma_guess, sigma_dev); J2 = (1./2) *sigma_dev.doubleDot(sigma_dev); projection_error = sqrt(3.*J2) - deviatoric_max; } projection_error = sigma_guess.trace() - hydrostatic_min; while(tolerance < std::abs(projection_error)) { Matrix delta_sigma(dim, dim); Matrix jacobian(dim, dim); Matrix jacobian_inv(dim, dim); jacobian = Matrix::eye(dim, dim); jacobian_inv.inverse(jacobian); delta_sigma += -projection_error * jacobian_inv; sigma_guess += delta_sigma; projection_error = sigma_guess.trace() - hydrostatic_min; } auto delta_sigma_final = sigma_trial - sigma_guess; auto delta_sigma_voigt = voigt_h::matrixToVoigt(delta_sigma_final); delta_inelastic_strain.mul(Ce, delta_sigma_voigt); return; } }*/ - - - // lambda function to compute gradient of yield surface in voigt notation - auto compute_gradient_f = [&sigma_guess, &scaling_matrix, &kronecker_delta, - &gradient_f](Real & alpha){ - - Matrix sigma_dev(dim, dim, 0); - - for (UInt i = 0; i < dim; ++i) - for (UInt j = 0; j < dim; ++j) - sigma_dev(i, j) = sigma_guess(i, j); - - sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); - - - Vector sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); - - // compute deviatoric invariant - Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); - - gradient_f.mul(scaling_matrix, sigma_dev_voigt, 3./ (2. * std::sqrt(3. * j2)) ); - gradient_f += alpha * kronecker_delta; - }; /* -------------------------------------------------- */ /* Generalized cutting plane algorithm (explicit) */ /* -------------------------------------------------- */ /*UInt iteration = 0; auto error = this->computeYieldFunction(sigma_guess); while( tolerance < error ) { compute_gradient_f(alpha); Vector tmp1(size); tmp1.mul(Ce, gradient_f); auto denominator = gradient_f.dot(tmp1); auto delta_plastic_multiplier = error; delta_plastic_multiplier /= denominator; plastic_multiplier_guess += delta_plastic_multiplier; delta_sigma.mul(Ce, gradient_f, delta_plastic_multiplier); // compute sigma_guess from voigt notation Matrix delta_sigma_mat(dim, dim); voigt_h::voigtToMatrix(delta_sigma, delta_sigma_mat); sigma_guess -= delta_sigma_mat; // compute increment inelastic strain auto sigma_trial_voigt = voigt_h::matrixToVoigt(sigma_trial); auto sigma_guess_voigt = voigt_h::matrixToVoigt(sigma_guess); auto tmp = sigma_trial_voigt - sigma_guess_voigt; delta_inelastic_strain.mul(Ce, tmp); error = this->computeYieldFunction(sigma_guess); iteration++; }*/ - /* ------------------------------------- */ - /* Closet point projection (implicit) */ - /* ------------------------------------- */ - - - // lambda function to compute hessian matrix of yield surface - auto compute_hessian_f = [&sigma_guess, &hessian_f, &scaling_matrix, - &kronecker_delta](){ - - Matrix sigma_dev(dim, dim, 0); - - for (UInt i = 0; i < dim; ++i) - for (UInt j = 0; j < dim; ++j) - sigma_dev(i, j) = sigma_guess(i, j); - - sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); - - auto sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); - - // compute deviatoric invariant J2 - Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); - - Vector temp(sigma_dev_voigt.size()); - temp.mul(scaling_matrix, sigma_dev_voigt); - - Matrix id(kronecker_delta.size(), kronecker_delta.size()); - id.outerProduct(kronecker_delta, kronecker_delta); - id *= -1./3.; - id += Matrix::eye(kronecker_delta.size(), 1.); - - Matrix tmp3(kronecker_delta.size(), kronecker_delta.size()); - tmp3.mul(scaling_matrix, id); - hessian_f.outerProduct(temp, temp); - hessian_f *= -9./(4.* pow(3.*j2, 3./2.)); - hessian_f += (3./(2.* pow(3.*j2, 1./2.)))*tmp3; - }; - - /* --------------------------- */ - /* init before iteration loop */ - /* --------------------------- */ - auto update_f = [&f, &sigma_guess, &sigma_trial, &plastic_multiplier_guess, &Ce, &De, - &yield_function, &gradient_f, &delta_inelastic_strain, - &compute_gradient_f](Real & k, Real & alpha){ - - // compute gradient - compute_gradient_f(alpha); - - // compute yield function - Matrix sigma_dev(dim, dim, 0); - - for (UInt i = 0; i < dim; ++i) - for (UInt j = 0; j < dim; ++j) - sigma_dev(i, j) = sigma_guess(i, j); - - sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); - - Real j2 = (1./ 2.) * sigma_dev.doubleDot(sigma_dev); - Real sigma_dev_eff = std::sqrt(3. * j2); - - Real modified_yield_stress = alpha * sigma_guess.trace() - k; - - yield_function = sigma_dev_eff + modified_yield_stress; - - // compute increment strain - auto sigma_trial_voigt = voigt_h::matrixToVoigt(sigma_trial); - auto sigma_guess_voigt = voigt_h::matrixToVoigt(sigma_guess); - auto tmp = sigma_trial_voigt - sigma_guess_voigt; - delta_inelastic_strain.mul(Ce, tmp); + /* ----------------------------------------------------------------------- */ + /* Projection on yield surface using Closet point projection (implicit) */ + /* ----------------------------------------------------------------------- */ - // compute objective function - f.mul(De, gradient_f, plastic_multiplier_guess); - f = tmp - f; - - // compute error - auto error = std::max(f.norm(), std::abs(yield_function)); - return error; - }; - - Real alpha_tmp{alpha}; - Real k_tmp{k}; - if(radial_return_mapping){ - alpha_tmp = 0; - k_tmp = std::abs(alpha*sigma_guess.trace() - k); - } + Real projection_error = computeObjectiveFunction(sigma_guess, sigma_trial, plastic_multiplier_guess, + gradient_f, delta_inelastic_strain, yield_function); - auto projection_error = update_f(k_tmp , alpha_tmp); - /* --------------------------- */ /* iteration loop */ /* --------------------------- */ Matrix xi(size, size); Matrix xi_inv(size, size); Vector tmp(size); Vector tmp1(size); Matrix tmp2(size, size); UInt iterations{0}; while(tolerance < projection_error and iterations < max_iterations) { // compute hessian at previous step - compute_hessian_f(); - + computeHessian(sigma_guess, hessian_f); + // compute inverse matrix Xi xi = Ce + plastic_multiplier_guess * hessian_f; // compute inverse matrix Xi xi_inv.inverse(xi); tmp.mul(xi_inv, gradient_f); auto denominator = gradient_f.dot(tmp); // compute plastic multplier guess - tmp1.mul(xi_inv, delta_inelastic_strain); plastic_multiplier_guess = gradient_f.dot(tmp1); plastic_multiplier_guess += yield_function; plastic_multiplier_guess /= denominator; // compute delta sigma - tmp2.outerProduct(tmp, tmp); tmp2 /= denominator; tmp2 = xi_inv - tmp2; delta_sigma.mul(tmp2, delta_inelastic_strain); delta_sigma -= tmp*yield_function/denominator; // compute sigma_guess Matrix delta_sigma_mat(dim, dim); voigt_h::voigtToMatrix(delta_sigma, delta_sigma_mat); sigma_guess += delta_sigma_mat; - projection_error = update_f(k_tmp, alpha_tmp); + projection_error = computeObjectiveFunction(sigma_guess, sigma_trial, plastic_multiplier_guess, + gradient_f, delta_inelastic_strain, yield_function); + iterations++; } + if ( iterations >= max_iterations and tolerance < projection_error) { + AKANTU_ERROR("Maximum iterations reached in projection without reaching the tolerance " + << " Increase the iterations or reduce the tolerance " + << " Error is " << projection_error); + } + } /* -------------------------------------------------------------------------- */ /// Infinitesimal deformations template inline void MaterialDruckerPrager::computeStressOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelastic_strain, const Matrix & previous_inelastic_strain, const Real & sigma_th, const Real & previous_sigma_th) { Real delta_sigma_th = sigma_th - previous_sigma_th; Matrix grad_delta_u(grad_u); grad_delta_u -= previous_grad_u; // Compute trial stress, sigma_tr Matrix sigma_tr(dim, dim); MaterialElastic::computeStressOnQuad(grad_delta_u, sigma_tr, delta_sigma_th); sigma_tr += previous_sigma; // Compute the yielding sress bool initial_yielding = (this->computeYieldFunction(sigma_tr) > 0); bool above_capped = is_capped and (sigma_tr.trace() < hydrostatic_min); // use closet point projection to compute the plastic multiplier and // gradient and inealstic strain at the surface for the given trial // stress state Matrix delta_inelastic_strain(dim, dim, 0.); if(initial_yielding or above_capped) { UInt size = voigt_h::size; // plastic multiplier Real dp{0.}; // gradient of yield surface in voigt notation Vector gradient_f(size, 0.); // inelastic strain in voigt notation Vector delta_inelastic_strain_voigt(size, 0.); // compute using closet-point projection this->computeGradientAndPlasticMultplier(sigma_tr, dp, gradient_f, delta_inelastic_strain_voigt); for(auto i: arange(dim, size)) delta_inelastic_strain_voigt[i] /= 2.; voigt_h::voigtToMatrix(delta_inelastic_strain_voigt, delta_inelastic_strain); } // Compute the increment in inelastic strain MaterialPlastic::computeStressAndInelasticStrainOnQuad( grad_delta_u, sigma, previous_sigma, inelastic_strain, previous_inelastic_strain, delta_inelastic_strain); } /* -------------------------------------------------------------------------- */ template -inline void MaterialDruckerPrager::projectStressOnThreshold(Matrix & sigma_guess, +inline void MaterialDruckerPrager::projectionOnThreshold(Matrix & sigma_guess, Vector & gradient_f, Real & plastic_multiplier_guess, Real & tolerance){ UInt size = voigt_h::size; // elastic stifnness tensor Matrix De(size, size, 0.); MaterialElastic::computeTangentModuliOnQuad(De); // elastic compliance tensor Matrix Ce(size, size, 0.); Ce.inverse(De); // variation in sigma at each iteration Vector delta_sigma(size, 0.); // krocker delta vector in voigt notation Vector kronecker_delta(size, 0.); for(auto i : arange(dim)) { kronecker_delta[i] = 1.; } // hessian matrix of yield surface Matrix hessian_f(size, size, 0.); // scaling matrix for computing gradient and hessian from voigt notation Matrix scaling_matrix(size, size, 0.); scaling_matrix.eye(1.); for(auto i : arange(dim, size)) { for(auto j : arange(dim, size)){ scaling_matrix(i, j) *= 2.; } } auto update_first_obj = [&sigma_guess]() { Matrix sigma_dev(dim, dim, 0); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma_guess(i, j); sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); Real J2 = (1./2) *sigma_dev.doubleDot(sigma_dev); auto error = sqrt(3.*J2); return error; }; auto update_sec_obj = [&sigma_guess](Real & k, Real & alpha) { auto error = alpha*sigma_guess.trace() - k; return error; }; // lambda function to compute gradient of yield surface in voigt notation auto compute_gradient_f2 = [&sigma_guess, &scaling_matrix, &kronecker_delta, &gradient_f](){ Matrix sigma_dev(dim, dim, 0); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma_guess(i, j); sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); Vector sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); // compute deviatoric invariant Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); gradient_f.mul(scaling_matrix, sigma_dev_voigt, 3./ (2. * std::sqrt(3. * j2)) ); }; auto projection_error = update_first_obj(); while( tolerance < projection_error ) { compute_gradient_f2(); Vector tmp1(size); tmp1.mul(Ce, gradient_f); auto denominator = gradient_f.dot(tmp1); auto delta_plastic_multiplier = projection_error; delta_plastic_multiplier /= denominator; plastic_multiplier_guess += delta_plastic_multiplier; delta_sigma.mul(Ce, gradient_f, delta_plastic_multiplier); // compute sigma_guess from voigt notation Matrix delta_sigma_mat(dim, dim); voigt_h::voigtToMatrix(delta_sigma, delta_sigma_mat); sigma_guess -= delta_sigma_mat; projection_error = update_first_obj(); } projection_error = update_sec_obj(k, alpha); while(tolerance < projection_error) { Matrix delta_sigma(dim, dim); Matrix jacobian(dim, dim); Matrix jacobian_inv(dim, dim); jacobian = this->alpha * Matrix::eye(dim, dim); jacobian_inv.inverse(jacobian); delta_sigma += -projection_error * jacobian_inv; sigma_guess += delta_sigma; projection_error = update_sec_obj(k, alpha); } } /* -------------------------------------------------------------------------- */ /// Finite deformations template inline void MaterialDruckerPrager::computeStressOnQuad( __attribute__((unused)) const Matrix & grad_u, __attribute__((unused)) const Matrix & previous_grad_u, __attribute__((unused)) Matrix & sigma, __attribute__((unused)) const Matrix & previous_sigma, __attribute__((unused)) Matrix & inelastic_strain, __attribute__((unused)) const Matrix & previous_inelastic_strain, __attribute__((unused)) const Real & sigma_th, __attribute__((unused)) const Real & previous_sigma_th, __attribute__((unused)) const Matrix & F_tensor) { } /* -------------------------------------------------------------------------- */ template inline void MaterialDruckerPrager::computeTangentModuliOnQuad( __attribute__((unused)) Matrix & tangent, __attribute__((unused)) const Matrix & grad_u, __attribute__((unused)) const Matrix & previous_grad_u, __attribute__((unused)) const Matrix & sigma_tensor, __attribute__((unused)) const Matrix & previous_sigma_tensor) const { } } diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc index 7b3c460b7..36e40bbe0 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc @@ -1,264 +1,316 @@ /** * @file test_plastic_materials.cc * * @author Guillaume Anciaux * * @date creation: Fri Nov 17 2017 * @date last modification: Wed Feb 21 2018 * * @brief Tests the plastic material * * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_linear_isotropic_hardening.hh" #include "material_drucker_prager.hh" #include "solid_mechanics_model.hh" #include "test_material_fixtures.hh" #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; using mat_types = ::testing::Types< // Traits, // Traits, Traits, Traits>; /* -------------------------------------------------------------------------- */ template <> void FriendMaterial>::testComputeStress() { Real E = 1.; // Real nu = .3; Real nu = 0.; Real rho = 1.; Real sigma_0 = 1.; Real h = 0.; Real bulk_modulus_K = E / 3. / (1 - 2. * nu); Real shear_modulus_mu = 0.5 * E / (1 + nu); setParam("E", E); setParam("nu", nu); setParam("rho", rho); setParam("sigma_y", sigma_0); setParam("h", h); auto rotation_matrix = getRandomRotation(); Real max_strain = 10.; Real strain_steps = 100; Real dt = max_strain / strain_steps; std::vector steps(strain_steps); std::iota(steps.begin(), steps.end(), 0.); Matrix previous_grad_u_rot(3, 3, 0.); Matrix previous_sigma(3, 3, 0.); Matrix previous_sigma_rot(3, 3, 0.); Matrix inelastic_strain_rot(3, 3, 0.); Matrix inelastic_strain(3, 3, 0.); Matrix previous_inelastic_strain(3, 3, 0.); Matrix previous_inelastic_strain_rot(3, 3, 0.); Matrix sigma_rot(3, 3, 0.); Real iso_hardening = 0.; Real previous_iso_hardening = 0.; // hydrostatic loading (should not plastify) for (auto && i : steps) { auto t = i * dt; auto grad_u = this->getHydrostaticStrain(t); auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix); this->computeStressOnQuad(grad_u_rot, previous_grad_u_rot, sigma_rot, previous_sigma_rot, inelastic_strain_rot, previous_inelastic_strain_rot, iso_hardening, previous_iso_hardening, 0., 0.); auto sigma = this->reverseRotation(sigma_rot, rotation_matrix); Matrix sigma_expected = t * 3. * bulk_modulus_K * Matrix::eye(3, 1.); Real stress_error = (sigma - sigma_expected).norm(); ASSERT_NEAR(stress_error, 0., 1e-13); ASSERT_NEAR(inelastic_strain_rot.norm(), 0., 1e-13); previous_grad_u_rot = grad_u_rot; previous_sigma_rot = sigma_rot; previous_inelastic_strain_rot = inelastic_strain_rot; previous_iso_hardening = iso_hardening; } // deviatoric loading (should plastify) // stress at onset of plastication Real beta = sqrt(42); Real t_P = sigma_0 / 2. / shear_modulus_mu / beta; Matrix sigma_P = sigma_0 / beta * this->getDeviatoricStrain(1.); for (auto && i : steps) { auto t = i * dt; auto grad_u = this->getDeviatoricStrain(t); auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix); Real iso_hardening{0.}; Real previous_iso_hardening{0.}; this->computeStressOnQuad(grad_u_rot, previous_grad_u_rot, sigma_rot, previous_sigma_rot, inelastic_strain_rot, previous_inelastic_strain_rot, iso_hardening, previous_iso_hardening, 0., 0.); auto sigma = this->reverseRotation(sigma_rot, rotation_matrix); auto inelastic_strain = this->reverseRotation(inelastic_strain_rot, rotation_matrix); if (t < t_P) { Matrix sigma_expected = shear_modulus_mu * (grad_u + grad_u.transpose()); Real stress_error = (sigma - sigma_expected).norm(); ASSERT_NEAR(stress_error, 0., 1e-13); ASSERT_NEAR(inelastic_strain_rot.norm(), 0., 1e-13); } else if (t > t_P + dt) { // skip the transition from non plastic to plastic auto delta_lambda_expected = dt / t * previous_sigma.doubleDot(grad_u + grad_u.transpose()) / 2.; auto delta_inelastic_strain_expected = delta_lambda_expected * 3. / 2. / sigma_0 * previous_sigma; auto inelastic_strain_expected = delta_inelastic_strain_expected + previous_inelastic_strain; ASSERT_NEAR((inelastic_strain - inelastic_strain_expected).norm(), 0., 1e-13); auto delta_sigma_expected = 2. * shear_modulus_mu * (0.5 * dt / t * (grad_u + grad_u.transpose()) - delta_inelastic_strain_expected); auto delta_sigma = sigma - previous_sigma; ASSERT_NEAR((delta_sigma_expected - delta_sigma).norm(), 0., 1e-13); } previous_sigma = sigma; previous_sigma_rot = sigma_rot; previous_grad_u_rot = grad_u_rot; previous_inelastic_strain = inelastic_strain; previous_inelastic_strain_rot = inelastic_strain_rot; } } template <> void FriendMaterial>::testComputeStress() { Real E = 1.; Real nu = 0; Real rho = 1.; Real fc = 1.; Real angle = 30.; setParam("E", E); setParam("nu", nu); setParam("rho", rho); setParam("phi", angle); setParam("fc", fc); - setParam("radial_return", false); updateInternalParameters(); Matrix previous_grad_u(3, 3, 0.); Matrix previous_sigma(3, 3, 0.); Matrix inelastic_strain(3, 3, 0.); Matrix previous_inelastic_strain(3, 3, 0.); Matrix sigma(3, 3, 0.); Real k = this->getK(); Real alpha = this->getAlpha(); Real max_strain = 10.*k/alpha; Real strain_steps = 100; Real dt = max_strain / strain_steps; std::vector steps(strain_steps); std::iota(steps.begin(), steps.end(), 0.); // hydrostatic loading (should plastify only when I1 > k/alpha) for (auto && i : steps) { auto t = i * dt; auto grad_u = this->getHydrostaticStrain(t); this->computeStressOnQuad(grad_u, previous_grad_u, sigma, previous_sigma, inelastic_strain, previous_inelastic_strain, 0., 0.); Matrix sigma_elastic(3, 3, 0); MaterialElastic<3>::computeStressOnQuad(grad_u, sigma_elastic, 0); Real I1_elastic = sigma_elastic.trace(); Real I1 = sigma.trace(); if (I1_elastic < k/alpha) { Real stress_error = (sigma - sigma_elastic).norm(); ASSERT_NEAR(stress_error, 0., 1e-13); ASSERT_NEAR(inelastic_strain.norm(), 0., 1e-13); } else if (I1_elastic > k/alpha) { ASSERT_NEAR(I1 , k/alpha, 1e-13); } previous_grad_u = grad_u; previous_sigma = sigma; previous_inelastic_strain = inelastic_strain; } + + previous_grad_u.zero(); + previous_sigma.zero(); + inelastic_strain.zero(); + previous_inelastic_strain.zero(); + sigma.zero(); + // deviatoric loading (should plastify only when J2 > k) + for (auto && i : steps) { + auto t = i * dt; + + auto grad_u = this->getDeviatoricStrain(t); + + this->computeStressOnQuad(grad_u, previous_grad_u, sigma, + previous_sigma, inelastic_strain, + previous_inelastic_strain, 0., 0.); + + Matrix sigma_elastic(3, 3, 0); + MaterialElastic<3>::computeStressOnQuad(grad_u, sigma_elastic, 0); + + Matrix sigma_elastic_dev(3, 3, 0); + this->computeDeviatoricStress(sigma_elastic, sigma_elastic_dev); + + // compute deviatoric invariant J2 + Real j2_elastic = (1./ 2.) * sigma_elastic_dev.doubleDot(sigma_elastic_dev); + Real sigma_elastic_dev_eff = std::sqrt(3. * j2_elastic); + + + Matrix sigma_dev(3, 3, 0); + this->computeDeviatoricStress(sigma, sigma_dev); + + // compute deviatoric invariant J2 + Real j2 = (1./ 2.) * sigma_dev.doubleDot(sigma_dev); + Real sigma_dev_eff = std::sqrt(3. * j2); + + Real I1 = sigma.trace(); + + if (sigma_elastic_dev_eff < k) { + Real stress_error = (sigma - sigma_elastic).norm(); + + ASSERT_NEAR(stress_error, 0., 1e-13); + ASSERT_NEAR(inelastic_strain.norm(), 0., 1e-13); + } + else if (sigma_elastic_dev_eff > k) { + + auto error = alpha*I1 + sigma_dev_eff - k; + ASSERT_NEAR(error, 0., 1e-13); + } + + previous_grad_u = grad_u; + previous_sigma = sigma; + previous_inelastic_strain = inelastic_strain; + } + } namespace { template class TestPlasticMaterialFixture : public ::TestMaterialFixture {}; TYPED_TEST_SUITE(TestPlasticMaterialFixture, mat_types, ); TYPED_TEST(TestPlasticMaterialFixture, ComputeStress) { this->material->testComputeStress(); } TYPED_TEST(TestPlasticMaterialFixture, DISABLED_EnergyDensity) { this->material->testEnergyDensity(); } TYPED_TEST(TestPlasticMaterialFixture, DISABLED_ComputeTangentModuli) { this->material->testComputeTangentModuli(); } TYPED_TEST(TestPlasticMaterialFixture, DISABLED_ComputeCelerity) { this->material->testCelerity(); } } // namespace /*****************************************************************/