diff --git a/packages/solid_mechanics.cmake b/packages/solid_mechanics.cmake index 0ea58434b..901d29ad7 100644 --- a/packages/solid_mechanics.cmake +++ b/packages/solid_mechanics.cmake @@ -1,134 +1,137 @@ #=============================================================================== # @file solid_mechanics.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for core # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 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 . # #=============================================================================== package_declare(solid_mechanics DEFAULT ON DESCRIPTION "Solid mechanics model" DEPENDS core lapack ) package_declare_sources(solid_mechanics model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.hh model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.hh model/solid_mechanics/solid_mechanics_model_io.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.hh model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.cc model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.hh model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.hh + model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc + model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh + model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.hh model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh model/solid_mechanics/materials/material_non_local.hh model/solid_mechanics/materials/material_non_local_tmpl.hh model/solid_mechanics/materials/material_non_local_includes.hh ) package_declare_material_infos(solid_mechanics LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_documentation_files(solid_mechanics manual-solidmechanicsmodel.tex manual-constitutive-laws.tex manual-lumping.tex manual-appendix-materials.tex figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/static.pdf figures/static.svg figures/hooke_law.pdf figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/problemDomain.pdf_tex figures/problemDomain.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/stress_strain_neo.pdf figures/visco_elastic_law.pdf figures/isotropic_hardening_plasticity.pdf figures/stress_strain_visco.pdf ) package_declare_extra_files_to_package(solid_mechanics SOURCES model/solid_mechanics/material_list.hh.in ) diff --git a/src/model/solid_mechanics/materials/material_core_includes.hh b/src/model/solid_mechanics/materials/material_core_includes.hh index 490c60aab..9f4eaa05a 100644 --- a/src/model/solid_mechanics/materials/material_core_includes.hh +++ b/src/model/solid_mechanics/materials/material_core_includes.hh @@ -1,67 +1,71 @@ /** * @file material_core_includes.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Feb 03 2016 * * @brief List of materials for core package * * * 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_CORE_INCLUDES_HH__ #define __AKANTU_MATERIAL_CORE_INCLUDES_HH__ /* -------------------------------------------------------------------------- */ /* Material list */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_CMAKE_LIST_MATERIALS // elastic materials #include "material_elastic.hh" #include "material_elastic_linear_anisotropic.hh" #include "material_elastic_orthotropic.hh" #include "material_neohookean.hh" // visco-elastic materials #include "material_standard_linear_solid_deviatoric.hh" // damage laws #include "material_marigo.hh" #include "material_mazars.hh" // small-deformation plasticity #include "material_linear_isotropic_hardening.hh" +// Drucker-Prager plasticity +#include "material_drucker_prager.hh" + #endif #define AKANTU_CORE_MATERIAL_LIST \ ((2, (elastic, MaterialElastic)))((2, (neohookean, MaterialNeohookean)))( \ (2, (elastic_orthotropic, MaterialElasticOrthotropic)))( \ (2, (elastic_anisotropic, MaterialElasticLinearAnisotropic)))( \ (2, (sls_deviatoric, MaterialStandardLinearSolidDeviatoric)))( \ (2, (marigo, MaterialMarigo)))((2, (mazars, MaterialMazars)))( \ (2, (plastic_linear_isotropic_hardening, \ - MaterialLinearIsotropicHardening))) + MaterialLinearIsotropicHardening)))( \ + (2, (plastic_drucker_prager, MaterialDruckerPrager))) #endif /* __AKANTU_MATERIAL_CORE_INCLUDES_HH__ */ 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 new file mode 100644 index 000000000..a5922cf40 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc @@ -0,0 +1,192 @@ +/** + * @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"); + this->registerParam("fc", fc, Real(1.), _pat_parsable | _pat_modifiable, + "Compressive strength"); + + // compute alpha and k parameters for Drucker-Prager + this->alpha = (6.*sin(this->phi))/(3.-sin(this->phi)); + Real cohesion = this->fc * (1. - sin(this->phi))/(2.*cos(this->phi)); + this->k = (6. * cohesion * cos(this->phi))/(3. - sin(this->phi)); +} + +/* -------------------------------------------------------------------------- */ +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)) const 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 new file mode 100644 index 000000000..24d094869 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh @@ -0,0 +1,136 @@ +/** + * @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-pruger 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(const 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 void computeDeviatoricStress(const Matrix & sigma, + Matrix & sigma_dev); + +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-10); + + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +private: + // friction angle + Real phi; + + // compressive strength + Real fc; + + // + Real alpha; + + // + Real k; +}; + +} // 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 new file mode 100644 index 000000000..509925d7d --- /dev/null +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh @@ -0,0 +1,330 @@ +/** + * @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 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 = this->alpha * sigma.trace() - this->k; + + return sigma_dev_eff + modified_yield_stress; +} + +/* -------------------------------------------------------------------------- */ +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); + + // plastic multiplier guess at each iteration, initial guess is zero + plastic_multiplier_guess = 0.; + + // gradient of yield surface in voigt notation + gradient_f.clear(); + + // plastic strain increment at each iteration + delta_inelastic_strain.clear(); + + // 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; + + // lambda function to compute gradient of yield surface in voigt notation + auto compute_gradient_f = [&sigma_guess, &scaling_matrix, &kronecker_delta, + &gradient_f, this](){ + + const UInt dimension = sigma_guess.cols(); + + Matrix sigma_dev(dimension, dimension, 0); + this->computeDeviatoricStress(sigma_guess, sigma_dev); + + 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 += this->alpha * kronecker_delta; + }; + + // lambda function to compute hessian matrix of yield surface + auto compute_hessian_f = [&sigma_guess, &hessian_f, &scaling_matrix, + &kronecker_delta, this](){ + + const UInt dimension = sigma_guess.cols(); + + Matrix sigma_dev(dimension, dimension, 0); + this->computeDeviatoricStress(sigma_guess, sigma_dev); + + 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, this](){ + + // compute gradient + compute_gradient_f(); + + // compute yield function + yield_function = this->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 + auto error = std::max(f.norm(), std::abs(yield_function)); + return error; + }; + + auto projection_error = update_f(); + + /* --------------------------- */ + /* iteration loop */ + /* --------------------------- */ + UInt iterations{0}; + while(tolerance < projection_error and iterations < max_iterations) { + + // compute hessian at previous step + compute_hessian_f(); + + // compute inverse matrix Xi + Matrix xi(size, size); + xi = Ce + plastic_multiplier_guess * hessian_f; + + // compute inverse matrix Xi + Matrix xi_inv(size, size); + xi_inv.inverse(xi); + + Vector tmp(size); + tmp.mul(xi_inv, gradient_f); + + auto denominator = gradient_f.dot(tmp); + + // compute plastic multplier guess + Vector tmp1(size); + 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 + Matrix tmp2(size, size); + 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(); + iterations++; + } +} + +/* -------------------------------------------------------------------------- */ +/// 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 function + bool initial_yielding = + (this->computeYieldFunction(sigma_tr) < 0); + + // 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) { + + 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.); + + this->computeGradientAndPlasticMultplier(sigma_tr, dp, gradient_f, + delta_inelastic_strain_voigt); + + 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); +} + +/* -------------------------------------------------------------------------- */ +/// 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 { + +} + +}