diff --git a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc index 841e5c890..0de48b797 100644 --- a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc @@ -1,571 +1,576 @@ /** * @file material_viscoelastic_maxwell.hh * * @author Emil Gallyamov * * @date creation: Tue May 08 2018 * @date last modification: Tue May 08 2018 * * @brief Material Visco-elastic, based on Maxwell chain, * see * [] R. de Borst and A.H. van den Boogaard "Finite-element modeling of * deformation and cracking in early-age concrete", J.Eng.Mech., 1994 * as well as * [] Manual of DIANA FEA Theory manual v.10.2 Section 37.6 * * @section LICENSE * * 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 "material_viscoelastic_maxwell.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialViscoelasticMaxwell::MaterialViscoelasticMaxwell( SolidMechanicsModel & model, const ID & id) : MaterialElastic(model, id), C(voigt_h::size, voigt_h::size), sigma_v("sigma_v", *this), dissipated_energy("dissipated_energy", *this) { AKANTU_DEBUG_IN(); this->registerParam("Einf", Einf, Real(1.), _pat_parsable | _pat_modifiable, "Stiffness of the elastic element"); this->registerParam("previous_dt", previous_dt, Real(0.), _pat_readable, "Time step of previous solveStep"); this->registerParam("Eta", Eta, _pat_parsable | _pat_modifiable, "Viscosity of a Maxwell element"); this->registerParam("Ev", Ev, _pat_parsable | _pat_modifiable, "Stiffness of a Maxwell element"); this->use_previous_stress = true; this->use_previous_gradu = true; this->use_previous_stress_thermal = true; this->dissipated_energy.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::initMaterial() { AKANTU_DEBUG_IN(); // this->E = Einf + Ev1; this->E = std::min(this->Einf, this->Ev(0)); MaterialElastic::initMaterial(); AKANTU_DEBUG_ASSERT(this->Eta.size() == this->Ev.size(), "Eta and Ev have different dimensions! Please correct."); UInt stress_size = spatial_dimension * spatial_dimension; this->sigma_v.initialize(stress_size * this->Ev.size()); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell< spatial_dimension>::updateInternalParameters() { MaterialElastic::updateInternalParameters(); Real pre_mult = 1 / (1 + this->nu) / (1 - 2 * this->nu); UInt n = voigt_h::size; Real Miiii = pre_mult * (1 - this->nu); Real Miijj = pre_mult * this->nu; Real Mijij = pre_mult * 0.5 * (1 - 2 * this->nu); if (spatial_dimension == 1) { C(0, 0) = 1; } else { C(0, 0) = Miiii; } if (spatial_dimension >= 2) { C(1, 1) = Miiii; C(0, 1) = Miijj; C(1, 0) = Miijj; C(n - 1, n - 1) = Mijij; } if (spatial_dimension == 3) { C(2, 2) = Miiii; C(0, 2) = Miijj; C(1, 2) = Miijj; C(2, 0) = Miijj; C(2, 1) = Miijj; C(3, 3) = Mijij; C(4, 4) = Mijij; } } /* -------------------------------------------------------------------------- */ template <> void MaterialViscoelasticMaxwell<2>::updateInternalParameters() { MaterialElastic<2>::updateInternalParameters(); Real pre_mult = 1 / (1 + this->nu) / (1 - 2 * this->nu); UInt n = voigt_h::size; Real Miiii = pre_mult * (1 - this->nu); Real Miijj = pre_mult * this->nu; Real Mijij = pre_mult * 0.5 * (1 - 2 * this->nu); C(0, 0) = Miiii; C(1, 1) = Miiii; C(0, 1) = Miijj; C(1, 0) = Miijj; C(n - 1, n - 1) = Mijij; } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::setToSteadyState( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /* Array & stress_dev_vect = stress_dev(el_type, ghost_type); Array & history_int_vect = history_integral(el_type, ghost_type); Array::matrix_iterator stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator history_int = history_int_vect.begin(spatial_dimension, spatial_dimension); /// Loop on all quadrature points MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & dev_s = *stress_d; Matrix & h = *history_int; /// Compute the first invariant of strain Real Theta = grad_u.trace(); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { dev_s(i, j) = 2 * this->mu * (.5 * (grad_u(i, j) + grad_u(j, i)) - 1. / 3. * Theta * (i == j)); h(i, j) = 0.; } /// Save the deviator of stress ++stress_d; ++history_int; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; */ AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialThermal::computeStress(el_type, ghost_type); 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 sigma_v_it = this->sigma_v(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension, this->Eta.size()); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); auto & previous_grad_u = *previous_gradu_it; auto & previous_sigma = *previous_stress_it; computeStressOnQuad(grad_u, previous_grad_u, sigma, previous_sigma, *sigma_v_it, *sigma_th_it, *previous_sigma_th_it); ++sigma_th_it; ++previous_sigma_th_it; ++previous_stress_it; ++previous_gradu_it; ++sigma_v_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } + /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::computeStressOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Tensor3 & sigma_v, const Real & sigma_th, const Real & previous_sigma_th) { - Matrix delta_sigma(spatial_dimension, spatial_dimension); Matrix grad_delta_u(grad_u); grad_delta_u -= previous_grad_u; Real delta_sigma_th = sigma_th - previous_sigma_th; // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation - Vector voigt_delta_strain(voigt_h::size); - Vector voigt_delta_stress(voigt_h::size); + Vector voigt_current_strain(voigt_h::size); + Vector voigt_previous_strain(voigt_h::size); + Vector voigt_stress(voigt_h::size); Vector voigt_sigma_v(voigt_h::size); for (UInt I = 0; I < voigt_h::size; ++I) { Real voigt_factor = voigt_h::factors[I]; UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; - voigt_delta_strain(I) = - voigt_factor * (grad_delta_u(i, j) + grad_delta_u(j, i)) / 2.; + voigt_current_strain(I) = + voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; + voigt_previous_strain(I) = + voigt_factor * (previous_grad_u(i, j) + previous_grad_u(j, i)) / 2.; } - voigt_delta_stress = - this->Einf * this->C * voigt_delta_strain; - + voigt_stress = + this->Einf * this->C * voigt_current_strain; Real dt = this->model.getTimeStep(); for (UInt k = 0; k < Eta.size(); ++k) { Real lambda = this->Eta(k) / this->Ev(k); Real exp_dt_lambda = exp(-dt / lambda); Real E_additional = (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; for (UInt I = 0; I < voigt_h::size; ++I) { UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; voigt_sigma_v(I) = sigma_v(i, j, k); } - voigt_delta_stress =+ E_additional * this->C * voigt_delta_strain - (1 - exp_dt_lambda) * voigt_sigma_v; + voigt_stress += E_additional * this->C * + (voigt_current_strain - voigt_previous_strain) + + exp_dt_lambda * voigt_sigma_v; } for (UInt I = 0; I < voigt_h::size; ++I) { UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; - delta_sigma(i, j) = delta_sigma(j, i) = - voigt_delta_stress(I) + (i == j) * delta_sigma_th; + sigma(i, j) = sigma(j, i) = + voigt_stress(I) + (i == j) * sigma_th; } - - sigma = previous_sigma + delta_sigma; } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::afterSolveStep() { Material::afterSolveStep(); for (auto & el_type : this->element_filter.elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { auto previous_gradu_it = this->gradu.previous(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension); auto sigma_v_it = this->sigma_v(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension, this->Eta.size()); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); auto & previous_grad_u = *previous_gradu_it; updateIntVarOnQuad(grad_u, previous_grad_u, *sigma_v_it); ++previous_gradu_it; ++sigma_v_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; this->updateDissipatedEnergy(el_type, _not_ghost); } } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::updateIntVarOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Tensor3 & sigma_v) { Matrix grad_delta_u(grad_u); grad_delta_u -= previous_grad_u; Real dt = this->model.getTimeStep(); + Vector voigt_delta_strain(voigt_h::size); + for (UInt I = 0; I < voigt_h::size; ++I) { + Real voigt_factor = voigt_h::factors[I]; + UInt i = voigt_h::vec[I][0]; + UInt j = voigt_h::vec[I][1]; + + voigt_delta_strain(I) = + voigt_factor * (grad_delta_u(i, j) + grad_delta_u(j, i)) / 2.; + } + for (UInt k = 0; k < Eta.size(); ++k) { Real lambda = this->Eta(k) / this->Ev(k); Real exp_dt_lambda = exp(-dt / lambda); Real E_ef_v = (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; - // Wikipedia convention: - // 2*eps_ij (i!=j) = voigt_eps_I - // http://en.wikipedia.org/wiki/Voigt_notation - Vector voigt_delta_strain(voigt_h::size); Vector voigt_sigma_v(voigt_h::size); for (UInt I = 0; I < voigt_h::size; ++I) { Real voigt_factor = voigt_h::factors[I]; UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; - voigt_delta_strain(I) = - voigt_factor * (grad_delta_u(i, j) + grad_delta_u(j, i)) / 2.; - voigt_sigma_v(I) = sigma_v(i, j, k); } voigt_sigma_v = exp_dt_lambda * voigt_sigma_v + E_ef_v * this->C * voigt_delta_strain; for (UInt I = 0; I < voigt_h::size; ++I) { UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; sigma_v(i, j, k) = sigma_v(j, i, k) = voigt_sigma_v(I); } } } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::computeTangentModuli( const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real dt = this->model.getTimeStep(); Real E_ef = this->Einf; for (UInt k = 0; k < Eta.size(); ++k) { Real lambda = this->Eta(k) / this->Ev(k); Real exp_dt_lambda = exp(-dt / lambda); - E_ef =+ (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; + E_ef += (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; } this->previous_dt = dt; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); this->computeTangentModuliOnQuad(tangent); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; tangent_matrix *= E_ef; this->was_stiffness_assembled = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::computeTangentModuliOnQuad( Matrix & tangent) { tangent.copy(C); } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::savePreviousState() { for (auto & el_type : this->element_filter.elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { auto sigma_th_it = this->sigma_th(el_type, _not_ghost).begin(); auto previous_sigma_th_it = this->sigma_th.previous(el_type, _not_ghost).begin(); auto previous_gradu_it = this->gradu.previous(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension); auto previous_sigma_it = this->stress.previous(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension); auto sigma_v_it = this->sigma_v(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension, this->Eta.size()); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); auto & previous_grad_u = *previous_gradu_it; auto & previous_sigma = *previous_sigma_it; previous_grad_u.copy(*gradu_it); previous_sigma.copy(*stress_it); *previous_sigma_th_it = *sigma_th_it; ++previous_gradu_it, ++previous_sigma_it, ++previous_sigma_th_it, ++sigma_v_it, ++sigma_th_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::updateIntVariables() { for (auto & el_type : this->element_filter.elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { auto previous_gradu_it = this->gradu.previous(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension); auto previous_sigma_it = this->stress.previous(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension); auto sigma_v_it = this->sigma_v(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension, this->Eta.size()); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); auto & previous_grad_u = *previous_gradu_it; updateIntVarOnQuad(grad_u, previous_grad_u, *sigma_v_it); ++previous_gradu_it, ++sigma_v_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::updateDissipatedEnergy( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); // if(ghost_type == _ghost) return 0.; Real * dis_energy = dissipated_energy(el_type, ghost_type).storage(); /* Array & sigma_vec = stress_dev(el_type, ghost_type); Array & history_int_vect = history_integral(el_type, ghost_type); Array::matrix_iterator stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator history_int = history_int_vect.begin(spatial_dimension, spatial_dimension); Matrix q(spatial_dimension, spatial_dimension); Matrix q_rate(spatial_dimension, spatial_dimension); Matrix epsilon_d(spatial_dimension, spatial_dimension); Matrix epsilon_v(spatial_dimension, spatial_dimension); Real dt = this->model.getTimeStep(); Real gamma_v = Ev / this->E; Real alpha = 1. / (2. * this->mu * gamma_v); /// Loop on all quadrature points MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & dev_s = *stress_d; Matrix & h = *history_int; /// Compute the first invariant of strain this->template gradUToEpsilon(grad_u, epsilon_d); Real Theta = epsilon_d.trace(); epsilon_v.eye(Theta / Real(3.)); epsilon_d -= epsilon_v; q.copy(dev_s); q -= h; q *= gamma_v; q_rate.copy(dev_s); q_rate *= gamma_v; q_rate -= q; q_rate /= tau; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) *dis_energy += (epsilon_d(i, j) - alpha * q(i, j)) * q_rate(i, j) * dt; /// Save the deviator of stress ++stress_d; ++history_int; ++dis_energy; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; */ AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialViscoelasticMaxwell::getDissipatedEnergy() const { AKANTU_DEBUG_IN(); Real de = 0.; /// integrate the dissipated energy for each type of elements for (auto & type : this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { de += this->fem.integrate(dissipated_energy(type, _not_ghost), type, _not_ghost, this->element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return de; } /* -------------------------------------------------------------------------- */ template Real MaterialViscoelasticMaxwell::getDissipatedEnergy( ElementType type, UInt index) const { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = this->fem.getNbIntegrationPoints(type); auto it = this->dissipated_energy(type, _not_ghost).begin(nb_quadrature_points); UInt gindex = (this->element_filter(type, _not_ghost))(index); AKANTU_DEBUG_OUT(); return this->fem.integrate(it[index], type, gindex); } /* -------------------------------------------------------------------------- */ template Real MaterialViscoelasticMaxwell::getEnergy( const std::string & type) { if (type == "dissipated") return getDissipatedEnergy(); else if (type == "viscoelastic_maxwell") return getDissipatedEnergy(); else return MaterialElastic::getEnergy(type); } /* -------------------------------------------------------------------------- */ template Real MaterialViscoelasticMaxwell::getEnergy( const std::string & energy_id, ElementType type, UInt index) { if (energy_id == "dissipated") return getDissipatedEnergy(type, index); else if (energy_id == "viscoelastic_maxwell") return getDissipatedEnergy(type, index); else return MaterialElastic::getEnergy(energy_id, type, index); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(viscoelastic_maxwell, MaterialViscoelasticMaxwell); } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh index 65fafe6f0..f351d5098 100644 --- a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh @@ -1,193 +1,193 @@ /** * @file material_viscoelastic_maxwell.hh * * @author Emil Gallyamov * * @date creation: Tue May 08 2018 * @date last modification: Tue May 08 2018 * * @brief Material Visco-elastic, based on Maxwell chain, * see * [] R. de Borst and A.H. van den Boogaard "Finite-element modeling of * deformation and cracking in early-age concrete", J.Eng.Mech., 1994 * as well as * [] Manual of DIANA FEA Theory manual v.10.2 Section 37.6 * * @section LICENSE * * 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_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_VISCOELASTIC_MAXWELL_HH__ #define __AKANTU_MATERIAL_VISCOELASTIC_MAXWELL_HH__ namespace akantu { /** * Material Viscoelastic based on Maxwell chain * * * @verbatim E_0 ------|\/\/\|------- | | ---| |--- | | ----|\/\/\|--[|----- - | E_v1 \lambda_1| + | E_v1 \Eta_1| ---| |--- | | ----|\/\/\|--[|----- - | E_v2 \lambda_2 | + | E_v2 \Eta_2 | ---| |--- | | ----|\/\/\|--[|---- - E_vN \lambda_N + E_vN \Eta_N @endverbatim * * keyword : viscoelastic_maxwell * * parameters in the material files : * - N : number of Maxwell elements * - Einf : one spring element stiffness * - Ev1 : stiffness of the 1st viscous element - * - lambda1: relaxation time of the 1st Maxwell element + * - Eta1: viscosity of the 1st Maxwell element * ... * - Ev : stiffness of the Nst viscous element - * - lambda: relaxation time of the Nst Maxwell element + * - Eta: viscosity of the Nst Maxwell element */ template class MaterialViscoelasticMaxwell : public MaterialElastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialViscoelasticMaxwell(SolidMechanicsModel & model, const ID & id = ""); ~MaterialViscoelasticMaxwell() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter void initMaterial() override; /// recompute the lame coefficient and effective tangent moduli void updateInternalParameters() override; /// update internal variable on a converged Newton void afterSolveStep() override; /// update internal variable based on previous and current strain values void updateIntVariables(); /// update the internal variable sigma_v on quadrature point void updateIntVarOnQuad(const Matrix & grad_u, const Matrix & previous_grad_u, Tensor3 & sigma_v); /// set material to steady state void setToSteadyState(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// 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; /// save previous stress and strain values into "previous" arrays void savePreviousState() override; protected: /// update the dissipated energy, is called after the stress have been /// computed void updateDissipatedEnergy(ElementType el_type, GhostType ghost_type); /// compute stresses on a quadrature point void computeStressOnQuad(const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Tensor3 & sigma_v, const Real & sigma_th, const Real & previous_sigma_th); /// compute tangent moduli on a quadrature point void computeTangentModuliOnQuad(Matrix & tangent); bool hasStiffnessMatrixChanged() override { Real dt = this->model.getTimeStep(); return ((this->previous_dt == dt) ? (!(this->previous_dt == dt)) * (this->was_stiffness_assembled) : (!(this->previous_dt == dt))); // return (!(this->previous_dt == dt)); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// give the dissipated energy for the time step Real getDissipatedEnergy() const; Real getDissipatedEnergy(ElementType type, UInt index) const; /// get the energy using an energy type string for the time step Real getEnergy(const std::string & type) override; Real getEnergy(const std::string & energy_id, ElementType type, UInt index) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: using voigt_h = VoigtHelper; /// Vectors of viscosity, viscous elastic modulus, one spring element elastic modulus Vector Eta; Vector Ev; Real Einf; /// time step from previous solveStep Real previous_dt; /// Effective viscoelastic stiffness tensor in voigt notation Matrix C; /// Internal variable: viscous_stress InternalField sigma_v; /// Dissipated energy InternalField dissipated_energy; }; } // namespace akantu #endif /* __AKANTU_MATERIAL_VISCOELASTIC_MAXWELL_HH__ */