Page MenuHomec4science
Files F63959794

material_cohesive_linear.cc
No OneTemporary
Actions

File Metadata

Created
Thu, May 23, 15:27

material_cohesive_linear.cc
View Options

/**
* @file material_cohesive_linear.cc
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue May 08 2012
* @date last modification: Thu Aug 07 2014
*
* @brief Linear irreversible cohesive law of mixed mode loading with
* random stress definition for extrinsic type
*
* @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 <numeric>
/* -------------------------------------------------------------------------- */
#include "material_cohesive_linear.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "sparse_matrix.hh"
#include "dof_synchronizer.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
MaterialCohesiveLinear<spatial_dimension>::MaterialCohesiveLinear(SolidMechanicsModel & model,
const ID & id) :
MaterialCohesive(model,id),
sigma_c_eff("sigma_c_eff", *this),
delta_c_eff("delta_c_eff", *this),
insertion_stress("insertion_stress", *this) {
AKANTU_DEBUG_IN();
this->registerParam("beta" , beta , 0. ,
_pat_parsable | _pat_readable,
"Beta parameter" );
this->registerParam("G_c" , G_c , 0. ,
_pat_parsable | _pat_readable,
"Mode I fracture energy" );
this->registerParam("penalty", penalty, 0. ,
_pat_parsable | _pat_readable,
"Penalty coefficient" );
this->registerParam("volume_s", volume_s, 0. ,
_pat_parsable | _pat_readable,
"Reference volume for sigma_c scaling");
this->registerParam("m_s", m_s, 1. ,
_pat_parsable | _pat_readable,
"Weibull exponent for sigma_c scaling");
this->registerParam("kappa" , kappa , 1. ,
_pat_parsable | _pat_readable,
"Kappa parameter");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
MaterialCohesive::initMaterial();
/// compute scalars
beta2_kappa2 = beta * beta/kappa/kappa;
beta2_kappa = beta * beta/kappa;
if (Math::are_float_equal(beta, 0))
beta2_inv = 0;
else
beta2_inv = 1./beta/beta;
sigma_c_eff.initialize(1);
delta_c_eff.initialize(1);
insertion_stress.initialize(spatial_dimension);
if (!Math::are_float_equal(delta_c, 0.))
delta_c_eff.setDefaultValue(delta_c);
else
delta_c_eff.setDefaultValue(2 * G_c / sigma_c);
if (model->getIsExtrinsic()) scaleInsertionTraction();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::scaleInsertionTraction() {
AKANTU_DEBUG_IN();
// do nothing if volume_s hasn't been specified by the user
if (Math::are_float_equal(volume_s, 0.)) return;
const Mesh & mesh_facets = model->getMeshFacets();
const FEEngine & fe_engine = model->getFEEngine();
const FEEngine & fe_engine_facet = model->getFEEngine("FacetsFEEngine");
// loop over facet type
Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1);
Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1);
Real base_sigma_c = sigma_c;
for(;first != last; ++first) {
ElementType type_facet = *first;
const Array< std::vector<Element> > & facet_to_element
= mesh_facets.getElementToSubelement(type_facet);
UInt nb_facet = facet_to_element.getSize();
UInt nb_quad_per_facet = fe_engine_facet.getNbQuadraturePoints(type_facet);
// iterator to modify sigma_c for all the quadrature points of a facet
Array<Real>::vector_iterator sigma_c_iterator
= sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet);
for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) {
const std::vector<Element> & element_list = facet_to_element(f);
// compute bounding volume
Real volume = 0;
std::vector<Element>::const_iterator elem = element_list.begin();
std::vector<Element>::const_iterator elem_end = element_list.end();
for (; elem != elem_end; ++elem) {
if (*elem == ElementNull) continue;
// unit vector for integration in order to obtain the volume
UInt nb_quadrature_points = fe_engine.getNbQuadraturePoints(elem->type);
Vector<Real> unit_vector(nb_quadrature_points, 1);
volume += fe_engine.integrate(unit_vector, elem->type,
elem->element, elem->ghost_type);
}
// scale sigma_c
*sigma_c_iterator -= base_sigma_c;
*sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s);
*sigma_c_iterator += base_sigma_c;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::checkInsertion() {
AKANTU_DEBUG_IN();
const Mesh & mesh_facets = model->getMeshFacets();
CohesiveElementInserter & inserter = model->getElementInserter();
Real tolerance = Math::getTolerance();
Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1);
Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1);
for (; it != last; ++it) {
ElementType type_facet = *it;
ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet);
const Array<bool> & facets_check = inserter.getCheckFacets(type_facet);
Array<bool> & f_insertion = inserter.getInsertionFacets(type_facet);
Array<UInt> & f_filter = facet_filter(type_facet);
Array<Real> & sig_c_eff = sigma_c_eff(type_cohesive);
Array<Real> & del_c = delta_c_eff(type_cohesive);
Array<Real> & ins_stress = insertion_stress(type_cohesive);
Array<Real> & trac_old = tractions_old(type_cohesive);
const Array<Real> & f_stress = model->getStressOnFacets(type_facet);
const Array<Real> & sigma_lim = sigma_c(type_facet);
UInt nb_quad_facet = model->getFacetFEEngine().getNbQuadraturePoints(type_facet);
UInt nb_facet = f_filter.getSize();
if (nb_facet == 0) continue;
Array<Real>::const_iterator<Real> sigma_lim_it = sigma_lim.begin();
Matrix<Real> stress_tmp(spatial_dimension, spatial_dimension);
Matrix<Real> normal_traction(spatial_dimension, nb_quad_facet);
Vector<Real> stress_check(nb_quad_facet);
UInt sp2 = spatial_dimension * spatial_dimension;
const Array<Real> & tangents = model->getTangents(type_facet);
const Array<Real> & normals
= model->getFacetFEEngine().getNormalsOnQuadPoints(type_facet);
Array<Real>::const_vector_iterator normal_begin = normals.begin(spatial_dimension);
Array<Real>::const_vector_iterator tangent_begin = tangents.begin(tangents.getNbComponent());
Array<Real>::const_matrix_iterator facet_stress_begin =
f_stress.begin(spatial_dimension, spatial_dimension * 2);
std::list<Real> new_sigmas;
std::list< Vector<Real> > new_normal_traction;
std::list<Real> new_delta_c;
// loop over each facet belonging to this material
for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) {
UInt facet = f_filter(f);
// skip facets where check shouldn't be realized
if (!facets_check(facet)) continue;
// compute the effective norm on each quadrature point of the facet
for (UInt q = 0; q < nb_quad_facet; ++q) {
UInt current_quad = f * nb_quad_facet + q;
const Vector<Real> & normal = normal_begin[current_quad];
const Vector<Real> & tangent = tangent_begin[current_quad];
const Matrix<Real> & facet_stress_it = facet_stress_begin[current_quad];
// compute average stress on the current quadrature point
Matrix<Real> stress_1(facet_stress_it.storage(),
spatial_dimension,
spatial_dimension);
Matrix<Real> stress_2(facet_stress_it.storage() + sp2,
spatial_dimension,
spatial_dimension);
stress_tmp.copy(stress_1);
stress_tmp += stress_2;
stress_tmp /= 2.;
Vector<Real> normal_traction_vec(normal_traction(q));
// compute normal and effective stress
stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent,
normal_traction_vec);
}
// verify if the effective stress overcomes the threshold
if (stress_check.mean() > (*sigma_lim_it - tolerance)) {
f_insertion(facet) = true;
// store the new cohesive material parameters for each quadrature point
for (UInt q = 0; q < nb_quad_facet; ++q) {
Real new_sigma = stress_check(q);
Vector<Real> normal_traction_vec(normal_traction(q));
if (spatial_dimension != 3)
normal_traction_vec *= -1.;
new_sigmas.push_back(new_sigma);
new_normal_traction.push_back(normal_traction_vec);
Real new_delta;
// set delta_c in function of G_c or a given delta_c value
if (Math::are_float_equal(delta_c, 0.))
new_delta = 2 * G_c / new_sigma;
else
new_delta = (*sigma_lim_it) / new_sigma * delta_c;
new_delta_c.push_back(new_delta);
}
}
}
// update material data for the new elements
UInt old_nb_quad_points = sig_c_eff.getSize();
UInt new_nb_quad_points = new_sigmas.size();
sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points);
ins_stress.resize(old_nb_quad_points + new_nb_quad_points);
trac_old.resize(old_nb_quad_points + new_nb_quad_points);
del_c.resize(old_nb_quad_points + new_nb_quad_points);
std::list<Real>::iterator new_sig_it = new_sigmas.begin();
std::list<Real>::iterator new_del_it = new_delta_c.begin();
std::list< Vector<Real> >::iterator new_normal_trac_it = new_normal_traction.begin();
for (UInt q = 0; q < new_nb_quad_points; ++q, ++new_sig_it, ++new_del_it, ++new_normal_trac_it) {
sig_c_eff(old_nb_quad_points + q) = *new_sig_it;
del_c(old_nb_quad_points + q) = *new_del_it;
for (UInt dim = 0; dim < spatial_dimension; ++dim) {
ins_stress(old_nb_quad_points + q, dim) = (*new_normal_trac_it)(dim);
trac_old(old_nb_quad_points + q, dim) = (*new_normal_trac_it)(dim);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
inline Real MaterialCohesiveLinear<spatial_dimension>::computeEffectiveNorm(const Matrix<Real> & stress,
const Vector<Real> & normal,
const Vector<Real> & tangent,
Vector<Real> & normal_traction) {
normal_traction.mul<false>(stress, normal);
Real normal_contrib = normal_traction.dot(normal);
/// in 3D tangential components must be summed
Real tangent_contrib = 0;
if (spatial_dimension == 2) {
Real tangent_contrib_tmp = normal_traction.dot(tangent);
tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp;
}
else if (spatial_dimension == 3) {
for (UInt s = 0; s < spatial_dimension - 1; ++s) {
const Vector<Real> tangent_v(tangent.storage() + s * spatial_dimension,
spatial_dimension);
Real tangent_contrib_tmp = normal_traction.dot(tangent_v);
tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp;
}
}
tangent_contrib = std::sqrt(tangent_contrib);
normal_contrib = std::max(0., normal_contrib);
return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::computeTraction(const Array<Real> & normal,
ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
/// define iterators
Array<Real>::iterator< Vector<Real> > traction_it =
tractions(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > opening_it =
opening(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > contact_traction_it =
contact_tractions(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > contact_opening_it =
contact_opening(el_type, ghost_type).begin(spatial_dimension);
Array<Real>::const_iterator< Vector<Real> > normal_it =
normal.begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> >traction_end =
tractions(el_type, ghost_type).end(spatial_dimension);
Array<Real>::iterator<Real>sigma_c_it =
sigma_c_eff(el_type, ghost_type).begin();
Array<Real>::iterator<Real>delta_max_it =
delta_max(el_type, ghost_type).begin();
Array<Real>::iterator<Real>delta_c_it =
delta_c_eff(el_type, ghost_type).begin();
Array<Real>::iterator<Real>damage_it =
damage(el_type, ghost_type).begin();
Array<Real>::iterator<Vector<Real> > insertion_stress_it =
insertion_stress(el_type, ghost_type).begin(spatial_dimension);
Real * memory_space = new Real[2*spatial_dimension];
Vector<Real> normal_opening(memory_space, spatial_dimension);
Vector<Real> tangential_opening(memory_space + spatial_dimension,
spatial_dimension);
/// loop on each quadrature point
for (; traction_it != traction_end;
++traction_it, ++opening_it, ++normal_it, ++sigma_c_it,
++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it,
++insertion_stress_it, ++contact_opening_it) {
/// compute normal and tangential opening vectors
Real normal_opening_norm = opening_it->dot(*normal_it);
normal_opening = (*normal_it);
normal_opening *= normal_opening_norm;
tangential_opening = *opening_it;
tangential_opening -= normal_opening;
Real tangential_opening_norm = tangential_opening.norm();
/**
* compute effective opening displacement
* @f$ \delta = \sqrt{
* \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$
*/
Real delta = tangential_opening_norm * tangential_opening_norm * beta2_kappa2;
bool penetration = normal_opening_norm < -Math::getTolerance();
if (penetration) {
/// use penalty coefficient in case of penetration
*contact_traction_it = normal_opening;
*contact_traction_it *= penalty;
*contact_opening_it = normal_opening;
/// don't consider penetration contribution for delta
*opening_it = tangential_opening;
normal_opening.clear();
}
else {
delta += normal_opening_norm * normal_opening_norm;
contact_traction_it->clear();
contact_opening_it->clear();
}
delta = std::sqrt(delta);
/// update maximum displacement and damage
*delta_max_it = std::max(*delta_max_it, delta);
*damage_it = std::min(*delta_max_it / *delta_c_it, 1.);
/**
* Compute traction @f$ \mathbf{T} = \left(
* \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n
* \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1-
* \frac{\delta}{\delta_c} \right)@f$
*/
if (Math::are_float_equal(*damage_it, 1.))
traction_it->clear();
else if (Math::are_float_equal(*damage_it, 0.)) {
if (penetration)
traction_it->clear();
else
*traction_it = *insertion_stress_it;
}
else {
*traction_it = tangential_opening;
*traction_it *= beta2_kappa;
*traction_it += normal_opening;
AKANTU_DEBUG_ASSERT(*delta_max_it != 0.,
"Division by zero, tolerance might be too low");
*traction_it *= *sigma_c_it / *delta_max_it * (1. - *damage_it);
}
}
delete [] memory_space;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialCohesiveLinear);
__END_AKANTU__

Event Timeline

c4science · Help