diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc index 7ac2e8459..ec2a35d07 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc @@ -1,714 +1,714 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jan 14 2016 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 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 . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinear::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), opening_prec("opening_prec", *this), reduction_penalty("reduction_penalty", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable | _pat_readable, "Beta parameter"); this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, "Mode I fracture energy"); this->registerParam("penalty", penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty coefficient"); this->registerParam("volume_s", volume_s, Real(0.), _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); this->registerParam("m_s", m_s, Real(1.), _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); this->registerParam("kappa", kappa, Real(1.), _pat_parsable | _pat_readable, "Kappa parameter"); this->registerParam("contact_after_breaking", contact_after_breaking, false, _pat_parsable | _pat_readable, "Activation of contact when the elements are fully damaged"); this->registerParam("max_quad_stress_insertion", max_quad_stress_insertion, false, _pat_parsable | _pat_readable, "Insertion of cohesive element when stress is high enough just on one quadrature point"); this->registerParam("recompute", recompute, false, _pat_parsable | _pat_modifiable, "recompute solution"); this->use_previous_delta_max = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); /// compute scalars - beta2_kappa2 = beta * beta/kappa/kappa; - beta2_kappa = beta * beta/kappa; + 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); opening_prec.initialize(spatial_dimension); reduction_penalty.initialize(1); 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(); opening_prec.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::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 > & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.getSize(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet Array::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_list = facet_to_element(f); // compute bounding volume Real volume = 0; std::vector::const_iterator elem = element_list.begin(); std::vector::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.getNbIntegrationPoints(elem->type); Vector 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 void MaterialCohesiveLinear::checkInsertion(bool check_only) { 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 & facets_check = inserter.getCheckFacets(type_facet); Array & f_insertion = inserter.getInsertionFacets(type_facet); Array & f_filter = facet_filter(type_facet); Array & sig_c_eff = sigma_c_eff(type_cohesive); Array & del_c = delta_c_eff(type_cohesive); Array & ins_stress = insertion_stress(type_cohesive); Array & trac_old = tractions_old(type_cohesive); Array & open_prec = opening_prec(type_cohesive); Array & red_penalty = reduction_penalty(type_cohesive); const Array & f_stress = model->getStressOnFacets(type_facet); const Array & sigma_lim = sigma_c(type_facet); Real max_ratio = 0.; UInt index_f = 0; UInt index_filter = 0; UInt nn = 0; UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); UInt nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine").getNbIntegrationPoints(type_cohesive); UInt nb_facet = f_filter.getSize(); // if (nb_facet == 0) continue; Array::const_iterator sigma_lim_it = sigma_lim.begin(); Matrix stress_tmp(spatial_dimension, spatial_dimension); Matrix normal_traction(spatial_dimension, nb_quad_facet); Vector stress_check(nb_quad_facet); UInt sp2 = spatial_dimension * spatial_dimension; const Array & tangents = model->getTangents(type_facet); const Array & normals = model->getFEEngine("FacetsFEEngine").getNormalsOnIntegrationPoints(type_facet); Array::const_vector_iterator normal_begin = normals.begin(spatial_dimension); Array::const_vector_iterator tangent_begin = tangents.begin(tangents.getNbComponent()); Array::const_matrix_iterator facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::vector< Vector > new_normal_traction; std::vector 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 = facet * nb_quad_facet + q; const Vector & normal = normal_begin[current_quad]; const Vector & tangent = tangent_begin[current_quad]; const Matrix & facet_stress_it = facet_stress_begin[current_quad]; // compute average stress on the current quadrature point Matrix stress_1(facet_stress_it.storage(), spatial_dimension, spatial_dimension); Matrix stress_2(facet_stress_it.storage() + sp2, spatial_dimension, spatial_dimension); stress_tmp.copy(stress_1); stress_tmp += stress_2; stress_tmp /= 2.; Vector 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 Real final_stress = stress_check.mean(); if (max_quad_stress_insertion) final_stress = *std::max_element(stress_check.storage(), stress_check.storage() + nb_quad_facet); if (final_stress > (*sigma_lim_it - tolerance)) { if (model->isExplicit()){ f_insertion(facet) = true; if (!check_only) { // 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 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); } } }else{ Real ratio = final_stress/(*sigma_lim_it); if (ratio > max_ratio){ ++nn; max_ratio = ratio; index_f = f; index_filter = f_filter(f); } } } } /// Insertion of only 1 cohesive element in case of implicit approach. The one subjected to the highest stress. if (!model->isExplicit()){ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Array abs_max(comm.getNbProc()); abs_max(comm.whoAmI()) = max_ratio; comm.allGather(abs_max.storage(), 1); Array::scalar_iterator it = std::max_element(abs_max.begin(), abs_max.end()); Int pos = it - abs_max.begin(); if (pos != comm.whoAmI()) { AKANTU_DEBUG_OUT(); return; } if (nn) { f_insertion(index_filter) = true; if (!check_only) { // Array::iterator > normal_traction_it = // normal_traction.begin_reinterpret(nb_quad_facet, spatial_dimension, nb_facet); Array::const_iterator sigma_lim_it = sigma_lim.begin(); for (UInt q = 0; q < nb_quad_cohesive; ++q) { Real new_sigma = (sigma_lim_it[index_f]); Vector normal_traction_vec(spatial_dimension, 0.0); 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 = delta_c; else new_delta = 2 * G_c / (new_sigma); 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); open_prec.resize(old_nb_quad_points + new_nb_quad_points); red_penalty.resize(old_nb_quad_points + new_nb_quad_points); for (UInt q = 0; q < new_nb_quad_points; ++q) { sig_c_eff(old_nb_quad_points + q) = new_sigmas[q]; del_c(old_nb_quad_points + q) = new_delta_c[q]; red_penalty(old_nb_quad_points + q) = false; for (UInt dim = 0; dim < spatial_dimension; ++dim) { ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); open_prec(old_nb_quad_points + q, dim) = 0.; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::vector_iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// opening_prec is the opening of the previous step in the /// Newton-Raphson loop Array::vector_iterator opening_prec_it = opening_prec(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator traction_end = tractions(el_type, ghost_type).end(spatial_dimension); Array::scalar_iterator sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator delta_max_it = delta_max(el_type, ghost_type).begin(); Array::scalar_iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); Array::vector_iterator insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); if (! this->model->isExplicit()) this->delta_max(el_type, ghost_type).copy(this->delta_max.previous(el_type, ghost_type)); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++opening_prec_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it, ++delta_max_prev_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; this->computeTractionOnQuad(*traction_it, *delta_max_it, *delta_max_prev_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, *opening_it, *opening_prec_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *reduction_penalty_it, current_penalty, *contact_traction_it, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkDeltaMax(GhostType ghost_type) { AKANTU_DEBUG_IN(); /** * This function set a predefined value to the parameter * delta_max_prev of the elements that have been inserted in the * last loading step for which convergence has not been * reached. This is done before reducing the loading and re-doing * the step. Otherwise, the updating of delta_max_prev would be * done with reference to the non-convergent solution. In this * function also other variables related to the contact and * friction behavior are correctly set. */ Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); /** * the variable "recompute" is set to true to activate the * procedure that reduce the penalty parameter for * compression. This procedure is set true only during the phase of * load_reduction, that has to be set in the maiin file. The * penalty parameter will be reduced only for the elements having * reduction_penalty = true. */ recompute = true; for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; /// define iterators Array::scalar_iterator delta_max_it = delta_max(el_type, ghost_type).begin(); Array::scalar_iterator delta_max_end = delta_max(el_type, ghost_type).end(); Array::scalar_iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); Array::vector_iterator opening_prec_it = opening_prec(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_prec_prev_it = opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); /// loop on each quadrature point for (; delta_max_it != delta_max_end; ++delta_max_it, ++delta_max_prev_it, ++delta_c_it, ++opening_prec_it, ++opening_prec_prev_it) { if (*delta_max_prev_it == 0) /// elements inserted in the last incremental step, that did /// not converge *delta_max_it = *delta_c_it / 1000; else /// elements introduced in previous incremental steps, for /// which a correct value of delta_max_prev already exists *delta_max_it = *delta_max_prev_it; /// in case convergence is not reached, set opening_prec to the /// value referred to the previous incremental step *opening_prec_it = *opening_prec_prev_it; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::resetVariables(GhostType ghost_type) { AKANTU_DEBUG_IN(); /** * This function set the variables "recompute" and * "reduction_penalty" to false. It is called by solveStepCohesive * when convergence is reached. Such variables, in fact, have to be * false at the beginning of a new incremental step. */ Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); recompute = false; for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Array::scalar_iterator reduction_penalty_end = reduction_penalty(el_type, ghost_type).end(); /// loop on each quadrature point for (; reduction_penalty_it != reduction_penalty_end; ++reduction_penalty_it) { *reduction_penalty_it = false; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::matrix_iterator tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step Array::scalar_iterator delta_max_it = delta_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); Array::vector_iterator contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; this->computeTangentTractionOnQuad(*tangent_it, *delta_max_it, *delta_c_it, *sigma_c_it, *opening_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *reduction_penalty_it, current_penalty, *contact_opening_it); // check if the tangential stiffness matrix is symmetric // for (UInt h = 0; h < spatial_dimension; ++h){ // for (UInt l = h; l < spatial_dimension; ++l){ // if (l > h){ // Real k_ls = (*tangent_it)[spatial_dimension*h+l]; // Real k_us = (*tangent_it)[spatial_dimension*l+h]; // // std::cout << "k_ls = " << k_ls << std::endl; // // std::cout << "k_us = " << k_us << std::endl; // if (std::abs(k_ls) > 1e-13 && std::abs(k_us) > 1e-13){ // Real error = std::abs((k_ls - k_us) / k_us); // if (error > 1e-10){ // std::cout << "non symmetric cohesive matrix" << std::endl; // // std::cout << "error " << error << std::endl; // } // } // } // } // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialCohesiveLinear); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc index 16ff1414b..df9d11baf 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc @@ -1,588 +1,590 @@ /** * @file material_cohesive_linear_uncoupled.cc * * @author Mauro Corrado * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jan 14 2016 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 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 . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear_uncoupled.hh" #include "solid_mechanics_model_cohesive.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearUncoupled ::MaterialCohesiveLinearUncoupled(SolidMechanicsModel & model, const ID & id) : MaterialCohesiveLinear(model,id), - delta_nt_max("delta_nt_max", *this), - delta_st_max("delta_st_max", *this), - damage_nt("damage_nt", *this), - damage_st("damage_st", *this){ + delta_n_max("delta_n_max", *this), + delta_t_max("delta_t_max", *this), + damage_n("damage_n", *this), + damage_t("damage_t", *this){ + AKANTU_DEBUG_IN(); this->registerParam("roughness", R, Real(1.), _pat_parsable | _pat_readable, "Roughness to define coupling between mode II and mode I"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearUncoupled::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::initMaterial(); - delta_nt_max.initialize(1); - delta_st_max.initialize(1); - damage_nt.initialize(1); - damage_st.initialize(1); + delta_n_max.initialize(1); + delta_t_max.initialize(1); + damage_n.initialize(1); + damage_t.initialize(1); - delta_nt_max.initializeHistory(); - delta_st_max.initializeHistory(); + delta_n_max.initializeHistory(); + delta_t_max.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearUncoupled::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - delta_nt_max.resize(); - delta_st_max.resize(); - damage_nt.resize(); - damage_st.resize(); + delta_n_max.resize(); + delta_t_max.resize(); + damage_n.resize(); + damage_t.resize(); /// define iterators Array::vector_iterator traction_it = this->tractions(el_type, ghost_type).begin(spatial_dimension); + Array::vector_iterator traction_end = + this->tractions(el_type, ghost_type).end(spatial_dimension); + Array::vector_iterator opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension); /// opening_prec is the opening of the previous step in the /// Newton-Raphson loop Array::vector_iterator opening_prec_it = this->opening_prec(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_traction_it = this->contact_tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_opening_it = this->contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::const_vector_iterator normal_it = this->normal.begin(spatial_dimension); - Array::vector_iterator traction_end = - this->tractions(el_type, ghost_type).end(spatial_dimension); - Array::scalar_iterator sigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin(); - Array::scalar_iterator delta_nt_max_it = - this->delta_nt_max(el_type, ghost_type).begin(); + Array::scalar_iterator delta_n_max_it = + this->delta_n_max(el_type, ghost_type).begin(); - Array::scalar_iterator delta_st_max_it = - this->delta_st_max(el_type, ghost_type).begin(); - - Array::scalar_iterator delta_nt_max_prev_it = - this->delta_nt_max.previous(el_type, ghost_type).begin(); - - Array::scalar_iterator delta_st_max_prev_it = - this->delta_st_max.previous(el_type, ghost_type).begin(); + Array::scalar_iterator delta_t_max_it = + this->delta_t_max(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); - Array::scalar_iterator damage_nt_it = - this->damage_nt(el_type, ghost_type).begin(); + Array::scalar_iterator damage_n_it = + this->damage_n(el_type, ghost_type).begin(); - Array::scalar_iterator damage_st_it = - this->damage_st(el_type, ghost_type).begin(); + Array::scalar_iterator damage_t_it = + this->damage_t(el_type, ghost_type).begin(); Array::vector_iterator insertion_stress_it = this->insertion_stress(el_type, ghost_type).begin(spatial_dimension); Array::scalar_iterator reduction_penalty_it = this->reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); - if (! this->model->isExplicit()){ - this->delta_nt_max(el_type, ghost_type).copy(this->delta_nt_max.previous(el_type, ghost_type)); - this->delta_st_max(el_type, ghost_type).copy(this->delta_st_max.previous(el_type, ghost_type)); + if (!this->model->isExplicit()){ + this->delta_n_max(el_type, ghost_type).copy(this->delta_n_max.previous(el_type, ghost_type)); + this->delta_t_max(el_type, ghost_type).copy(this->delta_t_max.previous(el_type, ghost_type)); } /// loop on each quadrature point for (; traction_it != traction_end; - ++traction_it, ++opening_it, ++opening_prec_it, - ++normal_it, ++sigma_c_it, ++delta_c_it, ++delta_nt_max_it, - ++delta_st_max_it, ++damage_nt_it, ++damage_st_it, - ++contact_traction_it, ++contact_opening_it, ++insertion_stress_it, - ++reduction_penalty_it) { + ++traction_it, ++opening_it, ++opening_prec_it, + ++contact_traction_it, ++contact_opening_it, ++normal_it, + ++sigma_c_it, ++delta_n_max_it, ++delta_t_max_it, + ++delta_c_it, ++damage_n_it, ++damage_t_it, + ++insertion_stress_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; - Real delta_c_s = this->kappa / this->beta * (*delta_c_it); - Real delta_c_s2_R2 = delta_c_s * delta_c_s / (R * R); + Real delta_c2_R2 = *delta_c_it * (*delta_c_it) / (R * R); /// compute normal and tangential opening vectors normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening = *normal_it; normal_opening *= normal_opening_norm; Vector tangential_opening = *opening_it; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); /// compute effective opening displacement - Real delta_nt = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; - Real delta_st = tangential_opening_norm * tangential_opening_norm; + Real delta_n = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; + Real delta_t = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; 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 + * 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 + * 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_it->dot(*normal_it); - *opening_prec_it = *opening_it; - - if (!this->model->isExplicit() && !this->recompute) - if ((normal_opening_prec_norm * normal_opening_norm) < 0.0) { + Vector normal_opening_prec = *normal_it; + normal_opening_prec *= normal_opening_prec_norm; + if (!this->model->isExplicit()) // && !this->recompute) + if ((normal_opening_prec_norm * normal_opening_norm) < 0.0) *reduction_penalty_it = true; - } + *opening_prec_it = *opening_it; + if (penetration) { if (this->recompute && *reduction_penalty_it){ /// the penalty parameter is locally reduced - current_penalty = this->penalty / 100.; + current_penalty = this->penalty / 10.; } else current_penalty = this->penalty; /// use penalty coefficient in case of penetration *contact_traction_it = normal_opening; *contact_traction_it *= current_penalty; *contact_opening_it = normal_opening; /// don't consider penetration contribution for delta *opening_it = tangential_opening; normal_opening.clear(); } else { - delta_nt += normal_opening_norm * normal_opening_norm; - delta_st += normal_opening_norm * normal_opening_norm * delta_c_s2_R2; + delta_n += normal_opening_norm * normal_opening_norm; + delta_t += normal_opening_norm * normal_opening_norm * delta_c2_R2; contact_traction_it->clear(); contact_opening_it->clear(); } - delta_nt = std::sqrt(delta_nt); - delta_st = std::sqrt(delta_st); + delta_n = std::sqrt(delta_n); + delta_t = std::sqrt(delta_t); /// update maximum displacement and damage - *delta_nt_max_it = std::max(*delta_nt_max_it, delta_nt); - *damage_nt_it = std::min(*delta_nt_max_it / *delta_c_it, Real(1.)); + *delta_n_max_it = std::max(*delta_n_max_it, delta_n); + *damage_n_it = std::min(*delta_n_max_it / *delta_c_it, Real(1.)); - *delta_st_max_it = std::max(*delta_st_max_it, delta_st); - *damage_st_it = std::min(*delta_st_max_it / delta_c_s, Real(1.)); + *delta_t_max_it = std::max(*delta_t_max_it, delta_t); + *damage_t_it = std::min(*delta_t_max_it / *delta_c_it, 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$ - */ - - Vector traction_n(spatial_dimension); - Vector traction_s(spatial_dimension); + Vector normal_traction(spatial_dimension); + Vector shear_traction(spatial_dimension); /// NORMAL TRACTIONS - if (Math::are_float_equal(*damage_nt_it, 1.)) - traction_n.clear(); - else if (Math::are_float_equal(*damage_nt_it, 0.)) { + if (Math::are_float_equal(*damage_n_it, 1.)) + normal_traction.clear(); + else if (Math::are_float_equal(*damage_n_it, 0.)) { if (penetration) - traction_n.clear(); + normal_traction.clear(); else - traction_n = *insertion_stress_it; + normal_traction = *insertion_stress_it; } else { - traction_n = normal_opening; + // the following formulation holds both in loading and in + // unloading-reloading + normal_traction = normal_opening; - AKANTU_DEBUG_ASSERT(*delta_nt_max_it != 0., + AKANTU_DEBUG_ASSERT(*delta_n_max_it != 0., "Division by zero, tolerance might be too low"); - traction_n *= *sigma_c_it / (*delta_nt_max_it) * (1. - *damage_nt_it); + normal_traction *= *sigma_c_it / (*delta_n_max_it) * (1. - *damage_n_it); } /// SHEAR TRACTIONS - if (Math::are_float_equal(*damage_st_it, 1.)) - traction_s.clear(); - else if (Math::are_float_equal(*damage_st_it, 0.)) { - traction_s.clear(); + if (Math::are_float_equal(*damage_t_it, 1.)) + shear_traction.clear(); + else if (Math::are_float_equal(*damage_t_it, 0.)) { + shear_traction.clear(); } else { - traction_s = tangential_opening; + shear_traction = tangential_opening; - AKANTU_DEBUG_ASSERT(*delta_st_max_it != 0., + AKANTU_DEBUG_ASSERT(*delta_t_max_it != 0., "Division by zero, tolerance might be too low"); - traction_s *= *sigma_c_it * this->beta / (*delta_st_max_it) * (1. - *damage_st_it); + shear_traction *= *sigma_c_it * this->beta2_kappa / (*delta_t_max_it) * (1. - *damage_t_it); } - *traction_it = traction_n + traction_s; + *traction_it = normal_traction; + *traction_it += shear_traction; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearUncoupled::checkDeltaMax(GhostType ghost_type) { AKANTU_DEBUG_IN(); /** * This function set a predefined value to the parameter * delta_max_prev of the elements that have been inserted in the * last loading step for which convergence has not been * reached. This is done before reducing the loading and re-doing * the step. Otherwise, the updating of delta_max_prev would be * done with reference to the non-convergent solution. In this * function also other variables related to the contact and * friction behavior are correctly set. */ Mesh & mesh = this->fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); /** * the variable "recompute" is set to true to activate the - * procedure that reduce the penalty parameter for - * compression. This procedure is set true only during the phase of - * load_reduction, that has to be set in the maiin file. The + * procedure that reduces the penalty parameter for + * compression. This procedure is available only during the phase of + * load_reduction, that has to be set in the main file. The * penalty parameter will be reduced only for the elements having * reduction_penalty = true. */ this->recompute = true; for(; it != last_type; ++it) { Array & elem_filter = this->element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; + // std::cout << "element type: " << el_type << std::endl; + /// define iterators - Array::scalar_iterator delta_nt_max_it = - delta_nt_max(el_type, ghost_type).begin(); + Array::scalar_iterator delta_n_max_it = + delta_n_max(el_type, ghost_type).begin(); - Array::scalar_iterator delta_nt_max_end = - delta_nt_max(el_type, ghost_type).end(); + Array::scalar_iterator delta_n_max_end = + delta_n_max(el_type, ghost_type).end(); - Array::scalar_iterator delta_nt_max_prev_it = - delta_nt_max.previous(el_type, ghost_type).begin(); + Array::scalar_iterator delta_n_max_prev_it = + delta_n_max.previous(el_type, ghost_type).begin(); - Array::scalar_iterator delta_st_max_it = - delta_st_max(el_type, ghost_type).begin(); + Array::scalar_iterator delta_t_max_it = + delta_t_max(el_type, ghost_type).begin(); - Array::scalar_iterator delta_st_max_prev_it = - delta_st_max.previous(el_type, ghost_type).begin(); + Array::scalar_iterator delta_t_max_prev_it = + delta_t_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = this-> delta_c_eff(el_type, ghost_type).begin(); Array::vector_iterator opening_prec_it = this->opening_prec(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_prec_prev_it = this->opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); + Int it = 0; + /// loop on each quadrature point - for (; delta_nt_max_it != delta_nt_max_end; - ++delta_nt_max_it, ++delta_st_max_it, ++delta_c_it, - ++delta_nt_max_prev_it, ++delta_st_max_prev_it, + for (; delta_n_max_it != delta_n_max_end; + ++delta_n_max_it, ++delta_t_max_it, ++delta_c_it, + ++delta_n_max_prev_it, ++delta_t_max_prev_it, ++opening_prec_it, ++opening_prec_prev_it) { - if (*delta_nt_max_prev_it == 0) + ++it; + + if (*delta_n_max_prev_it == 0){ /// elements inserted in the last incremental step, that did /// not converge - *delta_nt_max_it = *delta_c_it / 1000; + *delta_n_max_it = *delta_c_it / 1000.; + } else /// elements introduced in previous incremental steps, for /// which a correct value of delta_max_prev already exists - *delta_nt_max_it = *delta_nt_max_prev_it; + *delta_n_max_it = *delta_n_max_prev_it; - if (*delta_st_max_prev_it == 0) - *delta_st_max_it = *delta_c_it * this->kappa / this->beta / 1000; + if (*delta_t_max_prev_it == 0){ + *delta_t_max_it = *delta_c_it * this->kappa / this->beta / 1000.; + } else - *delta_st_max_it = *delta_st_max_prev_it; + *delta_t_max_it = *delta_t_max_prev_it; /// in case convergence is not reached, set opening_prec to the /// value referred to the previous incremental step *opening_prec_it = *opening_prec_prev_it; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearUncoupled::computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::matrix_iterator tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step - Array::scalar_iterator delta_nt_max_it = - this->delta_nt_max.previous(el_type, ghost_type).begin(); + Array::scalar_iterator delta_n_max_it = + this->delta_n_max.previous(el_type, ghost_type).begin(); - Array::scalar_iterator delta_st_max_it = - this->delta_st_max.previous(el_type, ghost_type).begin(); + Array::scalar_iterator delta_t_max_it = + this->delta_t_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator sigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); - // Array::scalar_iterator damage_it = - // damage(el_type, ghost_type).begin(); - Array::vector_iterator contact_opening_it = this->contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::scalar_iterator reduction_penalty_it = this->reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++sigma_c_it, - ++delta_c_it, ++delta_nt_max_it, ++delta_st_max_it, + ++delta_c_it, ++delta_n_max_it, ++delta_t_max_it, ++contact_opening_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; - Real delta_c_s = this->kappa / this->beta * (*delta_c_it); - Real delta_c_s2_R2 = delta_c_s * delta_c_s / (R * R); + Real delta_c2_R2 = *delta_c_it * (*delta_c_it) / (R * R); /** * 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_it += *contact_opening_it; /// compute normal and tangential opening vectors normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening = *normal_it; normal_opening *= normal_opening_norm; Vector tangential_opening = *opening_it; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); - Real delta_nt = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; - Real delta_st = tangential_opening_norm * tangential_opening_norm; + Real delta_n = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; + Real delta_t = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm < 0.0; // if (this->contact_after_breaking == false && Math::are_float_equal(damage, 1.)) // penetration = false; Real derivative = 0; // derivative = d(t/delta)/ddelta - Real t = 0; - // Matrix tangent_n(spatial_dimension, spatial_dimension); - // Matrix tangent_s(spatial_dimension, spatial_dimension); + Real T = 0; /// TANGENT STIFFNESS FOR NORMAL TRACTIONS Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(*normal_it, *normal_it); if (penetration){ if (this->recompute && *reduction_penalty_it) - current_penalty = this->penalty / 100.; + current_penalty = this->penalty / 10.; else current_penalty = this->penalty; /// stiffness in compression given by the penalty parameter *tangent_it = n_outer_n; *tangent_it *= current_penalty; *opening_it = tangential_opening; normal_opening.clear(); } else{ - delta_nt += normal_opening_norm * normal_opening_norm; - delta_nt = std::sqrt(delta_nt); + delta_n += normal_opening_norm * normal_opening_norm; + delta_n = std::sqrt(delta_n); - delta_st += normal_opening_norm * normal_opening_norm * delta_c_s2_R2; + delta_t += normal_opening_norm * normal_opening_norm * delta_c2_R2; /** * 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_nt < Math::getTolerance()) - delta_nt = *delta_c_it / 1000.; - - if (delta_nt >= *delta_nt_max_it){ - if (delta_nt <= *delta_c_it){ - derivative = - (*sigma_c_it) / (delta_nt * delta_nt); - t = *sigma_c_it * (1 - delta_nt / *delta_c_it); + if (delta_n < Math::getTolerance()) + delta_n = *delta_c_it / 1000.; + + // loading + if (delta_n >= *delta_n_max_it){ + if (delta_n <= *delta_c_it){ + derivative = - (*sigma_c_it) / (delta_n * delta_n); + T = *sigma_c_it * (1 - delta_n / *delta_c_it); } else { derivative = 0.; - t = 0.; + T = 0.; } - } else if (delta_nt < *delta_nt_max_it){ - Real tmax = *sigma_c_it * (1 - *delta_nt_max_it / *delta_c_it); - t = tmax / *delta_nt_max_it * delta_nt; + // unloading-reloading + } else if (delta_n < *delta_n_max_it){ + Real T_max = *sigma_c_it * (1 - *delta_n_max_it / *delta_c_it); + derivative = 0.; + T = T_max / *delta_n_max_it * delta_n; } /// computation of the derivative of the constitutive law (dT/ddelta) Matrix nn(n_outer_n); - nn *= t / delta_nt; + nn *= T / delta_n; - Vector t_tilde(normal_opening); - t_tilde *= (1. - this->beta2_kappa2); + Vector Delta_tilde(normal_opening); + Delta_tilde *= (1. - this->beta2_kappa2); Vector mm(*opening_it); mm *= this->beta2_kappa2; - t_tilde += mm; + Delta_tilde += mm; - Vector t_hat(normal_opening); + Vector Delta_hat(normal_opening); Matrix prov(spatial_dimension, spatial_dimension); - prov.outerProduct(t_hat, t_tilde); - prov *= derivative / delta_nt; + prov.outerProduct(Delta_hat, Delta_tilde); + prov *= derivative / delta_n; prov += nn; *tangent_it = prov.transpose(); } derivative = 0.; - t = 0.; + T = 0.; + /// TANGENT STIFFNESS FOR SHEAR TRACTIONS - delta_st = std::sqrt(delta_st); + delta_t = std::sqrt(delta_t); /** * 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_st < Math::getTolerance()) - delta_st = delta_c_s / 1000.; - - if (delta_st >= *delta_st_max_it){ - if (delta_st <= delta_c_s){ - derivative = - this->beta * (*sigma_c_it) / (delta_st * delta_st); - t = this->beta * (*sigma_c_it) * (1 - delta_st / delta_c_s); + if (delta_t < Math::getTolerance()) + delta_t = *delta_c_it / 1000.; + + // loading + if (delta_t >= *delta_t_max_it){ + if (delta_t <= *delta_c_it){ + derivative = - (*sigma_c_it) / (delta_t * delta_t); + T = *sigma_c_it * (1 - delta_t / *delta_c_it); } else { derivative = 0.; - t = 0.; + T = 0.; } - } else if (delta_st < *delta_st_max_it){ - Real tmax = this->beta * (*sigma_c_it) * (1 - *delta_st_max_it / delta_c_s); - t = tmax / *delta_st_max_it * delta_st; + // unloading-reloading + } else if (delta_t < *delta_t_max_it){ + Real T_max = (*sigma_c_it) * (1 - *delta_t_max_it / *delta_c_it); + derivative = 0.; + T = T_max / *delta_t_max_it * delta_t; } /// computation of the derivative of the constitutive law (dT/ddelta) Matrix I(spatial_dimension, spatial_dimension); I.eye(); Matrix nn(n_outer_n); I -= nn; - I *= t / delta_st; + I *= T / delta_t; - Vector t_tilde(normal_opening); - t_tilde *= (delta_c_s2_R2 - 1.); + Vector Delta_tilde(normal_opening); + Delta_tilde *= (delta_c2_R2 - this->beta2_kappa2); Vector mm(*opening_it); - t_tilde += mm; + mm *= this->beta2_kappa2; + Delta_tilde += mm; - Vector t_hat(tangential_opening); + Vector Delta_hat(tangential_opening); + Delta_hat *= this->beta2_kappa; Matrix prov(spatial_dimension, spatial_dimension); - prov.outerProduct(t_hat, t_tilde); - prov *= derivative / delta_st; + prov.outerProduct(Delta_hat, Delta_tilde); + prov *= derivative / delta_t; prov += I; - Matrix prov_t = prov.transpose(); - *tangent_it += prov_t; + // Matrix prov_t = prov.transpose(); + *tangent_it += prov.transpose(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialCohesiveLinearUncoupled); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh index b81d0f7e6..393d151ab 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh @@ -1,115 +1,115 @@ /** * @file material_cohesive_linear.hh * * @author Mauro Corrado * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jan 14 2016 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 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 . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_UNCOUPLED_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_UNCOUPLED_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Cohesive material linear with two different laws for mode I and * mode II, for extrinsic case * * parameters in the material files : * - roughness : define the interaction between mode I and mode II (default: 0) */ template class MaterialCohesiveLinearUncoupled : public MaterialCohesiveLinear { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ // typedef MaterialCohesiveLinear MaterialParent; public: MaterialCohesiveLinearUncoupled(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material parameters virtual void initMaterial(); protected: /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost); /// check delta_max for cohesive elements in case of no convergence /// in the solveStep (only for extrinsic-implicit) virtual void checkDeltaMax(GhostType ghost_type = _not_ghost); /// compute tangent stiffness matrix virtual void computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// parameter to tune the interaction between mode II and mode I Real R; /// maximum normal displacement - CohesiveInternalField delta_nt_max; + CohesiveInternalField delta_n_max; /// maximum tangential displacement - CohesiveInternalField delta_st_max; + CohesiveInternalField delta_t_max; /// damage associated to normal tractions - CohesiveInternalField damage_nt; + CohesiveInternalField damage_n; /// damage associated to shear tractions - CohesiveInternalField damage_st; + CohesiveInternalField damage_t; }; __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_UNCOUPLED_HH__ */