diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc index 8242f9948..b0a258e0b 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc @@ -1,204 +1,212 @@ /** * @file material_cohesive_bilinear.cc * * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Feb 20 2018 * * @brief Bilinear cohesive constitutive law * - * @section LICENSE * * Copyright (©) 2010-2018 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 . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_bilinear.hh" //#include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveBilinear::MaterialCohesiveBilinear( SolidMechanicsModel & model, const ID & id) : MaterialCohesiveLinear(model, id) { AKANTU_DEBUG_IN(); this->registerParam("delta_0", delta_0, Real(0.), _pat_parsable | _pat_readable, "Elastic limit displacement"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::initMaterial() { AKANTU_DEBUG_IN(); this->sigma_c_eff.setRandomDistribution(this->sigma_c.getRandomParameter()); MaterialCohesiveLinear::initMaterial(); this->delta_max.setDefaultValue(delta_0); this->insertion_stress.setDefaultValue(0); this->delta_max.reset(); this->insertion_stress.reset(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::onElementsAdded(element_list, event); bool scale_traction = false; // don't scale sigma_c if volume_s hasn't been specified by the user if (!Math::are_float_equal(this->volume_s, 0.)) scale_traction = true; - for (const auto & el : element_list) { + + auto el_it = element_list.begin(); + auto el_end = element_list.end(); + + for (auto & el: element_list) { // filter not ghost cohesive elements - if (el.ghost_type != _not_ghost || Mesh::getKind(el.type) != _ek_cohesive) + if (el.ghost_type != _not_ghost || + Mesh::getKind(el.type) != _ek_cohesive) continue; - UInt index = el.element; ElementType type = el.type; - UInt nb_element = this->element_filter(type).size(); - - UInt nb_quad_per_element = this->fem_cohesive.getNbIntegrationPoints(type); - - auto & sigma_c_eff = this->sigma_c_eff(type); + auto & elem_filter = this->element_filter(type, el.ghost_type); + auto filter_begin = elem_filter.begin(); + auto filter_end = elem_filter.end(); - auto sigma_c_begin = sigma_c_eff.begin_reinterpret( - nb_quad_per_element, nb_element); - auto delta_c_begin = this->delta_c_eff(type).begin_reinterpret( - nb_quad_per_element, nb_element); + UInt index = el.element; + if (std::find(filter_begin, filter_end, index) == filter_end) + continue; + UInt nb_quad_per_element = this->fem_cohesive.getNbIntegrationPoints(type); + auto sigma_c_begin = make_view(this->sigma_c_eff(type), nb_quad_per_element).begin(); Vector sigma_c_vec = sigma_c_begin[index]; + + auto delta_c_begin = make_view(this->delta_c_eff(type), nb_quad_per_element).begin(); Vector delta_c_vec = delta_c_begin[index]; if (scale_traction) - scaleTraction(el, sigma_c_vec); + scaleTraction(*el_it, sigma_c_vec); /** * Recompute sigma_c as * @f$ {\sigma_c}_\textup{new} = * \frac{{\sigma_c}_\textup{old} \delta_c} {\delta_c - \delta_0} @f$ */ for (UInt q = 0; q < nb_quad_per_element; ++q) { delta_c_vec(q) = 2 * this->G_c / sigma_c_vec(q); if (delta_c_vec(q) - delta_0 < Math::getTolerance()) AKANTU_ERROR("delta_0 = " << delta_0 << " must be lower than delta_c = " << delta_c_vec(q) << ", modify your material file"); sigma_c_vec(q) *= delta_c_vec(q) / (delta_c_vec(q) - delta_0); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::scaleTraction( const Element & el, Vector & sigma_c_vec) { AKANTU_DEBUG_IN(); Real base_sigma_c = this->sigma_c_eff; const Mesh & mesh_facets = this->model->getMeshFacets(); const FEEngine & fe_engine = this->model->getFEEngine(); auto coh_element_to_facet_begin = mesh_facets.getSubelementToElement(el.type).begin(2); const Vector & coh_element_to_facet = coh_element_to_facet_begin[el.element]; // compute bounding volume Real volume = 0; // loop over facets for (UInt f = 0; f < 2; ++f) { const Element & facet = coh_element_to_facet(f); - const auto & facet_to_element = + const Array> & facet_to_element = mesh_facets.getElementToSubelement(facet.type, facet.ghost_type); - const auto & element_list = facet_to_element(facet.element); + + const std::vector & element_list = facet_to_element(facet.element); + + auto elem = element_list.begin(); + auto elem_end = element_list.end(); // loop over elements connected to each facet - for (const auto & elem : element_list) { + for (; elem != elem_end; ++elem) { // skip cohesive elements and dummy elements - if (elem == ElementNull || Mesh::getKind(elem.type) == _ek_cohesive) + if (*elem == ElementNull || Mesh::getKind(elem->type) == _ek_cohesive) continue; // unit vector for integration in order to obtain the volume - UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem.type); + UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); - volume += fe_engine.integrate(unit_vector, elem.type, elem.element, - elem.ghost_type); + volume += fe_engine.integrate(unit_vector, elem->type, elem->element, + elem->ghost_type); } } // scale sigma_c sigma_c_vec -= base_sigma_c; sigma_c_vec *= std::pow(this->volume_s / volume, 1. / this->m_s); sigma_c_vec += base_sigma_c; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::computeTraction(normal, el_type, ghost_type); // adjust damage auto delta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); auto delta_max_it = this->delta_max(el_type, ghost_type).begin(); auto damage_it = this->damage(el_type, ghost_type).begin(); auto damage_end = this->damage(el_type, ghost_type).end(); for (; damage_it != damage_end; ++damage_it, ++delta_max_it, ++delta_c_it) { *damage_it = std::max((*delta_max_it - delta_0) / (*delta_c_it - delta_0), Real(0.)); *damage_it = std::min(*damage_it, Real(1.)); } } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_bilinear, MaterialCohesiveBilinear); } // namespace akantu