diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake
index da63bcf54..de704fb3e 100644
--- a/packages/cohesive_element.cmake
+++ b/packages/cohesive_element.cmake
@@ -1,135 +1,139 @@
#===============================================================================
# @file cohesive_element.cmake
#
# @author Mauro Corrado <mauro.corrado@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
# @author Marco Vocialta <marco.vocialta@epfl.ch>
#
# @date creation: Tue Oct 16 2012
# @date last modification: Tue Jan 12 2016
#
# @brief package description for cohesive elements
#
# @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 <http://www.gnu.org/licenses/>.
#
#===============================================================================
package_declare(cohesive_element
DESCRIPTION "Use cohesive_element package of Akantu"
DEPENDS lapack solid_mechanics)
package_declare_sources(cohesive_element
fe_engine/cohesive_element.hh
fe_engine/fe_engine_template_cohesive.cc
fe_engine/shape_cohesive.hh
fe_engine/shape_cohesive_inline_impl.cc
mesh_utils/cohesive_element_inserter.cc
mesh_utils/cohesive_element_inserter.hh
mesh_utils/cohesive_element_inserter_inline_impl.cc
mesh_utils/cohesive_element_inserter_parallel.cc
model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc
model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh
model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc
model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field_tmpl.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.hh
- model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_tmpl.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.cc
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.hh
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.cc
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_includes.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_inline_impl.cc
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.cc
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc
)
package_declare_elements(cohesive_element
ELEMENT_TYPES
_cohesive_1d_2
_cohesive_2d_4
_cohesive_2d_6
_cohesive_3d_12
_cohesive_3d_16
_cohesive_3d_6
_cohesive_3d_8
KIND cohesive
GEOMETRICAL_TYPES
_gt_cohesive_1d_2
_gt_cohesive_2d_4
_gt_cohesive_2d_6
_gt_cohesive_3d_12
_gt_cohesive_3d_16
_gt_cohesive_3d_6
_gt_cohesive_3d_8
FE_ENGINE_LISTS
compute_normals_on_integration_points
contains
get_shapes_derivatives
gradient_on_integration_points
interpolate_on_integration_points
inverse_map
lagrange_base
)
package_declare_material_infos(cohesive_element
LIST AKANTU_COHESIVE_MATERIAL_LIST
INCLUDE material_cohesive_includes.hh
)
package_declare_documentation_files(cohesive_element
manual-cohesive_elements.tex
manual-cohesive_elements_insertion.tex
manual-cohesive_laws.tex
manual-appendix-materials-cohesive.tex
figures/cohesive2d.pdf
figures/cohesive_exponential.pdf
figures/linear_cohesive_law.pdf
figures/bilinear_cohesive_law.pdf
)
package_declare_documentation(cohesive_element
"This package activates the cohesive elements engine within Akantu."
"It depends on:"
"\\begin{itemize}"
" \\item A fortran compiler."
" \\item An implementation of BLAS/LAPACK."
"\\end{itemize}"
)
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 c3f4c5a98..d6ec686f0 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,422 +1,461 @@
/**
* @file material_cohesive_linear.cc
*
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @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
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive_linear.hh"
#include "dof_synchronizer.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
#include <algorithm>
#include <numeric>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
MaterialCohesiveLinear<spatial_dimension>::MaterialCohesiveLinear(
SolidMechanicsModel & model, const ID & id)
: MaterialCohesive(model, id), sigma_c_eff("sigma_c_eff", *this),
delta_c_eff("delta_c_eff", *this),
- insertion_stress("insertion_stress", *this) {
+ insertion_stress("insertion_stress", *this),
+ insertion_compression("insertion_compression", *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,
+ "contact_after_breaking", contact_after_breaking, true,
_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 <UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
MaterialCohesive::initMaterial();
sigma_c_eff.initialize(1);
delta_c_eff.initialize(1);
insertion_stress.initialize(spatial_dimension);
+ insertion_compression.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 <UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::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 <UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::scaleInsertionTraction() {
AKANTU_DEBUG_IN();
// do nothing if volume_s hasn't been specified by the user
if (Math::are_float_equal(volume_s, 0.))
return;
const Mesh & mesh_facets = model->getMeshFacets();
const 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<std::vector<Element>> & 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> & 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<Real> unit_vector(nb_quadrature_points, 1);
volume += fe_engine.integrate(unit_vector, elem->type, elem->element,
elem->ghost_type);
}
// scale sigma_c
*sigma_c_iterator -= base_sigma_c;
*sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s);
*sigma_c_iterator += base_sigma_c;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::checkInsertion(
bool check_only) {
AKANTU_DEBUG_IN();
- const Mesh & mesh_facets = model->getMeshFacets();
- CohesiveElementInserter & inserter = model->getElementInserter();
+ // get mesh and inserter
+ Mesh & mesh = this->model->getMesh();
+ const Mesh & mesh_facets = this->model->getMeshFacets();
+ CohesiveElementInserter & inserter = this->model->getElementInserter();
+
+ // get nodal fields
+ auto & shift = mesh.getNodalData<Real>("shift", spatial_dimension);
+ shift.resize(mesh.getNbNodes());
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 & ins_comp = insertion_compression(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);
+ const auto & connectivity = mesh_facets.getConnectivity(type_facet);
+ auto conn_it = make_view(connectivity, connectivity.getNbComponent()).begin();
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<Real> stress_tmp(spatial_dimension, spatial_dimension);
Matrix<Real> normal_traction(spatial_dimension, nb_quad_facet);
+ Matrix<Real> normal_compression(spatial_dimension, nb_quad_facet);
Vector<Real> 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<Real> new_sigmas;
std::vector<Vector<Real>> new_normal_traction;
+ std::vector<Vector<Real>> new_normal_compression;
std::vector<Real> new_delta_c;
// loop over each facet belonging to this material
for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) {
UInt facet = f_filter(f);
// skip facets where check shouldn't be realized
if (!facets_check(facet))
continue;
// compute the effective norm on each quadrature point of the facet
+ Vector<Real> avg_compression(spatial_dimension);
for (UInt q = 0; q < nb_quad_facet; ++q) {
UInt current_quad = facet * nb_quad_facet + q;
const Vector<Real> & normal = normal_begin[current_quad];
const Vector<Real> & tangent = tangent_begin[current_quad];
const Matrix<Real> & facet_stress_it = facet_stress_begin[current_quad];
// compute average stress on the current quadrature point
Matrix<Real> stress_1(facet_stress_it.storage(), spatial_dimension,
spatial_dimension);
Matrix<Real> stress_2(facet_stress_it.storage() + sp2,
spatial_dimension, spatial_dimension);
stress_tmp.copy(stress_1);
stress_tmp += stress_2;
stress_tmp /= 2.;
Vector<Real> normal_traction_vec(normal_traction(q));
-
// compute normal and effective stress
stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent,
normal_traction_vec);
+ // compute compression
+ normal_compression(q) = normal;
+ normal_compression(q) *= std::min(0., normal_traction_vec.dot(normal));
+ avg_compression += std::min(0., normal_traction_vec.dot(normal))*normal;
+ // remove compression from traction
+ normal_traction_vec -=
+ std::min(0., normal_traction_vec.dot(normal))*normal;
}
+ avg_compression /= nb_quad_facet;
// 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<Real> normal_traction_vec(normal_traction(q));
- if (spatial_dimension != 3)
+ if (spatial_dimension != 3) {
normal_traction_vec *= -1.;
-
+ normal_compression(q) *= -1.;
+ }
new_sigmas.push_back(new_sigma);
new_normal_traction.push_back(normal_traction_vec);
+ new_normal_compression.push_back(normal_compression(q));
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);
}
+ if (spatial_dimension != 3)
+ avg_compression *= -1.;
+
+ // set initial interpenetration in case of compression
+ Vector<UInt> conn = conn_it[f];
+ for (auto && node : conn) {
+ for (UInt dim = 0; dim < spatial_dimension; ++dim) {
+ shift(node, dim) = avg_compression(dim) / this->penalty;
+ }
+ }
}
}
// 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);
+ ins_comp.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);
+ ins_comp(old_nb_quad_points + q, dim) = new_normal_compression[q](dim);
trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::computeTraction(
const Array<Real> & 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_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);
+ auto insertion_compression_it =
+ insertion_compression(el_type, ghost_type).begin(spatial_dimension);
Vector<Real> normal_opening(spatial_dimension);
Vector<Real> 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) {
Real normal_opening_norm, tangential_opening_norm;
bool penetration;
+
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, tangential_opening_norm, *damage_it, penetration,
+ *insertion_stress_it, *insertion_compression_it, *sigma_c_it,
+ normal_opening, tangential_opening, normal_opening_norm,
+ tangential_opening_norm, *damage_it, penetration,
*contact_traction_it, *contact_opening_it);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialCohesiveLinear<spatial_dimension>::computeTangentTraction(
const ElementType & el_type, Array<Real> & tangent_matrix,
const Array<Real> & 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);
/// 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<Real> normal_opening(spatial_dimension);
Vector<Real> 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, tangential_opening_norm;
bool penetration;
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);
}
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.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh
index a954a6953..05e5e5dad 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh
@@ -1,188 +1,191 @@
/**
* @file material_cohesive_linear.hh
*
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Linear irreversible cohesive law of mixed mode loading with
* random stress definition for extrinsic type
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__
#define __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__
namespace akantu {
/**
* Cohesive material linear damage for extrinsic case
*
* parameters in the material files :
* - sigma_c : critical stress sigma_c (default: 0)
* - beta : weighting parameter for sliding and normal opening (default:
* 0)
* - G_cI : fracture energy for mode I (default: 0)
* - G_cII : fracture energy for mode II (default: 0)
* - penalty : stiffness in compression to prevent penetration
*/
template <UInt spatial_dimension>
class MaterialCohesiveLinear : public MaterialCohesive {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialCohesiveLinear(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material parameters
void initMaterial() override;
void updateInternalParameters() override;
/// check stress for cohesive elements' insertion
void checkInsertion(bool check_only = false) override;
/// compute effective stress norm for insertion check
Real computeEffectiveNorm(const Matrix<Real> & stress,
const Vector<Real> & normal,
const Vector<Real> & tangent,
Vector<Real> & normal_stress) const;
protected:
/// constitutive law
void computeTraction(const Array<Real> & normal, ElementType el_type,
GhostType ghost_type = _not_ghost) override;
/// compute tangent stiffness matrix
void computeTangentTraction(const ElementType & el_type,
Array<Real> & tangent_matrix,
const Array<Real> & normal,
GhostType ghost_type) override;
/**
* Scale insertion traction sigma_c according to the volume of the
* two elements surrounding a facet
*
* see the article: F. Zhou and J. F. Molinari "Dynamic crack
* propagation with cohesive elements: a methodology to address mesh
* dependency" International Journal for Numerical Methods in
* Engineering (2004)
*/
void scaleInsertionTraction();
/// compute the traction for a given quadrature point
inline void computeTractionOnQuad(
Vector<Real> & traction, Vector<Real> & opening,
const Vector<Real> & normal, Real & delta_max, const Real & delta_c,
- const Vector<Real> & insertion_stress, const Real & sigma_c,
+ const Vector<Real> & insertion_stress,
+ const Vector<Real> & insertion_compression, const Real & sigma_c,
Vector<Real> & normal_opening, Vector<Real> & tangential_opening,
- Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage,
- bool & penetration, Vector<Real> & contact_traction,
+ Real & normal_opening_norm, Real & tangential_opening_norm,
+ Real & damage, bool & penetration, Vector<Real> & contact_traction,
Vector<Real> & contact_opening);
inline void computeTangentTractionOnQuad(
Matrix<Real> & tangent, Real & delta_max, const Real & delta_c,
const Real & sigma_c, Vector<Real> & opening, const Vector<Real> & normal,
Vector<Real> & normal_opening, Vector<Real> & tangential_opening,
Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage,
bool & penetration, Vector<Real> & contact_opening);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get sigma_c_eff
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(InsertionTraction, sigma_c_eff, Real);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// beta parameter
Real beta;
/// beta square inverse to compute effective norm
Real beta2_inv;
/// mode I fracture energy
Real G_c;
/// kappa parameter
Real kappa;
/// constitutive law scalar to compute delta
Real beta2_kappa2;
/// constitutive law scalar to compute traction
Real beta2_kappa;
/// penalty coefficient
Real penalty;
/// reference volume used to scale sigma_c
Real volume_s;
/// weibull exponent used to scale sigma_c
Real m_s;
/// variable defining if we are recomputing the last loading step
/// after load_reduction
bool recompute;
/// critical effective stress
RandomInternalField<Real, CohesiveInternalField> sigma_c_eff;
/// effective critical displacement (each element can have a
/// different value)
CohesiveInternalField<Real> delta_c_eff;
/// stress at insertion
CohesiveInternalField<Real> insertion_stress;
+ /// compression at insertion
+ CohesiveInternalField<Real> insertion_compression;
/// variable saying if there should be penalty contact also after
/// breaking the cohesive elements
bool contact_after_breaking;
/// insertion of cohesive element when stress is high enough just on
/// one quadrature point
bool max_quad_stress_insertion;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
} // namespace akantu
#include "material_cohesive_linear_inline_impl.cc"
#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc
deleted file mode 100644
index eec417b08..000000000
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc
+++ /dev/null
@@ -1,277 +0,0 @@
-/**
- * @file material_cohesive_linear_friction.cc
- *
- * @author Mauro Corrado <mauro.corrado@epfl.ch>
- *
- * @date creation: Tue Jan 12 2016
- * @date last modification: Wed Feb 21 2018
- *
- * @brief Linear irreversible cohesive law of mixed mode loading with
- * random stress definition for extrinsic type
- *
- * @section LICENSE
- *
- * 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 <http://www.gnu.org/licenses/>.
- *
- */
-
-/* -------------------------------------------------------------------------- */
-#include "material_cohesive_linear_friction.hh"
-#include "solid_mechanics_model_cohesive.hh"
-
-namespace akantu {
-
-/* -------------------------------------------------------------------------- */
-template <UInt spatial_dimension>
-MaterialCohesiveLinearFriction<spatial_dimension>::
- MaterialCohesiveLinearFriction(SolidMechanicsModel & model, const ID & id)
- : MaterialParent(model, id), residual_sliding("residual_sliding", *this),
- friction_force("friction_force", *this) {
- AKANTU_DEBUG_IN();
-
- this->registerParam("mu", mu_max, Real(0.), _pat_parsable | _pat_readable,
- "Maximum value of the friction coefficient");
-
- this->registerParam("penalty_for_friction", friction_penalty, Real(0.),
- _pat_parsable | _pat_readable,
- "Penalty parameter for the friction behavior");
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-template <UInt spatial_dimension>
-void MaterialCohesiveLinearFriction<spatial_dimension>::initMaterial() {
- AKANTU_DEBUG_IN();
-
- MaterialParent::initMaterial();
-
- friction_force.initialize(spatial_dimension);
- residual_sliding.initialize(1);
- residual_sliding.initializeHistory();
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-template <UInt spatial_dimension>
-void MaterialCohesiveLinearFriction<spatial_dimension>::computeTraction(
- __attribute__((unused)) const Array<Real> & normal, ElementType el_type,
- GhostType ghost_type) {
- AKANTU_DEBUG_IN();
-
- residual_sliding.resize();
- friction_force.resize();
-
- /// define iterators
- auto traction_it =
- this->tractions(el_type, ghost_type).begin(spatial_dimension);
- auto traction_end =
- this->tractions(el_type, ghost_type).end(spatial_dimension);
- auto opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension);
- auto previous_opening_it =
- this->opening.previous(el_type, ghost_type).begin(spatial_dimension);
- auto contact_traction_it =
- this->contact_tractions(el_type, ghost_type).begin(spatial_dimension);
- auto contact_opening_it =
- this->contact_opening(el_type, ghost_type).begin(spatial_dimension);
- auto normal_it = this->normal.begin(spatial_dimension);
- auto sigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin();
- auto delta_max_it = this->delta_max(el_type, ghost_type).begin();
- auto delta_max_prev_it =
- this->delta_max.previous(el_type, ghost_type).begin();
- auto delta_c_it = this->delta_c_eff(el_type, ghost_type).begin();
- auto damage_it = this->damage(el_type, ghost_type).begin();
- auto insertion_stress_it =
- this->insertion_stress(el_type, ghost_type).begin(spatial_dimension);
- auto res_sliding_it = this->residual_sliding(el_type, ghost_type).begin();
- auto res_sliding_prev_it =
- this->residual_sliding.previous(el_type, ghost_type).begin();
- auto friction_force_it =
- this->friction_force(el_type, ghost_type).begin(spatial_dimension);
-
- Vector<Real> normal_opening(spatial_dimension);
- Vector<Real> tangential_opening(spatial_dimension);
-
- if (not this->model->isDefaultSolverExplicit())
- 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, ++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, ++res_sliding_it,
- ++res_sliding_prev_it, ++friction_force_it, ++previous_opening_it) {
-
- Real normal_opening_norm, tangential_opening_norm;
- bool penetration;
- 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, tangential_opening_norm, *damage_it, penetration,
- *contact_traction_it, *contact_opening_it);
-
- if (penetration) {
- /// the friction coefficient mu increases with the damage. It
- /// equals the maximum value when damage = 1.
- // Real damage = std::min(*delta_max_prev_it / *delta_c_it,
- // Real(1.));
- Real mu = mu_max; // * damage;
-
- /// the definition of tau_max refers to the opening
- /// (penetration) of the previous incremental step
- Real normal_opening_prev_norm =
- std::min(previous_opening_it->dot(*normal_it), Real(0.));
-
- // Vector<Real> normal_opening_prev = (*normal_it);
- // normal_opening_prev *= normal_opening_prev_norm;
- Real tau_max = mu * this->penalty * (std::abs(normal_opening_prev_norm));
- Real delta_sliding_norm =
- std::abs(tangential_opening_norm - *res_sliding_prev_it);
-
- /// tau is the norm of the friction force, acting tangentially to the
- /// surface
- Real tau = std::min(friction_penalty * delta_sliding_norm, tau_max);
-
- if ((tangential_opening_norm - *res_sliding_prev_it) < 0.0)
- tau = -tau;
-
- /// from tau get the x and y components of friction, to be added in the
- /// force vector
- Vector<Real> tangent_unit_vector(spatial_dimension);
- tangent_unit_vector = tangential_opening / tangential_opening_norm;
- *friction_force_it = tau * tangent_unit_vector;
-
- /// update residual_sliding
- *res_sliding_it =
- tangential_opening_norm - (std::abs(tau) / friction_penalty);
- } else {
- friction_force_it->clear();
- }
-
- *traction_it += *friction_force_it;
- }
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-template <UInt spatial_dimension>
-void MaterialCohesiveLinearFriction<spatial_dimension>::computeTangentTraction(
- const ElementType & el_type, Array<Real> & tangent_matrix,
- __attribute__((unused)) const Array<Real> & 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 = this->normal.begin(spatial_dimension);
-
- auto opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension);
- auto previous_opening_it =
- this->opening.previous(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
- */
- auto delta_max_it = this->delta_max.previous(el_type, ghost_type).begin();
- auto sigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin();
- auto delta_c_it = this->delta_c_eff(el_type, ghost_type).begin();
- auto damage_it = this->damage(el_type, ghost_type).begin();
-
- auto contact_opening_it =
- this->contact_opening(el_type, ghost_type).begin(spatial_dimension);
-
- auto res_sliding_prev_it =
- this->residual_sliding.previous(el_type, ghost_type).begin();
-
- Vector<Real> normal_opening(spatial_dimension);
- Vector<Real> tangential_opening(spatial_dimension);
-
- for (; tangent_it != tangent_end;
- ++tangent_it, ++normal_it, ++opening_it, ++previous_opening_it,
- ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it,
- ++contact_opening_it, ++res_sliding_prev_it) {
-
- Real normal_opening_norm, tangential_opening_norm;
- bool penetration;
- 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);
-
- if (penetration) {
- // Real damage = std::min(*delta_max_it / *delta_c_it, Real(1.));
- Real mu = mu_max; // * damage;
-
- Real normal_opening_prev_norm =
- std::min(previous_opening_it->dot(*normal_it), Real(0.));
- // Vector<Real> normal_opening_prev = (*normal_it);
- // normal_opening_prev *= normal_opening_prev_norm;
-
- Real tau_max = mu * this->penalty * (std::abs(normal_opening_prev_norm));
- Real delta_sliding_norm =
- std::abs(tangential_opening_norm - *res_sliding_prev_it);
-
- // tau is the norm of the friction force, acting tangentially to the
- // surface
- Real tau = std::min(friction_penalty * delta_sliding_norm, tau_max);
-
- if (tau < tau_max && tau_max > Math::getTolerance()) {
- Matrix<Real> I(spatial_dimension, spatial_dimension);
- I.eye(1.);
-
- Matrix<Real> n_outer_n(spatial_dimension, spatial_dimension);
- n_outer_n.outerProduct(*normal_it, *normal_it);
-
- Matrix<Real> nn(n_outer_n);
- I -= nn;
- *tangent_it += I * friction_penalty;
- }
- }
-
- // 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(cohesive_linear_friction, MaterialCohesiveLinearFriction);
-
-} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
index 780ce1922..14936e21f 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
@@ -1,105 +1,103 @@
/**
* @file material_cohesive_linear_friction.hh
*
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Linear irreversible cohesive law of mixed mode loading with
* random stress definition for extrinsic type
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive_linear.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__
#define __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
/**
* Cohesive material linear with friction force
*
* parameters in the material files :
* - mu : friction coefficient
* - penalty_for_friction : Penalty parameter for the friction behavior
*/
template <UInt spatial_dimension>
class MaterialCohesiveLinearFriction
: public MaterialCohesiveLinear<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
using MaterialParent = MaterialCohesiveLinear<spatial_dimension>;
public:
MaterialCohesiveLinearFriction(SolidMechanicsModel & model,
const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material parameters
void initMaterial() override;
protected:
- /// constitutive law
- void computeTraction(const Array<Real> & normal, ElementType el_type,
- GhostType ghost_type = _not_ghost) override;
-
- /// compute tangent stiffness matrix
- void computeTangentTraction(const ElementType & el_type,
- Array<Real> & tangent_matrix,
- const Array<Real> & normal,
- GhostType ghost_type) override;
+ /// check stress for cohesive elements insertion
+ void checkInsertion(bool check_only = false) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
- /// maximum value of the friction coefficient
- Real mu_max;
-
/// penalty parameter for the friction law
Real friction_penalty;
+ /// friction at insertion
+ RandomInternalField<Real, FacetInternalField> mu_insertion;
+
+ /// friction coefficient
+ CohesiveInternalField<Real> mu_eff;
+
/// history parameter for the friction law
CohesiveInternalField<Real> residual_sliding;
/// friction force
CohesiveInternalField<Real> friction_force;
};
} // namespace akantu
+#include "material_cohesive_linear_friction_tmpl.hh"
+
#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.cc
new file mode 100644
index 000000000..0ff6dac87
--- /dev/null
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.cc
@@ -0,0 +1,251 @@
+/**
+ * @file material_cohesive_linear_friction_coulomb.cc
+ *
+ * @author Mauro Corrado <mauro.corrado@epfl.ch>
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Mathias Lebihain <mathias.lebihain@epfl.ch>
+ *
+ * @date creation: Tue Jan 12 2016
+ * @date last modification: Wed Feb 21 2018
+ *
+ * @brief Linear irreversible cohesive law of mixed mode loading with
+ * random stress definition for extrinsic type
+ *
+ * @section LICENSE
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "material_cohesive_linear_friction_coulomb.hh"
+#include "solid_mechanics_model_cohesive.hh"
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+MaterialCohesiveLinearFrictionCoulomb<
+ dim>::MaterialCohesiveLinearFrictionCoulomb(SolidMechanicsModel & model,
+ const ID & id)
+ : MaterialParent(model, id) {
+ AKANTU_DEBUG_IN();
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+/** \todo Frictional contributions might be computed from the contact_tractions,
+ * and delta_max at the previous time step
+ */
+template <UInt dim>
+void MaterialCohesiveLinearFrictionCoulomb<dim>::computeTraction(
+ const Array<Real> & normal, ElementType el_type, GhostType ghost_type) {
+ AKANTU_DEBUG_IN();
+
+ Vector<Real> normal_opening(dim);
+ Vector<Real> tangential_opening(dim);
+
+ /// loop on each quadrature point
+ for (auto && data :
+ zip(make_view(this->insertion_stress(el_type, ghost_type), dim),
+ make_view(this->insertion_compression(el_type, ghost_type), dim),
+ this->sigma_c_eff(el_type, ghost_type),
+ this->delta_c_eff(el_type, ghost_type),
+ this->mu_eff(el_type, ghost_type), make_view(normal, dim),
+ this->delta_max(el_type, ghost_type),
+ this->damage(el_type, ghost_type),
+ make_view(this->tractions(el_type, ghost_type), dim),
+ make_view(this->contact_tractions(el_type, ghost_type), dim),
+ make_view(this->friction_force(el_type, ghost_type), dim),
+ make_view(this->opening(el_type, ghost_type), dim),
+ make_view(this->contact_opening(el_type, ghost_type), dim),
+ this->residual_sliding(el_type, ghost_type),
+ this->residual_sliding.previous(el_type, ghost_type))) {
+ auto && insertion_traction = std::get<0>(data);
+ auto && insertion_compression = std::get<1>(data);
+ auto && sigma_c = std::get<2>(data);
+ auto && delta_c = std::get<3>(data);
+ auto && mu = std::get<4>(data);
+ auto && normal = std::get<5>(data);
+ auto && delta_max = std::get<6>(data);
+ auto && damage = std::get<7>(data);
+ auto && traction = std::get<8>(data);
+ auto && contact_traction = std::get<9>(data);
+ auto && friction_force = std::get<10>(data);
+ auto && opening = std::get<11>(data);
+ auto && contact_opening = std::get<12>(data);
+ auto && res_sliding = std::get<13>(data);
+ auto && res_sliding_prev = std::get<14>(data);
+
+ Real normal_opening_norm, tangential_opening_norm;
+ bool penetration;
+
+ this->computeTractionOnQuad(
+ traction, opening, normal, delta_max, delta_c, insertion_traction,
+ insertion_compression, sigma_c, normal_opening, tangential_opening,
+ normal_opening_norm, tangential_opening_norm, damage, penetration,
+ contact_traction, contact_opening);
+
+ if (penetration) {
+ /// the definition of tau_max refers to
+ /// the current contact_traction computed in computeTractionOnQuad
+ Real tau_max = mu * contact_traction.norm();
+
+ // elastic sliding
+ Real delta_sliding_norm =
+ std::abs(tangential_opening_norm - res_sliding_prev);
+
+ /// tau is the norm of the friction force, acting tangentially to the
+ /// surface
+ Real tau = std::min(this->friction_penalty * delta_sliding_norm, tau_max);
+
+ if ((tangential_opening_norm - res_sliding_prev) < 0.0)
+ tau = -tau;
+
+ /// from tau get the x and y components of friction, to be added in the
+ /// force vector
+ auto tangent_unit_vector = tangential_opening / tangential_opening_norm;
+ friction_force = tau * tangent_unit_vector;
+
+ /// update residual_sliding
+ res_sliding =
+ tangential_opening_norm - (std::abs(tau) / this->friction_penalty);
+ } else {
+ friction_force.clear();
+ }
+ traction += friction_force;
+ }
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+/** \todo The frictional contribution to the tangent stiffness matrix might be
+ * missing from implicit computations
+ */
+// template <UInt dim>
+// void MaterialCohesiveLinearFrictionCoulomb<dim>::computeTangentTraction(
+// const ElementType & el_type, Array<Real> & tangent_matrix,
+// __attribute__((unused)) const Array<Real> & normal, GhostType ghost_type)
+// {
+// AKANTU_DEBUG_IN();
+//
+// /// define iterators
+// auto tangent_it = tangent_matrix.begin(dim, dim);
+// auto tangent_end = tangent_matrix.end(dim, dim);
+//
+// auto normal_it = this->normal.begin(dim);
+//
+// auto opening_it = this->opening(el_type, ghost_type).begin(dim);
+//
+// /**
+// * 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 = this->delta_max.previous(el_type, ghost_type).begin();
+// auto sigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin();
+// auto delta_c_it = this->delta_c_eff(el_type, ghost_type).begin();
+// auto damage_it = this->damage(el_type, ghost_type).begin();
+// auto contact_traction_it =
+// this->contact_tractions(el_type, ghost_type).begin(dim);
+// auto contact_opening_it =
+// this->contact_opening(el_type, ghost_type).begin(dim);
+// auto res_sliding_prev_it =
+// this->residual_sliding.previous(el_type, ghost_type).begin();
+//
+// Vector<Real> normal_opening(dim);
+// Vector<Real> tangential_opening(dim);
+//
+// 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,
+// ++res_sliding_prev_it,
+// ++contact_traction_it) {
+//
+// Real normal_opening_norm, tangential_opening_norm;
+// bool penetration;
+//
+// 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);
+//
+// if (penetration) {
+// // Real damage = std::min(*delta_max_it / *delta_c_it, Real(1.));
+// // Real mu = mu_max * damage;
+// Real mu = mu_max;
+//
+// // Real normal_opening_norm =
+// // std::min(contact_opening_it->dot(*normal_it), Real(0.));
+// // Real tau_max = mu * this->penalty * (std::abs(normal_opening_norm));
+//
+// Real tau_max = mu * (*contact_traction_it).norm();
+//
+// // elastic sliding
+// Real delta_sliding_norm =
+// std::abs(tangential_opening_norm - *res_sliding_prev_it);
+//
+// /// tau is the norm of the friction force, acting tangentially to the
+// /// surface
+// Real tau = std::min(this->friction_penalty * delta_sliding_norm,
+// tau_max);
+//
+// if (tau < tau_max && tau_max > Math::getTolerance()) {
+// Matrix<Real> I(dim, dim);
+// I.eye(1.);
+//
+// Matrix<Real> n_outer_n(dim, dim);
+// n_outer_n.outerProduct(*normal_it, *normal_it);
+//
+// Matrix<Real> nn(n_outer_n);
+// I -= nn;
+// *tangent_it += I * this->friction_penalty;
+// }
+// }
+//
+// // check if the tangential stiffness matrix is symmetric
+// // for (UInt h = 0; h < dim; ++h){
+// // for (UInt l = h; l < dim; ++l){
+// // if (l > h){
+// // Real k_ls = (*tangent_it)[dim*h+l];
+// // Real k_us = (*tangent_it)[dim*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(cohesive_linear_friction_coulomb,
+ MaterialCohesiveLinearFrictionCoulomb);
+
+} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.hh
similarity index 53%
copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.hh
index 780ce1922..82a5f51bc 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_coulomb.hh
@@ -1,105 +1,78 @@
/**
- * @file material_cohesive_linear_friction.hh
+ * @file material_cohesive_linear_friction_coulomb.hh
*
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Mathias Lebihain <mathias.lebihain@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Linear irreversible cohesive law of mixed mode loading with
* random stress definition for extrinsic type
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
-#include "material_cohesive_linear.hh"
+#include "material_cohesive_linear_friction.hh"
/* -------------------------------------------------------------------------- */
-#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__
-#define __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__
+#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_COULOMB_HH__
+#define __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_COULOMB_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
/**
- * Cohesive material linear with friction force
- *
- * parameters in the material files :
- * - mu : friction coefficient
- * - penalty_for_friction : Penalty parameter for the friction behavior
+ * Cohesive material linear with coulomb friction force
*/
template <UInt spatial_dimension>
-class MaterialCohesiveLinearFriction
- : public MaterialCohesiveLinear<spatial_dimension> {
+class MaterialCohesiveLinearFrictionCoulomb
+ : public MaterialCohesiveLinearFriction<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
- using MaterialParent = MaterialCohesiveLinear<spatial_dimension>;
+ using MaterialParent = MaterialCohesiveLinearFriction<spatial_dimension>;
public:
- MaterialCohesiveLinearFriction(SolidMechanicsModel & model,
- const ID & id = "");
+ MaterialCohesiveLinearFrictionCoulomb(SolidMechanicsModel & model,
+ const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
-public:
- /// initialize the material parameters
- void initMaterial() override;
-
protected:
/// constitutive law
void computeTraction(const Array<Real> & normal, ElementType el_type,
GhostType ghost_type = _not_ghost) override;
- /// compute tangent stiffness matrix
- void computeTangentTraction(const ElementType & el_type,
- Array<Real> & tangent_matrix,
- const Array<Real> & normal,
- GhostType ghost_type) override;
-
- /* ------------------------------------------------------------------------ */
- /* Accessors */
- /* ------------------------------------------------------------------------ */
-public:
- /* ------------------------------------------------------------------------ */
- /* Class Members */
- /* ------------------------------------------------------------------------ */
-protected:
- /// maximum value of the friction coefficient
- Real mu_max;
-
- /// penalty parameter for the friction law
- Real friction_penalty;
-
- /// history parameter for the friction law
- CohesiveInternalField<Real> residual_sliding;
-
- /// friction force
- CohesiveInternalField<Real> friction_force;
+ // /// compute tangent stiffness matrix
+ // void computeTangentTraction(const ElementType & el_type,
+ // Array<Real> & tangent_matrix,
+ // const Array<Real> & normal,
+ // GhostType ghost_type) override;
};
} // namespace akantu
-#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__ */
+#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_COULOMB_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.cc
new file mode 100644
index 000000000..8111920ce
--- /dev/null
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.cc
@@ -0,0 +1,165 @@
+/**
+ * @file material_cohesive_linear_slipweakening.cc
+ *
+ * @author Mathias Lebihain <mathias.lebihain@epfl.ch>
+ *
+ * @date creation: Fri Mar 06 2020
+ * @date last modification: Fri Mar 06 2020
+ *
+ * @brief Linear irreversible cohesive law of mixed mode loading with
+ * Mohr-Coulomb insertion criterion and slip-weakening friction for
+ * extrinsic type
+ *
+ * @section LICENSE
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "material_cohesive_linear_friction_slipweakening.hh"
+#include "solid_mechanics_model_cohesive.hh"
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+MaterialCohesiveLinearFrictionSlipWeakening<dim>::
+ MaterialCohesiveLinearFrictionSlipWeakening(SolidMechanicsModel & model,
+ const ID & id)
+ : MaterialParent(model, id), mu_dynamic("mu_dynamic", *this) {
+ AKANTU_DEBUG_IN();
+
+ this->registerParam("mu_static", this->mu_insertion,
+ _pat_parsable | _pat_readable,
+ "Static value of the friction coefficient");
+
+ this->registerParam("mu_dynamic", this->mu_dynamic,
+ _pat_parsable | _pat_readable,
+ "Dynamic value of the friction coefficient");
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+void MaterialCohesiveLinearFrictionSlipWeakening<dim>::initMaterial() {
+ AKANTU_DEBUG_IN();
+
+ MaterialParent::initMaterial();
+
+ mu_dynamic.initialize(1);
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+/** \todo Frictional contributions might be computed from the contact_tractions,
+ * and delta_max at the previous time step
+ */
+template <UInt dim>
+void MaterialCohesiveLinearFrictionSlipWeakening<dim>::computeTraction(
+ const Array<Real> & normal, ElementType el_type, GhostType ghost_type) {
+ AKANTU_DEBUG_IN();
+
+ Vector<Real> normal_opening(dim);
+ Vector<Real> tangential_opening(dim);
+
+ /// loop on each quadrature point
+ for (auto && data :
+ zip(make_view(this->insertion_stress(el_type, ghost_type), dim),
+ make_view(this->insertion_compression(el_type, ghost_type), dim),
+ this->sigma_c_eff(el_type, ghost_type),
+ this->delta_c_eff(el_type, ghost_type),
+ this->mu_eff(el_type, ghost_type),
+ this->mu_dynamic(el_type, ghost_type), make_view(normal, dim),
+ this->delta_max(el_type, ghost_type),
+ this->damage(el_type, ghost_type),
+ make_view(this->tractions(el_type, ghost_type), dim),
+ make_view(this->contact_tractions(el_type, ghost_type), dim),
+ make_view(this->friction_force(el_type, ghost_type), dim),
+ make_view(this->opening(el_type, ghost_type), dim),
+ make_view(this->contact_opening(el_type, ghost_type), dim),
+ this->residual_sliding(el_type, ghost_type),
+ this->residual_sliding.previous(el_type, ghost_type))) {
+
+ auto && insertion_traction = std::get<0>(data);
+ auto && insertion_compression = std::get<1>(data);
+ auto && sigma_c = std::get<2>(data);
+ auto && delta_c = std::get<3>(data);
+ auto && mu_static = std::get<4>(data);
+ auto && mu_dynamic = std::get<5>(data);
+ auto && normal = std::get<6>(data);
+ auto && delta_max = std::get<7>(data);
+ auto && damage = std::get<8>(data);
+ auto && traction = std::get<9>(data);
+ auto && contact_traction = std::get<10>(data);
+ auto && friction_force = std::get<11>(data);
+ auto && opening = std::get<12>(data);
+ auto && contact_opening = std::get<13>(data);
+ auto && res_sliding = std::get<14>(data);
+ auto && res_sliding_prev = std::get<15>(data);
+
+ Real normal_opening_norm, tangential_opening_norm;
+ bool penetration;
+
+ this->computeTractionOnQuad(
+ traction, opening, normal, delta_max, delta_c, insertion_traction,
+ insertion_compression, sigma_c, normal_opening, tangential_opening,
+ normal_opening_norm, tangential_opening_norm, damage, penetration,
+ contact_traction, contact_opening);
+
+ if (penetration) {
+ /// current friction coefficient
+ Real mu = mu_static + (mu_dynamic - mu_static) * damage;
+ /// the definition of tau_max refers to
+ /// the current contact_traction computed in computeTractionOnQuad
+ Real tau_max = mu * contact_traction.norm();
+
+ // elastic sliding
+ Real delta_sliding_norm =
+ std::abs(tangential_opening_norm - res_sliding_prev);
+
+ /// tau is the norm of the friction force, acting tangentially to the
+ /// surface
+ Real tau = std::min(this->friction_penalty * delta_sliding_norm, tau_max);
+
+ if ((tangential_opening_norm - res_sliding_prev) < 0.0)
+ tau = -tau;
+
+ /// from tau get the x and y components of friction, to be added in the
+ /// force vector
+ auto tangent_unit_vector = tangential_opening / tangential_opening_norm;
+ friction_force = tau * tangent_unit_vector;
+
+ /// update residual_sliding
+ res_sliding =
+ tangential_opening_norm - (std::abs(tau) / this->friction_penalty);
+ } else {
+ friction_force.clear();
+ }
+ traction += friction_force;
+ }
+
+ AKANTU_DEBUG_OUT();
+}
+/* -------------------------------------------------------------------------- */
+
+INSTANTIATE_MATERIAL(cohesive_linear_friction_slipweakening,
+ MaterialCohesiveLinearFrictionSlipWeakening);
+
+} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.hh
similarity index 61%
copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.hh
index 780ce1922..3ffbdd44f 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_slipweakening.hh
@@ -1,105 +1,88 @@
/**
- * @file material_cohesive_linear_friction.hh
+ * @file material_cohesive_linear_slipweakening.hh
*
- * @author Mauro Corrado <mauro.corrado@epfl.ch>
- * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Mathias Lebihain <mathias.lebihain@epfl.ch>
*
- * @date creation: Fri Jun 18 2010
- * @date last modification: Wed Feb 21 2018
+ * @date creation: Fri Mar 06 2020
+ * @date last modification: Fri Mar 06 2020
*
* @brief Linear irreversible cohesive law of mixed mode loading with
- * random stress definition for extrinsic type
+ * Mohr-Coulomb insertion criterion and slip-weakening friction for
+ * extrinsic type
*
* @section LICENSE
*
- * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
-#include "material_cohesive_linear.hh"
+#include "material_cohesive_linear_friction.hh"
/* -------------------------------------------------------------------------- */
-#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__
-#define __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__
-
+#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_SLIPWEAKENING_HH__
+#define __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_SLIPWEAKENING_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
/**
- * Cohesive material linear with friction force
+ * Cohesive material linear with slip-weakening friction force
*
* parameters in the material files :
- * - mu : friction coefficient
- * - penalty_for_friction : Penalty parameter for the friction behavior
+ * - mu_d : dynamic friction coefficient
*/
template <UInt spatial_dimension>
-class MaterialCohesiveLinearFriction
- : public MaterialCohesiveLinear<spatial_dimension> {
+class MaterialCohesiveLinearFrictionSlipWeakening
+ : public MaterialCohesiveLinearFriction<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
- using MaterialParent = MaterialCohesiveLinear<spatial_dimension>;
+ using MaterialParent = MaterialCohesiveLinearFriction<spatial_dimension>;
public:
- MaterialCohesiveLinearFriction(SolidMechanicsModel & model,
- const ID & id = "");
+ MaterialCohesiveLinearFrictionSlipWeakening(SolidMechanicsModel & model,
+ const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material parameters
void initMaterial() override;
protected:
/// constitutive law
void computeTraction(const Array<Real> & normal, ElementType el_type,
GhostType ghost_type = _not_ghost) override;
- /// compute tangent stiffness matrix
- void computeTangentTraction(const ElementType & el_type,
- Array<Real> & tangent_matrix,
- const Array<Real> & normal,
- GhostType ghost_type) override;
-
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
- /// maximum value of the friction coefficient
- Real mu_max;
-
- /// penalty parameter for the friction law
- Real friction_penalty;
-
- /// history parameter for the friction law
- CohesiveInternalField<Real> residual_sliding;
-
- /// friction force
- CohesiveInternalField<Real> friction_force;
+ /// dynamic friction coefficient
+ RandomInternalField<Real, CohesiveInternalField> mu_dynamic;
};
} // namespace akantu
-#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__ */
+#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_SLIPWEAKENING_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_tmpl.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_tmpl.hh
new file mode 100644
index 000000000..b5d3d429b
--- /dev/null
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction_tmpl.hh
@@ -0,0 +1,335 @@
+/**
+ * @file material_cohesive_linear_friction.cc
+ *
+ * @author Mauro Corrado <mauro.corrado@epfl.ch>
+ *
+ * @date creation: Tue Jan 12 2016
+ * @date last modification: Wed Feb 21 2018
+ *
+ * @brief Linear irreversible cohesive law of mixed mode loading with
+ * random stress definition for extrinsic type
+ *
+ * @section LICENSE
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "material_cohesive_linear_friction.hh"
+#include "solid_mechanics_model_cohesive.hh"
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+template <UInt spatial_dimension>
+MaterialCohesiveLinearFriction<spatial_dimension>::
+ MaterialCohesiveLinearFriction(SolidMechanicsModel & model, const ID & id)
+ : MaterialParent(model, id), mu_insertion("mu_insertion", *this),
+ mu_eff("mu_eff", *this), residual_sliding("residual_sliding", *this),
+ friction_force("friction_force", *this) {
+ AKANTU_DEBUG_IN();
+
+ this->registerParam("mu", mu_insertion, _pat_parsable | _pat_readable,
+ "Value of the friction coefficient at insertion");
+
+ this->registerParam("penalty_for_friction", friction_penalty, Real(0.),
+ _pat_parsable | _pat_readable,
+ "Penalty parameter for the friction behavior");
+
+ this->use_previous_contact_opening = true;
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+template <UInt spatial_dimension>
+void MaterialCohesiveLinearFriction<spatial_dimension>::initMaterial() {
+ AKANTU_DEBUG_IN();
+
+ MaterialParent::initMaterial();
+
+ mu_insertion.initialize(1);
+ mu_eff.initialize(1);
+ friction_force.initialize(spatial_dimension);
+ residual_sliding.initialize(1);
+ residual_sliding.initializeHistory();
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+template <UInt spatial_dimension>
+void MaterialCohesiveLinearFriction<spatial_dimension>::checkInsertion(
+ bool check_only) {
+ AKANTU_DEBUG_IN();
+
+ // get mesh and inserter
+ const Mesh & mesh_facets = this->model->getMeshFacets();
+ CohesiveElementInserter & inserter = this->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 = this->facet_filter(type_facet);
+ auto & sig_c_eff = this->sigma_c_eff(type_cohesive);
+ auto & mu_eff = this->mu_eff(type_cohesive);
+ auto & del_c_eff = this->delta_c_eff(type_cohesive);
+ auto & ins_stress = this->insertion_stress(type_cohesive);
+ auto & ins_comp = this->insertion_compression(type_cohesive);
+ auto & trac_old = this->tractions.previous(type_cohesive);
+
+ UInt nb_quad_facet = this->model->getFEEngine("FacetsFEEngine")
+ .getNbIntegrationPoints(type_facet);
+
+#ifndef AKANTU_NDEBUG
+ UInt nb_quad_cohesive = this->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;
+
+ // get cohesive stress
+ const auto & sigma_lim = this->sigma_c(type_facet);
+ // get friction at insertion
+ const auto & mu_insert = this->mu_insertion(type_facet);
+ // get stress on facets
+ const auto & f_stress = this->model->getStressOnFacets(type_facet);
+ // get normals and tangents to facets
+ const auto & tangents = this->model->getTangents(type_facet);
+ const auto & normals = this->model->getFEEngine("FacetsFEEngine")
+ .getNormalsOnIntegrationPoints(type_facet);
+ // associated iterators
+ auto sigma_lim_it = sigma_lim.begin();
+ auto mu_insert_it = mu_insert.begin();
+ auto facet_stress_begin =
+ f_stress.begin(spatial_dimension, spatial_dimension * 2);
+ auto normal_begin = normals.begin(spatial_dimension);
+ auto tangent_begin = tangents.begin(tangents.getNbComponent());
+
+ // loop utils
+ UInt sp2 = spatial_dimension * spatial_dimension;
+
+ Matrix<Real> stress_quad(spatial_dimension, spatial_dimension);
+ Matrix<Real> normal_traction_quad(spatial_dimension, nb_quad_facet);
+ Vector<Real> normal_compression_quad(nb_quad_facet);
+ Vector<Real> effective_stress_quad(nb_quad_facet);
+ Vector<Real> mu_insert_quad(nb_quad_facet);
+ Vector<Real> sigma_lim_quad(nb_quad_facet);
+
+ std::vector<Real> new_sigmas;
+ std::vector<Real> new_mus;
+ std::vector<Vector<Real>> new_normal_tractions;
+ std::vector<Vector<Real>> new_normal_compressions;
+ std::vector<Real> new_delta_cs;
+
+ // loop over each facet belonging to this material
+ for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it, ++mu_insert_it) {
+ // facet id
+ UInt facet = f_filter(f);
+
+ // skip facets where check shouldn't be realized
+ if (!facets_check(facet))
+ continue;
+
+ // initialization for average activation
+ Vector<Real> avg_normal(spatial_dimension);
+ Vector<Real> avg_tangent(spatial_dimension);
+ Vector<Real> avg_normal_traction(spatial_dimension);
+ Real avg_normal_compression = 0.;
+ Real mu_insert_avg = 0.;
+ Real sigma_lim_avg = 0.;
+
+ // compute the effective norm on each quadrature point of the facet
+ Real normal_compression_q;
+ for (UInt q = 0; q < nb_quad_facet; ++q) {
+ UInt current_quad = facet * nb_quad_facet + q;
+ const Vector<Real> & normal = normal_begin[current_quad];
+ const Vector<Real> & tangent = tangent_begin[current_quad];
+ const Matrix<Real> & facet_stress_it = facet_stress_begin[current_quad];
+
+ // compute average stress on the current quadrature point
+ Matrix<Real> stress_1(facet_stress_it.storage(), spatial_dimension,
+ spatial_dimension);
+ Matrix<Real> stress_2(facet_stress_it.storage() + sp2,
+ spatial_dimension, spatial_dimension);
+ stress_quad.copy(stress_1);
+ stress_quad += stress_2;
+ stress_quad /= 2.;
+
+ // compute normal and effective stress taking into account the sole
+ // tensile contribution of the facet tractions
+ Vector<Real> normal_traction_vec(normal_traction_quad(q));
+ effective_stress_quad(q) = this->computeEffectiveNorm(
+ stress_quad, normal, tangent, normal_traction_vec);
+
+ // normal traction on quad
+ normal_compression_q = normal_traction_vec.dot(normal);
+ // compressive contribution of the normal traction
+ normal_compression_quad(q) = std::min(0., normal_compression_q);
+ // traction minus compressive part
+ normal_traction_quad(q) =
+ normal_traction_vec - normal * normal_compression_quad(q);
+
+ // store quad and average infos
+ avg_normal_compression += normal_compression_q;
+ avg_normal_traction += normal_traction_vec;
+ avg_normal += normal;
+ avg_tangent += tangent;
+ sigma_lim_avg += *sigma_lim_it;
+ sigma_lim_quad(q) = *sigma_lim_it;
+ mu_insert_avg += *mu_insert_it;
+ mu_insert_quad(q) = *mu_insert_it;
+ }
+ // rescale averages
+ sigma_lim_avg /= nb_quad_facet;
+ mu_insert_avg /= nb_quad_facet;
+ avg_normal /= nb_quad_facet;
+ avg_tangent /= nb_quad_facet;
+ avg_normal_traction /= nb_quad_facet;
+ avg_normal_compression /= nb_quad_facet;
+
+ // compressive contribution of the average normal traction
+ avg_normal_compression = std::min(0., avg_normal_compression);
+ // average traction minus compressive part
+ avg_normal_traction -= avg_normal * avg_normal_compression;
+
+ // verify if the effective stress overcomes the threshold
+ Real crit;
+ if (this->max_quad_stress_insertion) {
+ Vector<Real> crit_quad(nb_quad_facet);
+ for (UInt q = 0; q < nb_quad_facet; ++q) {
+ crit_quad(q) = effective_stress_quad(q) -
+ (sigma_lim_quad(q) - mu_insert_quad(q) / this->beta *
+ normal_compression_quad(q));
+ }
+ crit = *std::max_element(crit_quad.storage(),
+ crit_quad.storage() + nb_quad_facet);
+ // we must then ensure that the average value we set to
+ // new_sigma is above zero
+ if (effective_stress_quad.mean() +
+ mu_insert_avg / this->beta * avg_normal_compression <
+ 0) {
+ crit = effective_stress_quad.mean() -
+ (sigma_lim_avg -
+ mu_insert_avg / this->beta * avg_normal_compression);
+ }
+ } else {
+ crit = effective_stress_quad.mean() -
+ (sigma_lim_avg -
+ mu_insert_avg / this->beta * avg_normal_compression);
+ }
+
+ if (crit > -Math::getTolerance()) {
+ // facet will be inserted
+ f_insertion(facet) = true;
+ // only if insertion is considered
+ if (check_only)
+ continue;
+
+ // add tangential friction part to traction
+ avg_normal_traction +=
+ mu_insert_avg * avg_normal_compression * avg_tangent;
+ if (spatial_dimension != 3) {
+ avg_normal_traction *= -1.;
+ avg_normal_compression *= -1.;
+ }
+
+ // store the new cohesive material parameters for each quadrature point
+ for (UInt q = 0; q < nb_quad_facet; ++q) {
+ // extract tangent
+ UInt current_quad = facet * nb_quad_facet + q;
+ const Vector<Real> & normal = normal_begin[current_quad];
+ const Vector<Real> & tangent = tangent_begin[current_quad];
+
+ // effective stress at insertion
+ Real new_sigma =
+ effective_stress_quad(q) +
+ mu_insert_quad(q) / this->beta * normal_compression_quad(q);
+
+ // friction coefficient
+ Real new_mu = mu_insert_quad(q);
+
+ // traction at insertion with friction component
+ Vector<Real> new_normal_traction(normal_traction_quad(q));
+ new_normal_traction +=
+ mu_insert_quad(q) * normal_compression_quad(q) * tangent;
+ Vector<Real> new_normal_compression =
+ normal_compression_quad(q) * normal;
+ if (spatial_dimension != 3) {
+ new_normal_traction *= -1.;
+ new_normal_compression *= -1.;
+ }
+
+ if (new_sigma < 0.05 * sigma_lim_quad(q)) {
+ new_mu = mu_insert_avg;
+ new_sigma = effective_stress_quad.mean() +
+ mu_insert_avg / this->beta * avg_normal_compression;
+ new_normal_traction = avg_normal_traction;
+ new_normal_compression = avg_normal_compression * avg_normal;
+ }
+
+ // set delta_c in function of G_c or a given delta_c value
+ Real new_delta_c;
+ if (Math::are_float_equal(this->delta_c, 0.))
+ new_delta_c = 2 * this->G_c / new_sigma;
+ else
+ new_delta_c = (*sigma_lim_it) / new_sigma * this->delta_c;
+
+ // push to insertion list
+ new_sigmas.push_back(new_sigma);
+ new_mus.push_back(new_mu);
+ new_normal_tractions.push_back(new_normal_traction);
+ new_normal_compressions.push_back(new_normal_compression);
+ new_delta_cs.push_back(new_delta_c);
+ }
+ }
+ }
+
+ // 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);
+ mu_eff.resize(old_nb_quad_points + new_nb_quad_points);
+ del_c_eff.resize(old_nb_quad_points + new_nb_quad_points);
+ ins_stress.resize(old_nb_quad_points + new_nb_quad_points);
+ ins_comp.resize(old_nb_quad_points + new_nb_quad_points);
+ trac_old.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];
+ mu_eff(old_nb_quad_points + q) = new_mus[q];
+ del_c_eff(old_nb_quad_points + q) = new_delta_cs[q];
+ for (UInt dim = 0; dim < spatial_dimension; ++dim) {
+ ins_stress(old_nb_quad_points + q, dim) = new_normal_tractions[q](dim);
+ ins_comp(old_nb_quad_points + q, dim) = new_normal_compressions[q](dim);
+ trac_old(old_nb_quad_points + q, dim) = new_normal_tractions[q](dim);
+ }
+ }
+ }
+
+ AKANTU_DEBUG_OUT();
+}
+
+} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc
index 87c4ee528..8da2df437 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc
@@ -1,266 +1,265 @@
/**
* @file material_cohesive_linear_inline_impl.cc
*
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Apr 22 2015
* @date last modification: Wed Feb 21 2018
*
* @brief Inline functions of the MaterialCohesiveLinear
*
* @section LICENSE
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive_linear.hh"
#include "solid_mechanics_model_cohesive.hh"
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__
#define __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template <UInt dim>
inline Real MaterialCohesiveLinear<dim>::computeEffectiveNorm(
const Matrix<Real> & stress, const Vector<Real> & normal,
const Vector<Real> & tangent, Vector<Real> & normal_traction) const {
normal_traction.mul<false>(stress, normal);
Real normal_contrib = normal_traction.dot(normal);
/// in 3D tangential components must be summed
Real tangent_contrib = 0;
if (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<Real> 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 <UInt dim>
inline void MaterialCohesiveLinear<dim>::computeTractionOnQuad(
Vector<Real> & traction, Vector<Real> & opening,
const Vector<Real> & normal, Real & delta_max, const Real & delta_c,
- const Vector<Real> & insertion_stress, const Real & sigma_c,
+ const Vector<Real> & insertion_stress,
+ const Vector<Real> & insertion_compression, const Real & sigma_c,
Vector<Real> & normal_opening, Vector<Real> & tangential_opening,
- Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage,
- bool & penetration, Vector<Real> & contact_traction,
+ Real & normal_opening_norm, Real & tangential_opening_norm,
+ Real & damage, bool & penetration, Vector<Real> & contact_traction,
Vector<Real> & 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 (this->contact_after_breaking == false &&
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.clear();
} else {
delta += normal_opening_norm * normal_opening_norm;
contact_traction.clear();
contact_opening.clear();
}
delta = std::sqrt(delta);
/// update maximum displacement and damage
delta_max = std::max(delta_max, delta);
damage = std::min(delta_max / delta_c, Real(1.));
/**
* Compute traction @f$ \mathbf{T} = \left(
* \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n
* \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1-
* \frac{\delta}{\delta_c} \right)@f$
*/
if (Math::are_float_equal(damage, 1.))
traction.clear();
else if (Math::are_float_equal(damage, 0.)) {
- if (penetration)
- traction.clear();
- else
- traction = insertion_stress;
+ traction = insertion_stress;
+ contact_traction = insertion_compression;
} 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 <UInt dim>
inline void MaterialCohesiveLinear<dim>::computeTangentTractionOnQuad(
Matrix<Real> & tangent, Real & delta_max, const Real & delta_c,
const Real & sigma_c, Vector<Real> & opening, const Vector<Real> & normal,
Vector<Real> & normal_opening, Vector<Real> & tangential_opening,
Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage,
bool & penetration, Vector<Real> & 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 (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<Real> 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<Real> I(spatial_dimension, spatial_dimension);
I.eye(this->beta2_kappa);
Matrix<Real> nn(n_outer_n);
nn *= (1. - this->beta2_kappa);
nn += I;
nn *= t / delta;
Vector<Real> t_tilde(normal_opening);
t_tilde *= (1. - this->beta2_kappa2);
Vector<Real> mm(opening);
mm *= this->beta2_kappa2;
t_tilde += mm;
Vector<Real> t_hat(normal_opening);
t_hat += this->beta2_kappa * tangential_opening;
Matrix<Real> prov(spatial_dimension, spatial_dimension);
prov.outerProduct(t_hat, t_tilde);
prov *= derivative / delta;
prov += nn;
Matrix<Real> prov_t = prov.transpose();
tangent += prov_t;
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
/* -------------------------------------------------------------------------- */
#endif //__AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__
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 1c31d8d0a..9b9a9fc60 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,558 +1,560 @@
/**
* @file material_cohesive.cc
*
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Seyedeh Mohadeseh Taheri Mousavi <mohadeseh.taherimousavi@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Feb 22 2012
* @date last modification: Mon Feb 19 2018
*
* @brief Specialization of the material class for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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<MyFEEngineCohesiveType>("CohesiveFEEngine")),
reversible_energy("reversible_energy", *this),
total_energy("total_energy", *this), opening("opening", *this),
tractions("tractions", *this),
contact_tractions("contact_tractions", *this),
contact_opening("contact_opening", *this), delta_max("delta max", *this),
use_previous_delta_max(false), use_previous_opening(false),
- damage("damage", *this), sigma_c("sigma_c", *this),
- normal(0, spatial_dimension, "normal") {
+ use_previous_contact_opening(false), damage("damage", *this),
+ sigma_c("sigma_c", *this), normal(0, spatial_dimension, "normal") {
AKANTU_DEBUG_IN();
this->model = dynamic_cast<SolidMechanicsModelCohesive *>(&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->opening.initialize(spatial_dimension);
this->opening.initializeHistory();
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();
+ if (this->use_previous_contact_opening)
+ this->contact_opening.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<Array<Real> &>(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$
Array<Real> * traction_cpy = new Array<Real>(
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<Real> 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<Real> & shapes_filtered =
shapes_filtered_begin[current_quad];
traction_tmp.copy(*traction_it);
traction_tmp += *contact_traction_it;
traction_cpy_it->mul<false, false>(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$
*/
Array<Real> * partial_int_t_N = new Array<Real>(
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;
Array<Real> * int_t_N = new Array<Real>(
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<Real> & shapes = fem_cohesive.getShapes(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
UInt nb_element = elem_filter.size();
if (!nb_element)
continue;
UInt size_of_shapes = shapes.getNbComponent();
Array<Real> * shapes_filtered = new Array<Real>(
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<Real> 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$
Array<Real> * tangent_stiffness_matrix = new Array<Real>(
nb_element * nb_quadrature_points,
spatial_dimension * spatial_dimension, "tangent_stiffness_matrix");
// Array<Real> * normal = new Array<Real>(nb_element *
// nb_quadrature_points, spatial_dimension, "normal");
normal.resize(nb_quadrature_points);
computeNormal(model->getCurrentPosition(), normal, type, ghost_type);
/// compute openings @f$\mathbf{\delta}@f$
// computeOpening(model->getDisplacement(), opening(type, ghost_type), type,
// ghost_type);
tangent_stiffness_matrix->clear();
computeTangentTraction(type, *tangent_stiffness_matrix, normal, ghost_type);
// delete normal;
UInt size_at_nt_d_n_a = spatial_dimension * nb_nodes_per_element *
spatial_dimension * nb_nodes_per_element;
Array<Real> * at_nt_d_n_a = new Array<Real>(
nb_element * nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A");
Array<Real>::iterator<Vector<Real>> shapes_filt_it =
shapes_filtered->begin(size_of_shapes);
Array<Real>::matrix_iterator D_it =
tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension);
Array<Real>::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<Real>::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<Real> N(spatial_dimension, spatial_dimension * size_of_shapes);
Matrix<Real> N_A(spatial_dimension,
spatial_dimension * nb_nodes_per_element);
Matrix<Real> 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.clear();
/**
* 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<false, false>(N, A);
D_N_A.mul<false, false>(*D_it, N_A);
(*At_Nt_D_N_A_it).mul<true, false>(D_N_A, N_A);
}
delete tangent_stiffness_matrix;
delete shapes_filtered;
Array<Real> * K_e = new Array<Real>(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 (auto & type : element_filter.elementTypes(spatial_dimension, ghost_type,
_ek_cohesive)) {
Array<UInt> & elem_filter = element_filter(type, ghost_type);
UInt nb_element = elem_filter.size();
if (nb_element == 0)
continue;
UInt nb_quadrature_points =
nb_element * fem_cohesive.getNbIntegrationPoints(type, ghost_type);
normal.resize(nb_quadrature_points);
/// compute normals @f$\mathbf{n}@f$
computeNormal(model->getCurrentPosition(), normal, 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(normal, type, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::computeNormal(const Array<Real> & position,
Array<Real> & normal, ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto & fem_cohesive =
this->model->getFEEngineClass<MyFEEngineCohesiveType>("CohesiveFEEngine");
normal.clear();
#define COMPUTE_NORMAL(type) \
fem_cohesive.getShapeFunctions() \
.computeNormalsOnIntegrationPoints<type, CohesiveReduceFunctionMean>( \
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<Real> & displacement,
Array<Real> & opening, ElementType type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto & fem_cohesive =
this->model->getFEEngineClass<MyFEEngineCohesiveType>("CohesiveFEEngine");
#define COMPUTE_OPENING(type) \
fem_cohesive.getShapeFunctions() \
.interpolateOnIntegrationPoints<type, CohesiveReduceFunctionOpening>( \
displacement, opening, spatial_dimension, ghost_type, \
element_filter(type, ghost_type));
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING);
#undef COMPUTE_OPENING
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MaterialCohesive::updateEnergies(ElementType type) {
AKANTU_DEBUG_IN();
if (Mesh::getKind(type) != _ek_cohesive)
return;
Vector<Real> b(spatial_dimension);
Vector<Real> 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 (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 (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost,
_ek_cohesive)) {
Array<Real> 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 (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<Real> 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) {
AKANTU_DEBUG_IN();
if (type == "reversible")
return getReversibleEnergy();
else if (type == "dissipated")
return getDissipatedEnergy();
else if (type == "cohesive contact")
return getContactEnergy();
AKANTU_DEBUG_OUT();
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 187730ead..239cf0dfe 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,237 +1,242 @@
/**
* @file material_cohesive.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Seyedeh Mohadeseh Taheri Mousavi <mohadeseh.taherimousavi@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Specialization of the material class for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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<IntegratorGauss, ShapeLagrange, _ek_cohesive>;
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<Real> & /*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(const ElementType & /*el_type*/,
Array<Real> & /*tangent_matrix*/,
const Array<Real> & /*normal*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the normal
void computeNormal(const Array<Real> & position, Array<Real> & normal,
ElementType type, GhostType ghost_type);
/// compute the opening
void computeOpening(const Array<Real> & displacement, Array<Real> & opening,
ElementType type, GhostType ghost_type);
template <ElementType type>
void computeNormal(const Array<Real> & position, Array<Real> & normal,
GhostType ghost_type);
/// assemble stiffness
void assembleStiffnessMatrix(GhostType ghost_type) override;
/// constitutive law
virtual void computeTraction(const Array<Real> & normal, ElementType el_type,
GhostType ghost_type = _not_ghost) = 0;
/// parallelism functions
inline UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void unpackData(CommunicationBuffer & buffer,
const Array<Element> & 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 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<UInt> &);
// AKANTU_GET_MACRO(ElementFilter, element_filter, const
// ElementTypeMapArray<UInt> &);
/// 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<UInt> facet_filter;
/// Link to the cohesive fem object in the model
FEEngine & fem_cohesive;
private:
/// reversible energy by quadrature point
CohesiveInternalField<Real> reversible_energy;
/// total energy by quadrature point
CohesiveInternalField<Real> total_energy;
protected:
/// opening in all elements and quadrature points
CohesiveInternalField<Real> opening;
/// traction in all elements and quadrature points
CohesiveInternalField<Real> tractions;
/// traction due to contact
CohesiveInternalField<Real> contact_tractions;
/// normal openings for contact tractions
CohesiveInternalField<Real> contact_opening;
/// maximum displacement
CohesiveInternalField<Real> 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;
+ /// tell if the previous contact_opening state is needed (in iterative schemes)
+ bool use_previous_contact_opening;
+
+ /// tell if the
+
/// damage
CohesiveInternalField<Real> damage;
/// pointer to the solid mechanics model for cohesive elements
SolidMechanicsModelCohesive * model;
/// critical stress
RandomInternalField<Real, FacetInternalField> sigma_c;
/// critical displacement
Real delta_c;
/// array to temporarily store the normals
Array<Real> normal;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_cohesive_inline_impl.cc"
} // namespace akantu
#include "cohesive_internal_field_tmpl.hh"
#endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_includes.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_includes.hh
index 6bb8476c1..c16b4fd48 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_includes.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_includes.hh
@@ -1,49 +1,52 @@
/**
* @file material_cohesive_includes.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Sun Sep 26 2010
* @date last modification: Fri Oct 13 2017
*
* @brief List of includes for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
// /* --------------------------------------------------------------------------
// */
// #ifndef AKANTU_CMAKE_LIST_MATERIALS
// #include "material_cohesive.hh"
// #include "material_cohesive_bilinear.hh"
// #include "material_cohesive_exponential.hh"
// #include "material_cohesive_linear.hh"
// #include "material_cohesive_linear_fatigue.hh"
// #include "material_cohesive_linear_friction.hh"
// #include "material_cohesive_linear_uncoupled.hh"
// #endif
#define AKANTU_COHESIVE_MATERIAL_LIST \
((2, (cohesive_linear, MaterialCohesiveLinear)))( \
(2, (cohesive_linear_fatigue, MaterialCohesiveLinearFatigue)))( \
- (2, (cohesive_linear_friction, MaterialCohesiveLinearFriction)))( \
+ (2, (cohesive_linear_friction_coulomb, \
+ MaterialCohesiveLinearFrictionCoulomb)))( \
+ (2, (cohesive_linear_friction_slipweakening, \
+ MaterialCohesiveLinearFrictionSlipWeakening)))( \
(2, (cohesive_linear_uncoupled, MaterialCohesiveLinearUncoupled)))( \
(2, (cohesive_bilinear, MaterialCohesiveBilinear)))( \
(2, (cohesive_exponential, MaterialCohesiveExponential)))
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc
index 8a549afc0..9b7577f9b 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc
@@ -1,708 +1,708 @@
/**
* @file solid_mechanics_model_cohesive.cc
*
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue May 08 2012
* @date last modification: Wed Feb 21 2018
*
* @brief Solid mechanics model for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model_cohesive.hh"
#include "aka_iterators.hh"
#include "cohesive_element_inserter.hh"
#include "element_synchronizer.hh"
#include "facet_synchronizer.hh"
#include "fe_engine_template.hh"
#include "global_ids_updater.hh"
#include "integrator_gauss.hh"
#include "material_cohesive.hh"
#include "mesh_accessor.hh"
#include "mesh_global_data_updater.hh"
#include "parser.hh"
#include "shape_cohesive.hh"
/* -------------------------------------------------------------------------- */
#include "dumpable_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
/* -------------------------------------------------------------------------- */
#include <algorithm>
/* -------------------------------------------------------------------------- */
namespace akantu {
class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater {
public:
CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model)
: model(model), mesh(model.getMesh()),
global_ids_updater(model.getMesh(), *model.cohesive_synchronizer) {}
/* ------------------------------------------------------------------------ */
std::tuple<UInt, UInt>
updateData(NewNodesEvent & nodes_event,
NewElementsEvent & elements_event) override {
auto cohesive_nodes_event =
dynamic_cast<CohesiveNewNodesEvent *>(&nodes_event);
if (not cohesive_nodes_event)
return std::make_tuple(nodes_event.getList().size(),
elements_event.getList().size());
/// update nodes type
auto & new_nodes = cohesive_nodes_event->getList();
auto & old_nodes = cohesive_nodes_event->getOldNodesList();
auto local_nb_new_nodes = new_nodes.size();
auto nb_new_nodes = local_nb_new_nodes;
if (mesh.isDistributed()) {
MeshAccessor mesh_accessor(mesh);
auto & nodes_flags = mesh_accessor.getNodesFlags();
auto nb_old_nodes = nodes_flags.size();
nodes_flags.resize(nb_old_nodes + local_nb_new_nodes);
for (auto && data : zip(old_nodes, new_nodes)) {
UInt old_node, new_node;
std::tie(old_node, new_node) = data;
nodes_flags(new_node) = nodes_flags(old_node);
}
model.updateCohesiveSynchronizers();
nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size());
}
Vector<UInt> nb_new_stuff = {nb_new_nodes, elements_event.getList().size()};
const auto & comm = mesh.getCommunicator();
comm.allReduce(nb_new_stuff, SynchronizerOperation::_sum);
if (nb_new_stuff(1) > 0) {
mesh.sendEvent(elements_event);
MeshUtils::resetFacetToDouble(mesh.getMeshFacets());
}
if (nb_new_stuff(0) > 0) {
mesh.sendEvent(nodes_event);
// mesh.sendEvent(global_ids_updater.getChangedNodeEvent());
}
return std::make_tuple(nb_new_stuff(0), nb_new_stuff(1));
}
private:
SolidMechanicsModelCohesive & model;
Mesh & mesh;
GlobalIdsUpdater global_ids_updater;
};
/* -------------------------------------------------------------------------- */
SolidMechanicsModelCohesive::SolidMechanicsModelCohesive(
Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id)
: SolidMechanicsModel(mesh, dim, id, memory_id,
ModelType::_solid_mechanics_model_cohesive),
tangents("tangents", id), facet_stress("facet_stress", id),
facet_material("facet_material", id) {
AKANTU_DEBUG_IN();
registerFEEngineObject<MyFEEngineCohesiveType>("CohesiveFEEngine", mesh,
Model::spatial_dimension);
auto && tmp_material_selector =
std::make_shared<DefaultMaterialCohesiveSelector>(*this);
tmp_material_selector->setFallback(this->material_selector);
this->material_selector = tmp_material_selector;
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("cohesive elements", id);
this->mesh.addDumpMeshToDumper("cohesive elements", mesh,
Model::spatial_dimension, _not_ghost,
_ek_cohesive);
#endif
if (this->mesh.isDistributed()) {
/// create the distributed synchronizer for cohesive elements
this->cohesive_synchronizer = std::make_unique<ElementSynchronizer>(
mesh, "cohesive_distributed_synchronizer");
auto & synchronizer = mesh.getElementSynchronizer();
this->cohesive_synchronizer->split(synchronizer, [](auto && el) {
return Mesh::getKind(el.type) == _ek_cohesive;
});
this->registerSynchronizer(*cohesive_synchronizer,
SynchronizationTag::_material_id);
this->registerSynchronizer(*cohesive_synchronizer,
SynchronizationTag::_smm_stress);
this->registerSynchronizer(*cohesive_synchronizer,
SynchronizationTag::_smm_boundary);
}
this->inserter = std::make_unique<CohesiveElementInserter>(
this->mesh, id + ":cohesive_element_inserter");
registerFEEngineObject<MyFEEngineFacetType>(
"FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default;
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::setTimeStep(Real time_step,
const ID & solver_id) {
SolidMechanicsModel::setTimeStep(time_step, solver_id);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper("cohesive elements").setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) {
AKANTU_DEBUG_IN();
const auto & smmc_options =
aka::as_type<SolidMechanicsModelCohesiveOptions>(options);
this->is_extrinsic = smmc_options.is_extrinsic;
inserter->setIsExtrinsic(is_extrinsic);
if (mesh.isDistributed()) {
auto & mesh_facets = inserter->getMeshFacets();
auto & synchronizer =
aka::as_type<FacetSynchronizer>(mesh_facets.getElementSynchronizer());
// synchronizeGhostFacetsConnectivity();
/// create the facet synchronizer for extrinsic simulations
if (is_extrinsic) {
facet_stress_synchronizer = std::make_unique<ElementSynchronizer>(
synchronizer, id + ":facet_stress_synchronizer");
facet_stress_synchronizer->swapSendRecv();
this->registerSynchronizer(*facet_stress_synchronizer,
SynchronizationTag::_smmc_facets_stress);
}
}
MeshAccessor mesh_accessor(mesh);
mesh_accessor.registerGlobalDataUpdater(
std::make_unique<CohesiveMeshGlobalDataUpdater>(*this));
ParserSection section;
bool is_empty;
std::tie(section, is_empty) = this->getParserSection();
if (not is_empty) {
auto inserter_section =
section.getSubSections(ParserType::_cohesive_inserter);
if (inserter_section.begin() != inserter_section.end()) {
inserter->parseSection(*inserter_section.begin());
}
}
SolidMechanicsModel::initFullImpl(options);
AKANTU_DEBUG_OUT();
} // namespace akantu
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initMaterials() {
AKANTU_DEBUG_IN();
// make sure the material are instantiated
if (not are_materials_instantiated)
instantiateMaterials();
/// find the first cohesive material
UInt cohesive_index = UInt(-1);
for (auto && material : enumerate(materials)) {
if (dynamic_cast<MaterialCohesive *>(std::get<1>(material).get())) {
cohesive_index = std::get<0>(material);
break;
}
}
if (cohesive_index == UInt(-1))
AKANTU_EXCEPTION("No cohesive materials in the material input file");
material_selector->setFallback(cohesive_index);
// set the facet information in the material in case of dynamic insertion
// to know what material to call for stress checks
const Mesh & mesh_facets = inserter->getMeshFacets();
facet_material.initialize(
mesh_facets, _spatial_dimension = spatial_dimension - 1,
_with_nb_element = true,
_default_value = material_selector->getFallbackValue());
for_each_element(
mesh_facets,
[&](auto && element) {
auto mat_index = (*material_selector)(element);
auto & mat = aka::as_type<MaterialCohesive>(*materials[mat_index]);
facet_material(element) = mat_index;
if (is_extrinsic) {
mat.addFacet(element);
}
},
_spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost);
SolidMechanicsModel::initMaterials();
if (is_extrinsic) {
this->initAutomaticInsertion();
} else {
this->insertIntrinsicElements();
}
AKANTU_DEBUG_OUT();
} // namespace akantu
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*/
void SolidMechanicsModelCohesive::initModel() {
AKANTU_DEBUG_IN();
SolidMechanicsModel::initModel();
/// add cohesive type connectivity
ElementType type = _not_defined;
for (auto && type_ghost : ghost_types) {
for (const auto & tmp_type :
mesh.elementTypes(spatial_dimension, type_ghost)) {
const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost);
if (connectivity.size() == 0)
continue;
type = tmp_type;
auto type_facet = Mesh::getFacetType(type);
auto type_cohesive = FEEngine::getCohesiveElementType(type_facet);
mesh.addConnectivityType(type_cohesive, type_ghost);
}
}
AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh");
getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost);
getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost);
if (is_extrinsic) {
getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost);
getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::insertIntrinsicElements() {
AKANTU_DEBUG_IN();
inserter->insertIntrinsicElements();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initAutomaticInsertion() {
AKANTU_DEBUG_IN();
this->inserter->limitCheckFacets();
this->updateFacetStressSynchronizer();
this->resizeFacetStress();
/// compute normals on facets
this->computeNormals();
this->initStressInterpolation();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::updateAutomaticInsertion() {
AKANTU_DEBUG_IN();
this->inserter->limitCheckFacets();
this->updateFacetStressSynchronizer();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initStressInterpolation() {
Mesh & mesh_facets = inserter->getMeshFacets();
/// compute quadrature points coordinates on facets
Array<Real> & position = mesh.getNodes();
ElementTypeMapArray<Real> quad_facets("quad_facets", id);
quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension,
_spatial_dimension = Model::spatial_dimension - 1);
// mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension,
// Model::spatial_dimension - 1);
getFEEngine("FacetsFEEngine")
.interpolateOnIntegrationPoints(position, quad_facets);
/// compute elements quadrature point positions and build
/// element-facet quadrature points data structure
ElementTypeMapArray<Real> elements_quad_facets("elements_quad_facets", id);
elements_quad_facets.initialize(
mesh, _nb_component = Model::spatial_dimension,
_spatial_dimension = Model::spatial_dimension);
// mesh.initElementTypeMapArray(elements_quad_facets,
// Model::spatial_dimension,
// Model::spatial_dimension);
for (auto elem_gt : ghost_types) {
for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) {
UInt nb_element = mesh.getNbElement(type, elem_gt);
if (nb_element == 0)
continue;
/// compute elements' quadrature points and list of facet
/// quadrature points positions by element
const auto & facet_to_element =
mesh_facets.getSubelementToElement(type, elem_gt);
auto & el_q_facet = elements_quad_facets(type, elem_gt);
auto facet_type = Mesh::getFacetType(type);
auto nb_quad_per_facet =
getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type);
auto nb_facet_per_elem = facet_to_element.getNbComponent();
// small hack in the loop to skip boundary elements, they are silently
// initialized to NaN to see if this causes problems
el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet,
std::numeric_limits<Real>::quiet_NaN());
for (auto && data :
zip(make_view(facet_to_element),
make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) {
const auto & global_facet = std::get<0>(data);
auto & el_q = std::get<1>(data);
if (global_facet == ElementNull)
continue;
Matrix<Real> quad_f =
make_view(quad_facets(global_facet.type, global_facet.ghost_type),
spatial_dimension, nb_quad_per_facet)
.begin()[global_facet.element];
el_q = quad_f;
// for (UInt q = 0; q < nb_quad_per_facet; ++q) {
// for (UInt s = 0; s < Model::spatial_dimension; ++s) {
// el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet +
// f * nb_quad_per_facet + q,
// s) = quad_f(global_facet * nb_quad_per_facet + q,
// s);
// }
// }
//}
}
}
}
/// loop over non cohesive materials
for (auto && material : materials) {
if (dynamic_cast<MaterialCohesive *>(material.get()))
continue;
/// initialize the interpolation function
material->initElementalFieldInterpolation(elements_quad_facets);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::assembleInternalForces() {
AKANTU_DEBUG_IN();
// f_int += f_int_cohe
for (auto & material : this->materials) {
try {
auto & mat = aka::as_type<MaterialCohesive>(*material);
mat.computeTraction(_not_ghost);
} catch (std::bad_cast & bce) {
}
}
SolidMechanicsModel::assembleInternalForces();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::computeNormals() {
AKANTU_DEBUG_IN();
Mesh & mesh_facets = this->inserter->getMeshFacets();
this->getFEEngine("FacetsFEEngine")
.computeNormalsOnIntegrationPoints(_not_ghost);
/**
* @todo store tangents while computing normals instead of
* recomputing them as follows:
*/
/* ------------------------------------------------------------------------ */
UInt tangent_components =
Model::spatial_dimension * (Model::spatial_dimension - 1);
tangents.initialize(mesh_facets, _nb_component = tangent_components,
_spatial_dimension = Model::spatial_dimension - 1);
// mesh_facets.initElementTypeMapArray(tangents, tangent_components,
// Model::spatial_dimension - 1);
for (auto facet_type :
mesh_facets.elementTypes(Model::spatial_dimension - 1)) {
const Array<Real> & normals =
this->getFEEngine("FacetsFEEngine")
.getNormalsOnIntegrationPoints(facet_type);
Array<Real> & tangents = this->tangents(facet_type);
Math::compute_tangents(normals, tangents);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::interpolateStress() {
ElementTypeMapArray<Real> by_elem_result("temporary_stress_by_facets", id);
for (auto & material : materials) {
auto * mat = dynamic_cast<MaterialCohesive *>(material.get());
if (mat == nullptr)
/// interpolate stress on facet quadrature points positions
material->interpolateStressOnFacets(facet_stress, by_elem_result);
}
this->synchronize(SynchronizationTag::_smmc_facets_stress);
}
/* -------------------------------------------------------------------------- */
UInt SolidMechanicsModelCohesive::checkCohesiveStress() {
AKANTU_DEBUG_IN();
if (not is_extrinsic) {
AKANTU_EXCEPTION(
"This function can only be used for extrinsic cohesive elements");
}
interpolateStress();
for (auto & mat : materials) {
auto * mat_cohesive = dynamic_cast<MaterialCohesive *>(mat.get());
if (mat_cohesive) {
/// check which not ghost cohesive elements are to be created
mat_cohesive->checkInsertion();
}
}
/// communicate data among processors
// this->synchronize(SynchronizationTag::_smmc_facets);
/// insert cohesive elements
UInt nb_new_elements = inserter->insertElements();
// if (nb_new_elements > 0) {
// this->reinitializeSolver();
// }
AKANTU_DEBUG_OUT();
return nb_new_elements;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::onElementsAdded(
const Array<Element> & element_list, const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
SolidMechanicsModel::onElementsAdded(element_list, event);
if (is_extrinsic)
resizeFacetStress();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::onNodesAdded(const Array<UInt> & new_nodes,
const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
SolidMechanicsModel::onNodesAdded(new_nodes, event);
const CohesiveNewNodesEvent * cohesive_event;
if ((cohesive_event = dynamic_cast<const CohesiveNewNodesEvent *>(&event)) ==
nullptr)
return;
const auto & old_nodes = cohesive_event->getOldNodesList();
auto copy = [this, &new_nodes, &old_nodes](auto & arr) {
UInt new_node, old_node;
auto view = make_view(arr, spatial_dimension);
auto begin = view.begin();
for (auto && pair : zip(new_nodes, old_nodes)) {
std::tie(new_node, old_node) = pair;
auto old_ = begin + old_node;
auto new_ = begin + new_node;
*new_ = *old_;
}
};
-
copy(*displacement);
+
copy(*blocked_dofs);
if (velocity)
copy(*velocity);
if (acceleration)
copy(*acceleration);
if (current_position)
copy(*current_position);
if (previous_displacement)
copy(*previous_displacement);
// if (external_force)
// copy(*external_force);
// if (internal_force)
// copy(*internal_force);
if (displacement_increment)
copy(*displacement_increment);
copy(getDOFManager().getSolution("displacement"));
// this->assembleMassLumped();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::afterSolveStep(bool converged) {
AKANTU_DEBUG_IN();
/*
* This is required because the Cauchy stress is the stress measure that
* is used to check the insertion of cohesive elements
*/
if (converged) {
for (auto & mat : materials) {
if (mat->isFiniteDeformation())
mat->computeAllCauchyStresses(_not_ghost);
}
}
SolidMechanicsModel::afterSolveStep(converged);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::printself(std::ostream & stream,
int indent) const {
std::string space(indent, AKANTU_INDENT);
stream << space << "SolidMechanicsModelCohesive ["
<< "\n";
SolidMechanicsModel::printself(stream, indent + 2);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::resizeFacetStress() {
AKANTU_DEBUG_IN();
this->facet_stress.initialize(getFEEngine("FacetsFEEngine"),
_nb_component =
2 * spatial_dimension * spatial_dimension,
_spatial_dimension = spatial_dimension - 1);
// for (auto && ghost_type : ghost_types) {
// for (const auto & type :
// mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) {
// UInt nb_facet = mesh_facets.getNbElement(type, ghost_type);
// UInt nb_quadrature_points = getFEEngine("FacetsFEEngine")
// .getNbIntegrationPoints(type,
// ghost_type);
// UInt new_size = nb_facet * nb_quadrature_points;
// facet_stress(type, ghost_type).resize(new_size);
// }
// }
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper(
const std::string & dumper_name, const std::string & field_id,
const std::string & group_name, const ElementKind & element_kind,
bool padding_flag) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = Model::spatial_dimension;
ElementKind _element_kind = element_kind;
if (dumper_name == "cohesive elements") {
_element_kind = _ek_cohesive;
} else if (dumper_name == "facets") {
spatial_dimension = Model::spatial_dimension - 1;
}
SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id,
group_name, spatial_dimension,
_element_kind, padding_flag);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::onDump() {
this->flattenAllRegisteredInternals(_ek_cohesive);
SolidMechanicsModel::onDump();
}
/* -------------------------------------------------------------------------- */
} // namespace akantu