Homec4science
Diffusion akantu (bugfix/67-change-in-stl-headers) Browse

rAKA/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc
213fa2e405ebbugfix/67-change-in-stl…

Last Change
Raw File

material_cohesive_exponential.cc
View Options

/**
* @file material_cohesive_exponential.cc
*
* @author Seyedeh Mohadeseh Taheri Mousavi <mohadeseh.taherimousavi@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Mon Jul 09 2012
* @date last modification: Fri Mar 21 2014
*
* @brief Exponential irreversible cohesive law of mixed mode loading
*
* @section LICENSE
*
* Copyright (©) 2014 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_cohesive_exponential.hh"
#include "solid_mechanics_model.hh"
#include "sparse_matrix.hh"
#include "dof_synchronizer.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
MaterialCohesiveExponential<spatial_dimension>::MaterialCohesiveExponential(SolidMechanicsModel & model,
const ID & id) :
MaterialCohesive(model,id) {
AKANTU_DEBUG_IN();
this->registerParam("beta" , beta , 0. , _pat_parsable, "Beta parameter" );
this->registerParam("exponential_penalty", exp_penalty, true,
_pat_parsable, "Is contact penalty following the exponential law?" );
this->registerParam("contact_tangent", contact_tangent, 1.0,
_pat_parsable, "Ratio of contact tangent over the initial exponential tangent" );
// this->initInternalArray(delta_max, 1, _ek_cohesive);
use_previous_delta_max=true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveExponential<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
MaterialCohesive::initMaterial();
if ((exp_penalty)&&(contact_tangent!=1)){
contact_tangent = 1;
AKANTU_DEBUG_WARNING("The parsed paramter <contact_tangent> is forced to 1.0 when the contact penalty follows the exponential law");
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveExponential<spatial_dimension>::computeTraction(const Array<Real> & normal,
ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
/// define iterators
Array<Real>::vector_iterator traction_it =
tractions(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::vector_iterator opening_it =
opening(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::const_vector_iterator normal_it =
normal.begin(spatial_dimension);
Array<Real>::vector_iterator traction_end =
tractions(el_type, ghost_type).end(spatial_dimension);
Array<Real>::iterator<Real> delta_max_it =
delta_max(el_type, ghost_type).begin();
Array<Real>::iterator<Real> delta_max_prev_it =
delta_max.previous(el_type, ghost_type).begin();
/// compute scalars
Real beta2 = beta*beta;
/// loop on each quadrature point
for (; traction_it != traction_end;
++traction_it, ++opening_it, ++normal_it, ++delta_max_it, ++delta_max_prev_it) {
/// compute normal and tangential opening vectors
Real normal_opening_norm = opening_it->dot(*normal_it);
Vector<Real> normal_opening(spatial_dimension);
normal_opening = (*normal_it);
normal_opening *= normal_opening_norm;
Vector<Real> tangential_opening(spatial_dimension);
tangential_opening = *opening_it;
tangential_opening -= normal_opening;
Real tangential_opening_norm = tangential_opening.norm();
/**
* compute effective opening displacement
* @f$ \delta = \sqrt{
* \beta^2 \Delta_t^2 + \Delta_n^2 } @f$
*/
Real delta = tangential_opening_norm;
delta *= delta * beta2;
delta += normal_opening_norm * normal_opening_norm;
delta = sqrt(delta);
if ((normal_opening_norm<0) && (std::abs(normal_opening_norm) > Math::getTolerance())) {
Vector<Real> op_n(*opening_it);
op_n *= *normal_it;
Vector<Real> delta_s(*opening_it);
delta_s -= op_n;
delta = tangential_opening_norm*beta;
computeCoupledTraction(*traction_it, *normal_it, delta, delta_s,
*delta_max_it, *delta_max_prev_it);
computeCompressiveTraction(*traction_it, *normal_it, normal_opening_norm,
*opening_it);
}
else computeCoupledTraction(*traction_it, *normal_it, delta, *opening_it,
*delta_max_it, *delta_max_prev_it);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveExponential<spatial_dimension>::computeCoupledTraction(Vector<Real> & tract,
const Vector<Real> & normal,
Real delta,
const Vector<Real> & opening,
Real & delta_max_new,
Real delta_max) {
AKANTU_DEBUG_IN();
/// full damage case
if (std::abs(delta) < Math::getTolerance()) {
/// set traction to zero
tract.clear();
} else { /// element not fully damaged
/**
* Compute traction loading @f$ \mathbf{T} =
* e \sigma_c \frac{\delta}{\delta_c} e^{-\delta/ \delta_c}@f$
*/
/**
* Compute traction unloading @f$ \mathbf{T} =
* \frac{t_{max}}{\delta_{max}} \delta @f$
*/
Real beta2 = beta*beta;
Vector<Real> op_n_n(opening);
op_n_n*=normal;
op_n_n*=(1-beta2);
op_n_n*=normal;
tract = beta2*opening;
tract += op_n_n;
delta_max_new = std::max(delta_max, delta);
tract *= exp(1)*sigma_c*exp(- delta_max_new / delta_c)/delta_c;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveExponential<spatial_dimension>::computeCompressiveTraction(Vector<Real> & tract,
const Vector<Real> & normal,
Real delta_n,
const Vector<Real> & opening) {
Vector<Real> temp_tract(normal);
if(exp_penalty) {
temp_tract *= delta_n*exp(1)*sigma_c*exp(- delta_n / delta_c)/delta_c;
}
else {
Real initial_tg = contact_tangent*exp(1)*sigma_c*delta_n/delta_c;
temp_tract *= initial_tg;
}
tract += temp_tract;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveExponential<spatial_dimension>::computeTangentTraction(const ElementType & el_type,
Array<Real> & tangent_matrix,
const Array<Real> & normal,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array<Real>::matrix_iterator tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension);
Array<Real>::const_vector_iterator normal_it = normal.begin(spatial_dimension);
Array<Real>::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::vector_iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator<Real>delta_max_it = delta_max.previous(el_type, ghost_type).begin();
Real beta2 = beta*beta;
/**
* compute tangent matrix @f$ \frac{\partial \mathbf{t}}
* {\partial \delta} = \hat{\mathbf{t}} \otimes
* \frac{\partial (t/\delta)}{\partial \delta}
* \frac{\hat{\mathbf{t}}}{\delta}+ \frac{t}{\delta} [ \beta^2 \mathbf{I} +
* (1-\beta^2) (\mathbf{n} \otimes \mathbf{n})] @f$
**/
/**
* In which @f$
* \frac{\partial(t/ \delta)}{\partial \delta} =
* \left\{\begin{array} {l l}
* -e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad if
* \delta \geq \delta_{max} \\
* 0 & \quad if \delta < \delta_{max}, \delta_n > 0
* \end{array}\right. @f$
**/
for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++traction_it, ++ delta_max_it) {
Real normal_opening_norm = opening_it->dot(*normal_it);
Vector<Real> normal_opening(spatial_dimension);
normal_opening = (*normal_it);
normal_opening *= normal_opening_norm;
Vector<Real> tangential_opening(spatial_dimension);
tangential_opening = *opening_it;
tangential_opening -= normal_opening;
Real tangential_opening_norm = tangential_opening.norm();
Real delta = tangential_opening_norm;
delta *= delta * beta2;
delta += normal_opening_norm * normal_opening_norm;
delta = sqrt(delta);
if ((normal_opening_norm<0) && (std::abs(normal_opening_norm) > Math::getTolerance())) {
Vector<Real> op_n(*opening_it);
op_n *= *normal_it;
Vector<Real> delta_s(*opening_it);
delta_s -= op_n;
delta = tangential_opening_norm*beta;
computeCoupledTangent(*tangent_it, *traction_it, *normal_it, delta,
delta_s, *delta_max_it);
computeCompressivePenalty(*tangent_it, *normal_it, normal_opening_norm);
}
else computeCoupledTangent(*tangent_it, *traction_it, *normal_it, delta,
*opening_it, *delta_max_it);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveExponential<spatial_dimension>::computeCoupledTangent(Matrix<Real> & tangent,
Vector<Real> tract,
const Vector<Real> & normal,
Real delta,
const Vector<Real> & opening,
Real delta_max_new) {
AKANTU_DEBUG_IN();
Real beta2 = beta*beta;
Matrix<Real> I(spatial_dimension, spatial_dimension);
I.eye(beta2);
Matrix<Real> nn(spatial_dimension, spatial_dimension);
nn.outerProduct(normal, normal);
nn *= (1-beta2);
I += nn;
if(std::abs(delta) < Math::getTolerance()){
I *= exp(1)*sigma_c/delta_c;
tangent += I;
} else {
Real traction_norm = tract.norm();
Vector<Real> t_hat(opening);
Vector<Real> beta2delta(opening);
beta2delta *= beta2;
t_hat *= normal;
t_hat *= (beta2-1);
t_hat *= normal;
t_hat += beta2delta;
Vector<Real> deriv_t_hat(t_hat);
deriv_t_hat /= delta;
Real derivatives = 0;
if ((delta > delta_max_new) || (std::abs(delta - delta_max_new) < 1e-12)) {
derivatives = -exp(1- delta/delta_c)*sigma_c/(delta_c*delta_c);
}
deriv_t_hat *= derivatives;
Matrix<Real> t_var_t(spatial_dimension, spatial_dimension);
t_var_t.outerProduct(t_hat, deriv_t_hat);
I *= traction_norm / delta;
tangent += t_var_t;
tangent += I;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveExponential<spatial_dimension>::computeCompressivePenalty(Matrix<Real> & tangent,
const Vector<Real> & normal,
Real delta_n) {
if (!exp_penalty) delta_n=0;
Matrix<Real> n_outer_n(spatial_dimension,spatial_dimension);
n_outer_n.outerProduct(normal,normal);
Real normal_tg = contact_tangent*exp(1)*sigma_c*exp(-delta_n/delta_c)*(1-delta_n/delta_c)/delta_c;
n_outer_n *= normal_tg;
tangent += n_outer_n;
}
INSTANSIATE_MATERIAL(MaterialCohesiveExponential);
__END_AKANTU__
c4science · Help