Homec4science
Diffusion akantu (bugfixes/36-pretty-print-of-akantu-array) Browse

rAKA/src/material_damage/material_orthotropic_damage_tmpl.hh
59c892f82354bugfixes/36-pretty…

Last Change
Raw File

material_orthotropic_damage_tmpl.hh
View Options

/**
* @file material_orthotropic_damage_tmpl.hh
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Sun Mar 8 12:54:30 2015
*
* @brief Specialization of the material class for the orthotropic
* damage material
*
* @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 "material_orthotropic_damage.hh"
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
MaterialOrthotropicDamage<spatial_dimension, Parent>::MaterialOrthotropicDamage(SolidMechanicsModel & model,
const ID & id) :
Material(model, id), Parent<spatial_dimension>(model, id),
damage("damage", *this),
dissipated_energy("damage dissipated energy", *this),
int_sigma("integral of sigma", *this),
damage_dir_vecs("damage_principal_directions", *this) {
AKANTU_DEBUG_IN();
this->registerParam("eta" , eta , 2. , _pat_parsable | _pat_modifiable, "Damage sensitivity parameter" );
this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value" );
this->is_non_local = false;
this->use_previous_stress = true;
this->use_previous_gradu = true;
/// use second order tensor for description of damage state
this->damage .initialize(spatial_dimension * spatial_dimension);
this->dissipated_energy.initialize(1);
this->int_sigma .initialize(1);
this->damage_dir_vecs.initialize(spatial_dimension * spatial_dimension);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
void MaterialOrthotropicDamage<spatial_dimension, Parent>::initMaterial() {
AKANTU_DEBUG_IN();
Parent<spatial_dimension>::initMaterial();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the dissipated energy in each element by a trapezoidal approximation
* of
* @f$ Ed = \int_0^{\epsilon}\sigma(\omega)d\omega - \frac{1}{2}\sigma:\epsilon@f$
*/
template<UInt spatial_dimension, template<UInt> class Parent>
void MaterialOrthotropicDamage<spatial_dimension, Parent>::updateEnergies(ElementType el_type, GhostType ghost_type) {
Parent<spatial_dimension>::updateEnergies(el_type, ghost_type);
this->computePotentialEnergy(el_type, ghost_type);
Array<Real>::matrix_iterator epsilon_p =
this->gradu.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator sigma_p =
this->stress.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension);
Array<Real>::const_scalar_iterator epot = this->potential_energy(el_type, ghost_type).begin();
Array<Real>::scalar_iterator ints = this->int_sigma(el_type, ghost_type).begin();
Array<Real>::scalar_iterator ed = this->dissipated_energy(el_type, ghost_type).begin();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
Matrix<Real> delta_gradu_it(*gradu_it);
delta_gradu_it -= *epsilon_p;
Matrix<Real> sigma_h(sigma);
sigma_h += *sigma_p;
Real dint = .5 * sigma_h.doubleDot(delta_gradu_it);
*ints += dint;
*ed = *ints - *epot;
++epsilon_p;
++sigma_p;
++epot;
++ints;
++ed;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
void MaterialOrthotropicDamage<spatial_dimension, Parent>::computeTangentModuli(const ElementType & el_type,
Array<Real> & tangent_matrix,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
Parent<spatial_dimension>::computeTangentModuli(el_type, tangent_matrix, ghost_type);
/// get the damage array for current element type
Array<Real> & dam = this->damage(el_type);
Array<Real>::matrix_iterator dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension);
/// get the directions of damage for the current element type
Array<Real> & dam_dirs = this->damage_dir_vecs(el_type);
Array<Real>::matrix_iterator damage_directions_it =
dam_dirs.begin(this->spatial_dimension,
this->spatial_dimension);
/// for the computation of the Cauchy stress the matrices (1-D) and
/// (1-D)^(1/2) are needed. For the formulation see Engineering
/// Damage Mechanics by Lemaitre and Desmorat.
Matrix<Real> one_minus_D(this->spatial_dimension, this->spatial_dimension);
Matrix<Real> sqrt_one_minus_D(this->spatial_dimension, this->spatial_dimension);
Matrix<Real> one_minus_D_rot(spatial_dimension, spatial_dimension);
Matrix<Real> sqrt_one_minus_D_rot(spatial_dimension, spatial_dimension);
Matrix<Real> rotation_tmp(spatial_dimension, spatial_dimension);
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix);
if ( !(Math::are_float_equal((*dam_it).trace(), 0)) )
computeTangentModuliOnQuad(tangent,
tangent,
*dam_it,
*damage_directions_it,
one_minus_D,
sqrt_one_minus_D,
one_minus_D_rot,
sqrt_one_minus_D_rot,
rotation_tmp);
++dam_it;
++damage_directions_it;
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
void MaterialOrthotropicDamage<spatial_dimension, Parent>
::computeTangentModuliOnQuad(Matrix<Real> & tangent,
const Matrix<Real> C,
const Matrix<Real> & dam,
const Matrix<Real> & dam_directions,
Matrix<Real> & O_prime,
Matrix<Real> & S_prime,
Matrix<Real> & O,
Matrix<Real> & S,
Matrix<Real> & rotation_tmp) {
/// effect of damage on stiffness matrix: See Ragueneau et al. 2008, p. 423, ep. 7
Real trace_D = dam.trace();
this->computeOneMinusD(O_prime, dam);
this->computeSqrtOneMinusD(O_prime, S_prime);
this->rotateIntoComputationFrame(O_prime, O, dam_directions, rotation_tmp);
this->rotateIntoComputationFrame(S_prime, S, dam_directions, rotation_tmp);
/// compute new stiffness matrix in damage coordinate system
if (spatial_dimension == 1)
tangent *= (1-dam(0, 0));
if (spatial_dimension == 2) {
Real min_val = std::min( (this->eta / spatial_dimension * trace_D),
this->max_damage );
/// first row
tangent(0, 0) = (C(0,0)*S(0,0)*S(0,0) + C(1,0)*S(0,1)*S(0,1) - (min_val/2. - 1./2)*(C(0,0) + C(1,0)) + (O(0,0)*(C(0,0)*O(0,0) + C(1,0)*O(1,1)))/(trace_D - 2.));
tangent(0, 1) = (C(0,1)*S(0,0)*S(0,0) + C(1,1)*S(0,1)*S(0,1) - (min_val/2. - 1./2)*(C(0,1) + C(1,1)) + (O(0,0)*(C(0,1)*O(0,0) + C(1,1)*O(1,1)))/(trace_D - 2.));
tangent(0, 2) = (2.*C(2,2)*S(0,0)*S(0,1) + (2.*C(2,2)*O(0,0)*O(0,1))/(trace_D - 2.));
/// second row
tangent(1, 0) = (C(0,0)*S(0,1)*S(0,1) + C(1,0)*S(1,1)*S(1,1) - (min_val/2. - 1./2)*(C(0,0) + C(1,0)) + (O(1,1)*(C(0,0)*O(0,0) + C(1,0)*O(1,1)))/(trace_D - 2.));
tangent(1, 1) = (C(0,1)*S(0,1)*S(0,1) + C(1,1)*S(1,1)*S(1,1) - (min_val/2. - 1./2)*(C(0,1) + C(1,1)) + (O(1,1)*(C(0,1)*O(0,0) + C(1,1)*O(1,1)))/(trace_D - 2.));
tangent(1, 2) = (2.*C(2,2)*S(0,1)*S(1,1) + (2.*C(2,2)*O(0,1)*O(1,1))/(trace_D - 2.));
/// third row
tangent(2, 0) = ((O(0,1)*(C(0,0)*O(0,0) + C(1,0)*O(1,1)))/(trace_D - 2.) + C(0,0)*S(0,0)*S(0,1) + C(1,0)*S(0,1)*S(1,1));
tangent(2,1) = ((O(0,1)*(C(0,1)*O(0,0) + C(1,1)*O(1,1)))/(trace_D - 2.) + C(0,1)*S(0,0)*S(0,1) + C(1,1)*S(0,1)*S(1,1));
tangent(2,2) = ((2.*C(2,2)*O(0,1)*O(0,1))/(trace_D - 2.) + C(2,2)*S(0,1)*S(0,1) + C(2,2)*S(0,0)*S(1,1));
// /// first row
// tangent(0, 0) = (C(0,0)*S(0,0)*S(0,0) + C(1,0)*S(0,1)*S(0,1) - ((eta_effective*trace_D)/4. - 1./2)*(C(0,0) + C(1,0)) + (O(0,0)*(C(0,0)*O(0,0) + C(1,0)*O(1,1)))/(trace_D - 2.));
// tangent(0, 1) = (C(0,1)*S(0,0)*S(0,0) + C(1,1)*S(0,1)*S(0,1) - ((eta_effective*trace_D)/4. - 1./2)*(C(0,1) + C(1,1)) + (O(0,0)*(C(0,1)*O(0,0) + C(1,1)*O(1,1)))/(trace_D - 2.));
// tangent(0, 2) = (2.*C(2,2)*S(0,0)*S(0,1) + (2.*C(2,2)*O(0,0)*O(0,1))/(trace_D - 2.));
// /// second row
// tangent(1, 0) = (C(0,0)*S(0,1)*S(0,1) + C(1,0)*S(1,1)*S(1,1) - ((eta_effective*trace_D)/4. - 1./2)*(C(0,0) + C(1,0)) + (O(1,1)*(C(0,0)*O(0,0) + C(1,0)*O(1,1)))/(trace_D - 2.));
// tangent(1, 1) = (C(0,1)*S(0,1)*S(0,1) + C(1,1)*S(1,1)*S(1,1) - ((eta_effective*trace_D)/4. - 1./2)*(C(0,1) + C(1,1)) + (O(1,1)*(C(0,1)*O(0,0) + C(1,1)*O(1,1)))/(trace_D - 2.));
// tangent(1, 2) = (2.*C(2,2)*S(0,1)*S(1,1) + (2.*C(2,2)*O(0,1)*O(1,1))/(trace_D - 2.));
// /// third row
// tangent(2, 0) = ((O(0,1)*(C(0,0)*O(0,0) + C(1,0)*O(1,1)))/(trace_D - 2.) + C(0,0)*S(0,0)*S(0,1) + C(1,0)*S(0,1)*S(1,1));
// tangent(2,1) = ((O(0,1)*(C(0,1)*O(0,0) + C(1,1)*O(1,1)))/(trace_D - 2.) + C(0,1)*S(0,0)*S(0,1) + C(1,1)*S(0,1)*S(1,1));
// tangent(2,2) = ((2.*C(2,2)*O(0,1)*O(0,1))/(trace_D - 2.) + C(2,2)*S(0,1)*S(0,1) + C(2,2)*S(0,0)*S(1,1));
}
if (spatial_dimension == 3) {
}
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
inline void
MaterialOrthotropicDamage<spatial_dimension, Parent>::computeDamageAndStressOnQuad(Matrix<Real> & sigma, Matrix<Real> & one_minus_D, Matrix<Real> & sqrt_one_minus_D, Matrix<Real> & damage, Matrix<Real> & first_term, Matrix<Real> & third_term) {
// /// set all terms to zero
// first_term.clear();
// third_term.clear();
// /// hydrostatic sensitivity parameter
// Real eta = 2.;
// /// Definition of Cauchy stress based on second order damage tensor:
// /// "Anisotropic damage modelling of biaxial behaviour and rupture
// /// of concrete strucutres", Ragueneau et al., 2008, Eq. 7
// first_term.mul<false, false>(sqrt_one_minus_D, sigma);
// first_term *= sqrt_one_minus_D;
// Real second_term = 0;
// for (UInt i = 0; i < this->spatial_dimension; ++i) {
// for (UInt j = 0; j < this->spatial_dimension; ++j)
// second_term += sigma(i,j) * one_minus_D(i,j);
// }
// second_term /= (this->spatial_dimension - damage.trace());
// for (UInt i = 0; i < this->spatial_dimension; ++i) {
// for (UInt j = 0; j < this->spatial_dimension; ++j)
// one_minus_D(i,j) *= second_term;
// }
// third_term.eye(1./this->spatial_dimension * sigma.trace() * (1 - eta/(this->spatial_dimension) * damage.trace()));
// sigma = first_term - one_minus_D + third_term;
/// hydrostatic sensitivity parameter
// Real eta = 2.;
/// Definition of Cauchy stress based on second order damage tensor:
/// "Anisotropic damage modelling of biaxial behaviour and rupture
/// of concrete strucutres", Ragueneau et al., 2008, Eq. 7
first_term.mul<false, false>(sqrt_one_minus_D, sigma);
first_term *= sqrt_one_minus_D;
Real second_term = 0;
for (UInt i = 0; i < this->spatial_dimension; ++i) {
for (UInt j = 0; j < this->spatial_dimension; ++j)
second_term += sigma(i,j) * one_minus_D(i,j);
}
second_term /= (this->spatial_dimension - damage.trace());
// for (UInt i = 0; i < this->spatial_dimension; ++i) {
// for (UInt j = 0; j < this->spatial_dimension; ++j)
// one_minus_D(i,j) *= second_term;
// }
one_minus_D *= second_term;
third_term.eye(1./this->spatial_dimension * sigma.trace() * (1 - eta/(this->spatial_dimension) * damage.trace()));
sigma.copy(first_term);
sigma -= one_minus_D;
sigma += third_term;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
inline void MaterialOrthotropicDamage<spatial_dimension, Parent>
::rotateIntoComputationFrame(const Matrix<Real> & to_rotate,
Matrix<Real> & rotated,
const Matrix<Real> & damage_directions,
Matrix<Real> & rotation_tmp) {
/// copy matrix
// rotated.copy(to_rotate);
// to_rotate.mul<false, true>(rotated, damage_directions);
// rotated.mul<false, false>(damage_directions, to_rotate);
rotation_tmp.mul<false, true>(to_rotate, damage_directions);
rotated.mul<false, false>(damage_directions, rotation_tmp);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
inline void MaterialOrthotropicDamage<spatial_dimension, Parent>
::rotateIntoNewFrame(const Matrix<Real> & to_rotate,
Matrix<Real> & rotated,
const Matrix<Real> & damage_directions,
Matrix<Real> & rotation_tmp) {
/// copy matrix
// rotated.copy(to_rotate);
// to_rotate.mul<false, false>(rotated, damage_directions);
// rotated.mul<true, false>(damage_directions, to_rotate);
rotation_tmp.mul<false, false>(to_rotate, damage_directions);
rotated.mul<true, false>(damage_directions, rotation_tmp);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
inline void MaterialOrthotropicDamage<spatial_dimension, Parent>
::computeOneMinusD(Matrix<Real> & one_minus_D, const Matrix<Real> & damage) {
/// compute one minus
one_minus_D.eye();
one_minus_D -= damage;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template<UInt> class Parent>
inline void
MaterialOrthotropicDamage<spatial_dimension, Parent>::computeSqrtOneMinusD(const Matrix<Real> & one_minus_D, Matrix<Real> & sqrt_one_minus_D) {
/// To compute (1-D)^1/2 we need to check that we are in the
/// principal coordinate system of the damage
#ifndef AKANTU_NDEBUG
for (UInt i = 0; i < this->spatial_dimension; ++i) {
for (UInt j = 0; j < this->spatial_dimension; ++j) {
if (i != j)
AKANTU_DEBUG_ASSERT(Math::are_float_equal(one_minus_D(i,j), 0),
"The damage tensor has off-diagonal parts");
}
}
#endif //AKANTU_NDEBUG
/// compute (1-D)^1/2
sqrt_one_minus_D.copy(one_minus_D);
for (UInt i = 0; i < this->spatial_dimension; ++i)
sqrt_one_minus_D(i,i) = std::sqrt(sqrt_one_minus_D(i,i));
}
/* -------------------------------------------------------------------------- */
// /* -------------------------------------------------------------------------- */
// template<UInt spatial_dimension, template<UInt> class Parent>
// Real MaterialOrthotropicDamage<spatial_dimension, Parent>::getDissipatedEnergy() const {
// AKANTU_DEBUG_IN();
// Real de = 0.;
// const Mesh & mesh = this->model->getFEEngine().getMesh();
// /// integrate the dissipated energy for each type of elements
// Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost);
// Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost);
// for(; it != end; ++it) {
// de += this->model->getFEEngine().integrate(dissipated_energy(*it, _not_ghost), *it,
// _not_ghost, this->element_filter(*it, _not_ghost));
// }
// AKANTU_DEBUG_OUT();
// return de;
// }
// /* -------------------------------------------------------------------------- */
// template<UInt spatial_dimension, template<UInt> class Parent>
// Real MaterialOrthotropicDamage<spatial_dimension, Parent>::getEnergy(std::string type) {
// if(type == "dissipated") return getDissipatedEnergy();
// else return Parent<spatial_dimension>::getEnergy(type);
// }
// /* -------------------------------------------------------------------------- */
__END_AKANTU__
c4science · Help