Revision 3f51786289e66cdebc74007458f43157b17ca0b8 authored by Dmitri Naumov on 01 October 2021, 11:32:23 UTC, committed by Dmitri Naumov on 01 October 2021, 18:38:49 UTC
Update the ctest config for updated time step output. The reference files are still the same, no changes required.
1 parent f63a01b
CohesiveZoneModeI.cpp
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#include "CohesiveZoneModeI.h"
#include "LogPenalty.h"
#include "MathLib/MathTools.h"
namespace MaterialLib
{
namespace Fracture
{
namespace CohesiveZoneModeI
{
namespace
{
double computeDamage(double const damage_prev,
double const w_n,
double const w_np,
double const w_nf)
{
return std::min(
1.0,
std::max(damage_prev, std::max(0.0, (w_n - w_np)) / (w_nf - w_np)));
}
} // namespace
template <int DisplacementDim>
void CohesiveZoneModeI<DisplacementDim>::computeConstitutiveRelation(
double const t,
ParameterLib::SpatialPosition const& x,
double const aperture0,
Eigen::Ref<Eigen::VectorXd const>
sigma0,
Eigen::Ref<Eigen::VectorXd const>
/*w_prev*/,
Eigen::Ref<Eigen::VectorXd const>
w,
Eigen::Ref<Eigen::VectorXd const>
/*sigma_prev*/,
Eigen::Ref<Eigen::VectorXd>
sigma,
Eigen::Ref<Eigen::MatrixXd>
C,
typename FractureModelBase<DisplacementDim>::MaterialStateVariables&
material_state_variables)
{
assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
&material_state_variables) != nullptr);
auto& state =
static_cast<StateVariables<DisplacementDim>&>(material_state_variables);
// reset damage in each iteration
state.setInitialConditions();
auto const mp = evaluatedMaterialProperties(t, x);
C.setZero();
// Separately compute shear and normal stresses because of the penalty for
// the normal component.
const int index_ns = DisplacementDim - 1;
double const w_n = w[index_ns];
for (int i = 0; i < index_ns; i++)
{
C(i, i) = mp.Ks;
}
sigma.noalias() = C * w;
double const aperture = w_n + aperture0;
sigma.coeffRef(index_ns) =
mp.Kn * w_n * logPenalty(aperture0, aperture, _penalty_aperture_cutoff);
C(index_ns, index_ns) =
mp.Kn *
logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
sigma.noalias() += sigma0;
// Exit if fracture is closing
if (sigma[index_ns] < 0)
{
return; /// Undamaged stiffness used in compression.
}
//
// Continue with fracture opening.
//
// effective opening calculated by shift to state corresponding to sigma0=0
const double w_n_effective = w_n + sigma0[index_ns] / mp.Kn;
double degradation = 0.0;
if (mp.w_nf == 0.0) // can be used to model tension-cutoff
{
state.damage_prev = 1.0;
state.damage = 1.0;
}
else
{
state.damage =
computeDamage(state.damage_prev, w_n_effective, mp.w_np, mp.w_nf);
degradation = ((1 - state.damage) * mp.w_np) /
(mp.w_np + state.damage * (mp.w_nf - mp.w_np));
}
// Degrade stiffness tensor in tension.
C = C * degradation;
// Create effective relative displacement vector
Eigen::VectorXd w_eff = w;
// Only origin of normal entry is shifted, shear unaffected
w_eff[index_ns] = w_n_effective;
sigma = C * w_eff;
if (state.damage > state.damage_prev)
{
// If damage is increasing, provide extension to consistent tangent.
Eigen::Matrix<double, DisplacementDim, 1> dd_dw =
Eigen::Matrix<double, DisplacementDim, 1>::Zero();
double const tmp = mp.w_np + state.damage * (mp.w_nf - mp.w_np);
dd_dw[index_ns] =
(mp.w_np * mp.w_nf) / ((mp.w_nf - mp.w_np) * (tmp * tmp));
C -= C * w_eff * (dd_dw).transpose();
}
// TODO (nagel) Initial stress not considered, yet.
// sigma.noalias() += sigma0;
}
template class CohesiveZoneModeI<2>;
template class CohesiveZoneModeI<3>;
} // namespace CohesiveZoneModeI
} // namespace Fracture
} // namespace MaterialLib
Computing file changes ...