diff --git a/src/model/solid_mechanics/materials/material_elastic.hh b/src/model/solid_mechanics/materials/material_elastic.hh
index 977a6e6aa..a88779574 100644
--- a/src/model/solid_mechanics/materials/material_elastic.hh
+++ b/src/model/solid_mechanics/materials/material_elastic.hh
@@ -1,126 +1,127 @@
/**
* @file material_elastic.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Lucas Frerot <lucas.frerot@epfl.ch>
*
* @date Wed Aug 04 10:58:42 2010
*
* @brief Material isotropic elastic
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material_thermal.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_ELASTIC_HH__
#define __AKANTU_MATERIAL_ELASTIC_HH__
__BEGIN_AKANTU__
/**
* Material elastic isotropic
*
* parameters in the material files :
* - E : Young's modulus (default: 0)
* - nu : Poisson's ratio (default: 1/2)
* - Plane_Stress : if 0: plane strain, else: plane stress (default: 0)
*/
template<UInt spatial_dimension>
class MaterialElastic : public MaterialThermal<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialElastic(SolidMechanicsModel & model, const ID & id = "");
virtual ~MaterialElastic() {};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
virtual void initMaterial();
/// constitutive law for all element of a type
virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);
/// compute the tangent stiffness matrix for an element type
void computeTangentModuli(const ElementType & el_type,
Array<Real> & tangent_matrix,
GhostType ghost_type = _not_ghost);
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed() const;
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed() const;
protected:
/// constitutive law for a given quadrature point
inline void computeStressOnQuad(const Matrix<Real> & grad_u,
- Matrix<Real> & sigma, Real sigma_th = 0);
+ Matrix<Real> & sigma,
+ const Real sigma_th = 0);
/// compute the tangent stiffness matrix for an element
void computeTangentModuliOnQuad(Matrix<Real> & tangent);
/// recompute the lame coefficient if E or nu changes
virtual void updateInternalParameters();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the stable time step
inline Real getStableTimeStep(Real h, const Element & element);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// First Lamé coefficient
Real lambda;
/// Second Lamé coefficient (shear modulus)
Real mu;
/// Bulk modulus
Real kpa;
/// Plane stress or plane strain
bool plane_stress;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_elastic_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_ELASTIC_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc b/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc
index f14d3ca0a..9960dd5a7 100644
--- a/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc
@@ -1,105 +1,105 @@
/**
* @file material_elastic_inline_impl.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Lucas Frerot <lucas.frerot@epfl.ch>
*
* @date Wed Aug 04 10:58:42 2010
*
* @brief Implementation of the inline functions of the material elastic
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
inline void MaterialElastic<spatial_dimension>::computeStressOnQuad(const Matrix<Real> & grad_u,
Matrix<Real> & sigma,
- Real sigma_th) {
+ const Real sigma_th) {
Real trace = grad_u.trace();/// trace = (\nabla u)_{kk}
/// \sigma_{ij} = \lambda * (\nabla u)_{kk} * \delta_{ij} + \mu * (\nabla u_{ij} + \nabla u_{ji})
for (UInt i = 0; i < spatial_dimension; ++i) {
for (UInt j = 0; j < spatial_dimension; ++j) {
sigma(i, j) = (i == j)*lambda*trace + mu*(grad_u(i, j) + grad_u(j, i)) + (i == j) * sigma_th;
}
}
}
/* -------------------------------------------------------------------------- */
template<>
inline void MaterialElastic<1>::computeStressOnQuad(const Matrix<Real> & grad_u,
Matrix<Real> & sigma,
Real sigma_th) {
sigma(0, 0) = this->E * grad_u(0, 0) + sigma_th;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
inline void MaterialElastic<spatial_dimension>::computeTangentModuliOnQuad(Matrix<Real> & tangent) {
UInt n = tangent.cols();
//Real Ep = E/((1+nu)*(1-2*nu));
Real Miiii = lambda + 2*mu;
Real Miijj = lambda;
Real Mijij = mu;
if(spatial_dimension == 1)
tangent(0, 0) = this->E;
else
tangent(0, 0) = Miiii;
// test of dimension should by optimized out by the compiler due to the template
if(spatial_dimension >= 2) {
tangent(1, 1) = Miiii;
tangent(0, 1) = Miijj;
tangent(1, 0) = Miijj;
tangent(n - 1, n - 1) = Mijij;
}
if(spatial_dimension == 3) {
tangent(2, 2) = Miiii;
tangent(0, 2) = Miijj;
tangent(1, 2) = Miijj;
tangent(2, 0) = Miijj;
tangent(2, 1) = Miijj;
tangent(3, 3) = Mijij;
tangent(4, 4) = Mijij;
}
}
/* -------------------------------------------------------------------------- */
template<>
inline void MaterialElastic<1>::computeTangentModuliOnQuad(Matrix<Real> & tangent) {
tangent(0, 0) = E;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
inline Real MaterialElastic<spatial_dimension>::getStableTimeStep(Real h,
__attribute__ ((unused)) const Element & element) {
return (h/getPushWaveSpeed());
}
diff --git a/src/model/solid_mechanics/materials/material_plasticityinc/material_plasticityinc_inline_impl.cc b/src/model/solid_mechanics/materials/material_plasticityinc/material_plasticityinc_inline_impl.cc
index ecc0fc10e..a5821eaf7 100644
--- a/src/model/solid_mechanics/materials/material_plasticityinc/material_plasticityinc_inline_impl.cc
+++ b/src/model/solid_mechanics/materials/material_plasticityinc/material_plasticityinc_inline_impl.cc
@@ -1,165 +1,165 @@
/**
* @file material_plasticity_inc.cc
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Lucas Frerot <lucas.frerot@epfl.ch>
*
* @date Tue Jul 09 18:15:37 20130
*
* @brief Implementation of the inline functions of the material plasticity
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
#include <cmath>
#include "material_plasticityinc.hh"
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template<UInt dim>
inline void MaterialPlasticityinc<dim>::computeStressOnQuad(Matrix<Real> & grad_u,
Matrix<Real> & grad_delta_u,
Matrix<Real> & sigma,
Matrix<Real> & inelas_strain,
Real & iso_hardening,
Real sigma_th_cur,
Real sigma_th_prev) {
//Infinitesimal plasticity
Matrix<Real> sigma_tr(dim, dim);
Matrix<Real> sigma_tr_dev(dim, dim);
Matrix<Real> d_inelas_strain(dim,dim);
//Real r=iso_hardening;
Real dp=0.0;
Real d_dp=0.0;
UInt n=0;
Real delta_sigma_th = sigma_th_cur - sigma_th_prev;
//Compute trial stress, sigma_tr
- MaterialElastic<dim>::computeStressOnQuad(sigma, grad_delta_u, delta_sigma_th);
+ MaterialElastic<dim>::computeStressOnQuad(grad_delta_u, sigma_tr, delta_sigma_th);
sigma_tr += sigma;
// Compute deviatoric trial stress, sigma_tr_dev
sigma_tr_dev = sigma_tr;
sigma_tr_dev -= Matrix<Real>::eye(dim, sigma_tr.trace() / 3.0);
// Compute effective deviatoric trial stress
Real s = sigma_tr_dev.doubleDot(sigma_tr_dev);
Real sigma_tr_dev_eff = std::sqrt(3./2. * s);
//Loop for correcting stress based on yield function
while ((sigma_tr_dev_eff-iso_hardening-sigmay) > 0) {
//r = r + h * dp;
d_dp = (sigma_tr_dev_eff - 3. * this->mu *dp - iso_hardening - sigmay) / (3. * this->mu + h);
iso_hardening = iso_hardening + h * d_dp;
dp = dp + d_dp;
++n;
if ((d_dp < 1e-5) || (n>50))
break;
}
//Update internal variable
if (std::abs(sigma_tr_dev_eff) >
sigma_tr_dev.norm<L_inf>() * Math::getTolerance()) {
d_inelas_strain = sigma_tr_dev;
d_inelas_strain *= 3./2. * dp / sigma_tr_dev_eff;
}
- Matrix<Real> grad_u_elsatic(dim, dim);
- grad_u_elsatic = grad_delta_u;
- grad_u_elsatic -= d_inelas_strain;
+ Matrix<Real> grad_u_elastic(dim, dim);
+ grad_u_elastic = grad_delta_u;
+ grad_u_elastic -= d_inelas_strain;
Matrix<Real> sigma_elastic(dim, dim);
- MaterialElastic<dim>::computeStressOnQuad(sigma_elastic, grad_u_elsatic, delta_sigma_th);
+ MaterialElastic<dim>::computeStressOnQuad(grad_u_elastic, sigma_elastic, delta_sigma_th);
sigma += sigma_elastic;
inelas_strain += d_inelas_strain;
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
inline void MaterialPlasticityinc<dim>::computeTangentModuliOnQuad(Matrix<Real> & tangent,
Matrix<Real> & grad_delta_u,
Matrix<Real> & sigma_tensor,
Matrix<Real> & previous_sigma_tensor,
Real & iso_hardening) {
UInt cols = tangent.cols();
UInt rows = tangent.rows();
Matrix<Real> sigma_dev(dim, dim);
Matrix<Real> sigma_tr(dim, dim);
Matrix<Real> sigma_tr_dev(dim, dim);
// Real r=iso_hardening;
//Compute trial stress, sigma_tr
- MaterialElastic<dim>::computeStressOnQuad(sigma_tr, grad_delta_u);
+ MaterialElastic<dim>::computeStressOnQuad(grad_delta_u, sigma_tr);
sigma_tr += previous_sigma_tensor;
// Compute deviatoric trial stress, sigma_tr_dev
sigma_tr_dev = sigma_tr;
sigma_tr_dev -= Matrix<Real>::eye(dim, sigma_tr.trace() / 3.0);
// Compute effective deviatoric trial stress
Real s = sigma_tr_dev.doubleDot(sigma_tr_dev);
Real sigma_tr_dev_eff=std::sqrt(3./2. * s);
// Compute deviatoric stress, sigma_dev
sigma_dev = sigma_tensor;
sigma_dev -= Matrix<Real>::eye(dim, sigma_tensor.trace() / 3.0);
// Compute effective deviatoric stress
s = sigma_dev.doubleDot(sigma_dev);
Real sigma_dev_eff = std::sqrt(3./2. * s);
Real xr = 0.0;
if(sigma_tr_dev_eff > sigma_dev_eff * Math::getTolerance())
xr = sigma_dev_eff / sigma_tr_dev_eff;
Real q = 1.5 * (1. / (1. + 3. * this->mu / h) - xr);
for (UInt m = 0; m < rows; ++m) {
UInt i = VoigtHelper<dim>::vec[m][0];
UInt j = VoigtHelper<dim>::vec[m][1];
for (UInt n = 0; n < cols; ++n) {
- UInt k = VoigtHelper<dim>::vec[m][0];
- UInt l = VoigtHelper<dim>::vec[m][1];
+ UInt k = VoigtHelper<dim>::vec[n][0];
+ UInt l = VoigtHelper<dim>::vec[n][1];
if (((sigma_tr_dev_eff-iso_hardening-sigmay) > 0) && (xr > 0)) {
- tangent(m,n) =
- 2. * this->mu * q * (sigma_tr_dev (i,j) / sigma_tr_dev_eff) * (sigma_tr_dev (k,l) / sigma_tr_dev_eff) +
- (i==k) * (j==l) * 2. * this->mu * xr +
- (i==j) * (k==l) * (this->kpa - 2./3. * this->mu * xr);
- if ((m == n) && (m>=dim))
- tangent(m, n) = tangent(m, n) - this->mu * xr;
+ tangent(m,n) =
+ 2. * this->mu * q * (sigma_tr_dev (i,j) / sigma_tr_dev_eff) * (sigma_tr_dev (k,l) / sigma_tr_dev_eff) +
+ (i==k) * (j==l) * 2. * this->mu * xr +
+ (i==j) * (k==l) * (this->kpa - 2./3. * this->mu * xr);
+ if ((m == n) && (m>=dim))
+ tangent(m, n) = tangent(m, n) - this->mu * xr;
} else {
- tangent(m,n) =
- (i==k) * (j==l) * 2. * this->mu +
- (i==j) * (k==l) * this->lambda;
- if ((m==n) && (m>=dim))
- tangent(m,n) = tangent(m,n) - this->mu;
+ tangent(m,n) =
+ (i==k) * (j==l) * 2. * this->mu +
+ (i==j) * (k==l) * this->lambda;
+ if ((m==n) && (m>=dim))
+ tangent(m,n) = tangent(m,n) - this->mu;
}
//correct tangent stiffness for shear component
}
}
}