diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc index b557e9222..163017817 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc @@ -1,432 +1,440 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Wed Feb 21 2018 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * * 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 "dof_synchronizer.hh" #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace 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) { 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(); sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); if (not 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 void MaterialCohesiveLinear::updateInternalParameters() { /// 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; } } /* -------------------------------------------------------------------------- */ 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 auto & fe_engine = model->getFEEngine(); const auto & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); Real base_sigma_c = sigma_c; for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { const Array> & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.size(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet auto 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; auto elem = element_list.begin(); auto 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(); for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const auto & facets_check = inserter.getCheckFacets(type_facet); auto & f_insertion = inserter.getInsertionFacets(type_facet); auto & f_filter = facet_filter(type_facet); auto & sig_c_eff = sigma_c_eff(type_cohesive); auto & del_c = delta_c_eff(type_cohesive); auto & ins_stress = insertion_stress(type_cohesive); auto & trac_old = tractions.previous(type_cohesive); const auto & f_stress = model->getStressOnFacets(type_facet); const auto & sigma_lim = sigma_c(type_facet); UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); #ifndef AKANTU_NDEBUG UInt nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive); AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, "The cohesive element and the corresponding facet do " "not have the same numbers of integration points"); #endif UInt nb_facet = f_filter.size(); // if (nb_facet == 0) continue; auto 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 auto & tangents = model->getTangents(type_facet); const auto & normals = model->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(type_facet); auto normal_begin = normals.begin(spatial_dimension); auto tangent_begin = tangents.begin(tangents.getNbComponent()); auto facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::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) { f_insertion(facet) = true; if (check_only) { continue; } // 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); } } } // update material data for the new elements UInt old_nb_quad_points = sig_c_eff.size(); 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); 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]; 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); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); auto contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); auto normal_opening_norm_it = normal_opening_norm(el_type, ghost_type).begin(); + auto tangential_opening_norm_it = + tangential_opening_norm(el_type, ghost_type).begin(); auto normal_it = normal.begin(spatial_dimension); auto traction_end = tractions(el_type, ghost_type).end(spatial_dimension); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_max_it = delta_max(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(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, ++normal_opening_norm_it) { - Real tangential_opening_norm{0}; + ++contact_opening_it, ++normal_opening_norm_it, + ++tangential_opening_norm_it) { bool penetration{false}; this->computeTractionOnQuad( *traction_it, *opening_it, *normal_it, *delta_max_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, normal_opening, tangential_opening, - *normal_opening_norm_it, tangential_opening_norm, *damage_it, + *normal_opening_norm_it, *tangential_opening_norm_it, *damage_it, penetration, *contact_traction_it, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTangentTraction( ElementType el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); + auto normal_opening_norm_it = + normal_opening_norm(el_type, ghost_type).begin(); + auto tangential_opening_norm_it = + tangential_opening_norm(el_type, ghost_type).begin(); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step auto delta_max_it = delta_max.previous(el_type, ghost_type).begin(); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); 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) { - Real normal_opening_norm{0}; - Real tangential_opening_norm{0}; + 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, + ++normal_opening_norm_it, ++tangential_opening_norm_it) { + // Real normal_opening_norm{0}; + // Real tangential_opening_norm{0}; bool penetration{false}; 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, *contact_opening_it); + *normal_it, normal_opening, tangential_opening, *normal_opening_norm_it, + *tangential_opening_norm_it, *damage_it, penetration, + *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear, MaterialCohesiveLinear); } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh index 9d9ac04bc..ddfccea14 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh @@ -1,269 +1,267 @@ /** * @file material_cohesive_linear_inline_impl.hh * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Apr 22 2015 * @date last modification: Wed Feb 21 2018 * * @brief Inline functions of the MaterialCohesiveLinear * * * Copyright (©) 2015-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_linear.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH_ #define AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH_ /* -------------------------------------------------------------------------- */ namespace 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, Vector & opening, const Vector & normal, Real & delta_max, const Real & delta_c, const Vector & insertion_stress, const Real & sigma_c, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, 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 / delta_c < -Math::getTolerance(); // penetration = normal_opening_norm < 0.; - if (not this->contact_after_breaking and - Math::are_float_equal(damage, 1.)) { + if (not this->contact_after_breaking and Math::are_float_equal(damage, 1.)) { penetration = false; } if (penetration) { /// use penalty coefficient in case of penetration contact_traction = normal_opening; contact_traction *= this->penalty; contact_opening = normal_opening; /// don't consider penetration contribution for delta opening = tangential_opening; normal_opening.zero(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction.zero(); contact_opening.zero(); } 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.zero(); } else if (Math::are_float_equal(damage, 0.)) { if (penetration) { traction.zero(); } 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, 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. */ 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(); Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm < 0.0; - if (not this->contact_after_breaking and - Math::are_float_equal(damage, 1.)) { + if (not this->contact_after_breaking and Math::are_float_equal(damage, 1.)) { penetration = false; } Real derivative = 0; // derivative = d(t/delta)/ddelta Real t = 0; Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); if (penetration) { /// stiffness in compression given by the penalty parameter tangent += n_outer_n; tangent *= 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. */ 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 prov_t = prov.transpose(); tangent += prov_t; } /* -------------------------------------------------------------------------- */ } // namespace akantu /* -------------------------------------------------------------------------- */ -#endif //AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH_ +#endif // AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH_ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc index 1335f665b..7f22d48b4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc @@ -1,574 +1,578 @@ /** * @file material_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Mon Feb 19 2018 * * @brief Specialization of the material class for cohesive elements * * * 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.hh" #include "aka_random_generator.hh" #include "dof_synchronizer.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id, this->getMemoryID()), fem_cohesive( model.getFEEngineClass("CohesiveFEEngine")), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), eigen_opening("eigen_opening", *this), tractions("tractions", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), normal_opening_norm("normal_opening_norm", *this), + tangential_opening_norm("tangential_opening_norm", *this), delta_max("delta max", *this), use_previous_delta_max(false), use_previous_opening(false), damage("damage", *this), sigma_c("sigma_c", *this), normals("normals", *this) { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, Real(0.), _pat_parsable | _pat_readable, "Critical displacement"); this->element_filter.initialize(this->model->getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_cohesive); // this->model->getMesh().initElementTypeMapArray( // this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) { this->facet_filter.initialize(this->model->getMeshFacets(), _spatial_dimension = spatial_dimension - 1, _element_kind = _ek_regular); } // this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, // spatial_dimension - // 1); this->reversible_energy.initialize(1); this->total_energy.initialize(1); this->tractions.initialize(spatial_dimension); this->tractions.initializeHistory(); this->contact_tractions.initialize(spatial_dimension); this->contact_opening.initialize(spatial_dimension); this->normal_opening_norm.initialize(1); - //this->normal_opening_norm.initializeHistory(); + // this->normal_opening_norm.initializeHistory(); + + this->tangential_opening_norm.initialize(1); this->opening.initialize(spatial_dimension); - //this->opening.initializeHistory(); + this->opening.initializeHistory(); this->normals.initialize(spatial_dimension); this->eigen_opening.initialize(spatial_dimension); this->eigen_opening.setDefaultValue(0.); this->delta_max.initialize(1); this->damage.initialize(1); if (this->model->getIsExtrinsic()) { this->sigma_c.initialize(1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() = default; /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (this->use_previous_delta_max) { this->delta_max.initializeHistory(); } if (this->use_previous_opening) { this->opening.initializeHistory(); this->normal_opening_norm.initializeHistory(); + this->tangential_opening_norm.initializeHistory(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Tractions", tractions); #endif auto & internal_force = const_cast &>(model->getInternalForce()); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } const auto & shapes = fem_cohesive.getShapes(type, ghost_type); auto & traction = tractions(type, ghost_type); auto & contact_traction = contact_tractions(type, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); /// compute @f$t_i N_a@f$ auto * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); auto traction_it = traction.begin(spatial_dimension, 1); auto contact_traction_it = contact_traction.begin(spatial_dimension, 1); auto shapes_filtered_begin = shapes.begin(1, size_of_shapes); auto traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ auto * partial_int_t_N = new Array( nb_element, spatial_dimension * size_of_shapes, "int_t_N"); fem_cohesive.integrate(*traction_cpy, *partial_int_t_N, spatial_dimension * size_of_shapes, type, ghost_type, elem_filter); delete traction_cpy; auto * int_t_N = new Array( nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N"); Real * int_t_N_val = int_t_N->storage(); Real * partial_int_t_N_val = partial_int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val); std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val + size_of_shapes * spatial_dimension); for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) { int_t_N_val[n] *= -1.; } int_t_N_val += nb_nodes_per_element * spatial_dimension; partial_int_t_N_val += size_of_shapes * spatial_dimension; } delete partial_int_t_N; /// assemble model->getDOFManager().assembleElementalArrayLocalArray( *int_t_N, internal_force, type, ghost_type, 1, elem_filter); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & shapes = fem_cohesive.getShapes(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0U) { continue; } UInt size_of_shapes = shapes.getNbComponent(); auto * shapes_filtered = new Array(nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { auto * shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } Matrix A(spatial_dimension * size_of_shapes, spatial_dimension * nb_nodes_per_element); for (UInt i = 0; i < spatial_dimension * size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension * size_of_shapes) = -1; } /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} /// @f$ auto * tangent_stiffness_matrix = new Array( nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ // computeOpening(model->getDisplacement(), opening(type, ghost_type), type, // ghost_type); tangent_stiffness_matrix->zero(); computeTangentTraction(type, *tangent_stiffness_matrix, normals(type, ghost_type), ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension * nb_nodes_per_element * spatial_dimension * nb_nodes_per_element; auto * at_nt_d_n_a = new Array(nb_element * nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator> shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N(spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for (; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.zero(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt n = 0; n < size_of_shapes; ++n) { N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); } } /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} *{\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; auto * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive.integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, type, ghost_type, elem_filter); delete at_nt_d_n_a; model->getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _unsymmetric, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Openings", opening); #endif for (const auto & type : element_filter.elementTypes( spatial_dimension, ghost_type, _ek_cohesive)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(type, ghost_type), type, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normals(type, ghost_type), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); normal.zero(); #define COMPUTE_NORMAL(type) \ fem_cohesive.getShapeFunctions() \ .computeNormalsOnIntegrationPoints( \ position, normal, ghost_type, element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); #define COMPUTE_OPENING(type) \ fem_cohesive.getShapeFunctions() \ .interpolateOnIntegrationPoints( \ displacement, opening, spatial_dimension, ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING opening -= eigen_opening(type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::updateEnergies(ElementType type) { AKANTU_DEBUG_IN(); if (Mesh::getKind(type) != _ek_cohesive) { return; } Vector b(spatial_dimension); Vector h(spatial_dimension); auto erev = reversible_energy(type).begin(); auto etot = total_energy(type).begin(); auto traction_it = tractions(type).begin(spatial_dimension); auto traction_old_it = tractions.previous(type).begin(spatial_dimension); auto opening_it = opening(type).begin(spatial_dimension); auto opening_old_it = opening.previous(type).begin(spatial_dimension); auto traction_end = tractions(type).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } /// update old values AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { erev += fem_cohesive.integrate(reversible_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { Array dissipated_energy(total_energy(type, _not_ghost)); dissipated_energy -= reversible_energy(type, _not_ghost); edis += fem_cohesive.integrate(dissipated_energy, type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { auto & el_filter = element_filter(type, _not_ghost); UInt nb_quad_per_el = fem_cohesive.getNbIntegrationPoints(type, _not_ghost); UInt nb_quad_points = el_filter.size() * nb_quad_per_el; Array contact_energy(nb_quad_points); auto contact_traction_it = contact_tractions(type, _not_ghost).begin(spatial_dimension); auto contact_opening_it = contact_opening(type, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt q = 0; q < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++q) { contact_energy(q) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive.integrate(contact_energy, type, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(const std::string & type) { if (type == "reversible") { return getReversibleEnergy(); } if (type == "dissipated") { return getDissipatedEnergy(); } if (type == "cohesive contact") { return getContactEnergy(); } return 0.; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh index c65a03496..670913275 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh @@ -1,257 +1,259 @@ /** * @file material_cohesive.hh * * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Feb 21 2018 * * @brief Specialization of the material class for cohesive elements * * * 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.hh" /* -------------------------------------------------------------------------- */ #include "cohesive_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_COHESIVE_HH_ #define AKANTU_MATERIAL_COHESIVE_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } namespace akantu { class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineCohesiveType = FEEngineTemplate; public: MaterialCohesive(SolidMechanicsModel & model, const ID & id = ""); ~MaterialCohesive() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter void initMaterial() override; /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleInternalForces(GhostType ghost_type = _not_ghost) override; /// check stress for cohesive elements' insertion, by default it /// also updates the cohesive elements' data virtual void checkInsertion(bool /*check_only*/ = false) { AKANTU_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(const ElementType /*type*/, Array & /*result*/) {} /// compute the stresses void computeAllStresses(GhostType /*ghost_type*/ = _not_ghost) override{}; // add the facet to be handled by the material UInt addFacet(const Element & element); protected: virtual void computeTangentTraction(ElementType /*el_type*/, Array & /*tangent_matrix*/, const Array & /*normal*/, GhostType /*ghost_type*/ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } /// compute the normal void computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type); /// compute the opening void computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type); template void computeNormal(const Array & position, Array & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type) override; /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; protected: void updateEnergies(ElementType el_type) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the eigen opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE(EigenOpening, eigen_opening, Real); /// get the contact opening - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ContactOpening, contact_opening, - Real); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ContactOpening, contact_opening, Real); /// get the normal opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NormalOpeningNorm, normal_opening_norm, - Real); + Real); /// get the eigen opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(DeltaMax, delta_max, Real); /// get the normals AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Normals, normals, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const // ElementTypeMapArray &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy Real getEnergy(const std::string & type) override; /// return the energy (identified by id) for the provided element Real getEnergy(const std::string & energy_id, ElementType type, UInt index) override { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray facet_filter; /// Link to the cohesive fem object in the model FEEngine & fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField reversible_energy; /// total energy by quadrature point CohesiveInternalField total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField opening; /// eigen opening (fix value removed from the computed opening) CohesiveInternalField eigen_opening; /// traction in all elements and quadrature points CohesiveInternalField tractions; /// traction due to contact CohesiveInternalField contact_tractions; /// normal openings for contact tractions CohesiveInternalField contact_opening; /// normal openings CohesiveInternalField normal_opening_norm; + /// tangential openings + CohesiveInternalField tangential_opening_norm; + /// maximum displacement CohesiveInternalField delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; /// tell if the previous opening state is needed (in iterative schemes) bool use_previous_opening; /// damage CohesiveInternalField damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField sigma_c; /// critical displacement Real delta_c; /// array to temporarily store the normals CohesiveInternalField normals; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ -#include "material_cohesive_inline_impl.hh" #include "cohesive_internal_field_tmpl.hh" +#include "material_cohesive_inline_impl.hh" #endif /* AKANTU_MATERIAL_COHESIVE_HH_ */