diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc index 4c9d37833..d97dd491a 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc @@ -1,313 +1,324 @@ /** * @file material_cohesive_linear_inline_impl.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Wed Apr 22 2015 * @date last modification: Thu Jan 14 2016 * * @brief Inline functions of the MaterialCohesiveLinear * * @section LICENSE * * Copyright (©) 2015 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template inline Real MaterialCohesiveLinear::computeEffectiveNorm(const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_traction) const { normal_traction.mul(stress, normal); Real normal_contrib = normal_traction.dot(normal); /// in 3D tangential components must be summed Real tangent_contrib = 0; if (dim == 2) { Real tangent_contrib_tmp = normal_traction.dot(tangent); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } else if (dim == 3) { for (UInt s = 0; s < dim - 1; ++s) { const Vector tangent_v(tangent.storage() + s * dim, dim); 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(Real(0.), normal_contrib); return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv); } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTractionOnQuad( Vector & traction, Real & delta_max, __attribute__((unused)) const Real & delta_max_prev, const Real & delta_c, const Vector & insertion_stress, const Real & sigma_c, Vector & opening, Vector & opening_prec, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, bool & reduction_penalty, Real & current_penalty, Vector & contact_traction, Vector & contact_opening) { /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); - normal_opening = normal; normal_opening *= normal_opening_norm; + tangential_opening = opening; tangential_opening -= normal_opening; - 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 * this->beta2_kappa2; - penetration = normal_opening_norm < -Math::getTolerance(); + penetration = normal_opening_norm < 0.0; if (this->contact_after_breaking == false && Math::are_float_equal(damage, 1.)) penetration = false; /** * if during the convergence loop a cohesive element continues to * jumps from penetration to opening, and convergence is not * reached, its penalty parameter will be reduced in the * recomputation of the same incremental step. recompute is set * equal to true when convergence is not reached in the * solveStepCohesive function and the execution of the program * goes back to the main file where the variable load_reduction * is set equal to true. */ Real normal_opening_prec_norm = opening_prec.dot(normal); opening_prec = opening; if (!this->model->isExplicit() && !this->recompute) - if ((normal_opening_prec_norm * normal_opening_norm) < 0) { + if ((normal_opening_prec_norm * normal_opening_norm) < 0.0) { reduction_penalty = true; } if (penetration) { if (this->recompute && reduction_penalty){ /// the penalty parameter is locally reduced current_penalty = this->penalty / 100.; } else current_penalty = this->penalty; /// use penalty coefficient in case of penetration contact_traction = normal_opening; contact_traction *= current_penalty; contact_opening = normal_opening; /// don't consider penetration contribution for delta opening = tangential_opening; normal_opening.clear(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction.clear(); contact_opening.clear(); } delta = std::sqrt(delta); /// update maximum displacement and damage delta_max = std::max(delta_max, delta); damage = std::min(delta_max / delta_c, Real(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, 1.)) traction.clear(); else if (Math::are_float_equal(damage, 0.)) { if (penetration) traction.clear(); else traction = insertion_stress; } else { traction = tangential_opening; traction *= this->beta2_kappa; traction += normal_opening; AKANTU_DEBUG_ASSERT(delta_max != 0., "Division by zero, tolerance might be too low"); traction *= sigma_c / delta_max * (1. - damage); } } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTangentTractionOnQuad( Matrix & tangent, Real & delta_max, const Real & delta_c, const Real & sigma_c, Vector & opening, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, bool & reduction_penalty, Real & current_penalty, Vector & contact_opening) { - /// During the update of the residual the interpenetrations are - /// stored in the array "contact_opening", therefore, in the case - /// of penetration, in the array "opening" there are only the - /// tangential components. + /** + * During the update of the residual the interpenetrations are + * stored in the array "contact_opening", therefore, in the case + * of penetration, in the array "opening" there are only the + * tangential components. + */ opening += contact_opening; /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; - tangential_opening_norm = tangential_opening.norm(); - penetration = normal_opening_norm < -Math::getTolerance(); + + Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; + + penetration = normal_opening_norm < 0.0; if (contact_after_breaking == false && Math::are_float_equal(damage, 1.)) penetration = false; Real derivative = 0; // derivative = d(t/delta)/ddelta Real t = 0; - Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; - Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); if (penetration){ - if (recompute && reduction_penalty) + if (this->recompute && reduction_penalty) current_penalty = this->penalty / 100.; else current_penalty = this->penalty; /// stiffness in compression given by the penalty parameter tangent += n_outer_n; tangent *= current_penalty; opening = tangential_opening; normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; } else{ delta += normal_opening_norm * normal_opening_norm; } delta = std::sqrt(delta); - /// Delta has to be different from 0 to have finite values of - /// tangential stiffness. At the element insertion, delta = - /// 0. Therefore, a fictictious value is defined, for the - /// evaluation of the first value of K. + /** + * Delta has to be different from 0 to have finite values of + * tangential stiffness. At the element insertion, delta = + * 0. Therefore, a fictictious value is defined, for the + * evaluation of the first value of K. + */ if (delta < Math::getTolerance()) delta = (delta_c)/1000.; if (delta >= delta_max){ if (delta <= delta_c){ derivative = -sigma_c / (delta * delta); t = sigma_c * (1 - delta / delta_c); } else { derivative = 0.; t = 0.; } } else if (delta < delta_max){ Real tmax = sigma_c * (1 - delta_max / delta_c); t = tmax / delta_max * delta; } /// computation of the derivative of the constitutive law (dT/ddelta) Matrix I(spatial_dimension, spatial_dimension); I.eye(this->beta2_kappa); Matrix nn(n_outer_n); nn *= (1. - this->beta2_kappa); nn += I; nn *= t/delta; Vector t_tilde(normal_opening); t_tilde *= (1. - this->beta2_kappa2); Vector mm(opening); mm *= this->beta2_kappa2; t_tilde += mm; Vector t_hat(normal_opening); t_hat += this->beta2_kappa * tangential_opening; Matrix prov(spatial_dimension, spatial_dimension); prov.outerProduct(t_hat, t_tilde); prov *= derivative/delta; prov += nn; +// Matrix tmp(spatial_dimension, spatial_dimension); +// for (UInt i = 0; i < spatial_dimension; ++i) { +// for (UInt j = 0; j < spatial_dimension; ++j) { +// tmp(j,i) = prov(i,j); +// } +// } + Matrix prov_t = prov.transpose(); tangent += prov_t; + // tangent += tmp; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #endif //__AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__