From b61bdb31d0bf3801fc391d8d13ce9bbc7d6210cf Mon Sep 17 00:00:00 2001
From: "Dmitry Yu. Naumov" <github@naumov.de>
Date: Fri, 9 Nov 2018 18:47:11 +0100
Subject: [PATCH] [MatL] Extract calculateDegradedStress().
The multiple inheritance for the phase field extension
is not actually needed and can be completely avoided.
---
.../LinearElasticIsotropicPhaseField.cpp | 21 ---
.../LinearElasticIsotropicPhaseField.h | 175 +++++++-----------
MaterialLib/SolidModels/PhaseFieldExtension.h | 89 ---------
ProcessLib/PhaseField/PhaseFieldFEM.h | 24 ++-
4 files changed, 78 insertions(+), 231 deletions(-)
delete mode 100644 MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp
delete mode 100644 MaterialLib/SolidModels/PhaseFieldExtension.h
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp
deleted file mode 100644
index 2cfed161f5d..00000000000
--- a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp
+++ /dev/null
@@ -1,21 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, 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 "LinearElasticIsotropicPhaseField.h"
-
-namespace MaterialLib
-{
-namespace Solids
-{
-
-template class LinearElasticIsotropicPhaseField<2>;
-template class LinearElasticIsotropicPhaseField<3>;
-
-} // namespace Solids
-} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
index 003c68e70ff..ae497050da9 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
@@ -11,135 +11,84 @@
#include "MathLib/KelvinVector.h"
-#include "LinearElasticIsotropic.h"
-#include "PhaseFieldExtension.h"
-
namespace MaterialLib
{
namespace Solids
{
+namespace Phasefield
+{
+/** Decompose the stiffness into tensile and compressive part.
+ * Judging by the physical observations, compression perpendicular
+ * to a crack does not cause crack propagation. Thus,
+ * the phase-field parameter is only involved into the tensile part
+ * to degrade the elastic strain energy.
+ */
template <int DisplacementDim>
-class LinearElasticIsotropicPhaseField final
- : public LinearElasticIsotropic<DisplacementDim>,
- public PhaseFieldExtension<DisplacementDim>
+std::tuple<MathLib::KelvinVector::KelvinVectorType<
+ DisplacementDim> /* sigma_real */,
+ MathLib::KelvinVector::KelvinVectorType<
+ DisplacementDim> /* sigma_tensile */,
+ MathLib::KelvinVector::KelvinMatrixType<
+ DisplacementDim> /* C_tensile */,
+ MathLib::KelvinVector::KelvinMatrixType<
+ DisplacementDim> /* C_compressive */,
+ double /* strain_energy_tensile */,
+ double /* elastic_energy */
+ >
+calculateDegradedStress(
+ double const degradation,
+ double const bulk_modulus,
+ double const mu,
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
{
-public:
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ // calculation of deviatoric parts
+ auto const& P_dev = Invariants::deviatoric_projection;
+ KelvinVector const epsd_curr = P_dev * eps;
- explicit LinearElasticIsotropicPhaseField(
- typename LinearElasticIsotropic<DisplacementDim>::MaterialProperties&&
- material_properties)
- : LinearElasticIsotropic<DisplacementDim>(
- std::move(material_properties))
- {
- }
-
- std::unique_ptr<
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
- createMaterialStateVariables() const override
- {
- return LinearElasticIsotropic<
- DisplacementDim>::createMaterialStateVariables();
- }
-
- boost::optional<std::tuple<KelvinVector,
- std::unique_ptr<typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables>,
- KelvinMatrix>>
- integrateStress(
- double const t, ProcessLib::SpatialPosition const& x, double const dt,
- KelvinVector const& eps_prev, KelvinVector const& eps,
- KelvinVector const& sigma_prev,
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
- material_state_variables,
- double const T) const override
- {
- return LinearElasticIsotropic<DisplacementDim>::integrateStress(
- t, x, dt, eps_prev, eps, sigma_prev, material_state_variables, T);
- }
+ // Hydrostatic part for the stress and the tangent.
+ double const eps_curr_trace = Invariants::trace(eps);
- double computeFreeEnergyDensity(
- double const t,
- ProcessLib::SpatialPosition const& x,
- double const dt,
- KelvinVector const& eps,
- KelvinVector const& sigma,
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
- material_state_variables) const override
- {
- return LinearElasticIsotropic<DisplacementDim>::
- computeFreeEnergyDensity(
- t, x, dt, eps, sigma, material_state_variables);
- }
+ KelvinMatrix C_tensile = KelvinMatrix::Zero();
+ KelvinMatrix C_compressive = KelvinMatrix::Zero();
- /** Decompose the stiffness into tensile and compressive part.
- * Judging by the physical observations, compression perpendicular
- * to a crack does not cause crack propagation. Thus,
- * the phase-field parameter is only involved into the tensile part
- * to degrade the elastic strain energy.
- */
- bool calculateDegradedStress(double const t,
- ProcessLib::SpatialPosition const& x,
- KelvinVector const& eps,
- double& strain_energy_tensile,
- KelvinVector& sigma_tensile,
- KelvinVector& sigma_compressive,
- KelvinMatrix& C_tensile,
- KelvinMatrix& C_compressive,
- KelvinVector& sigma_real,
- double const degradation,
- double& elastic_energy) const override
+ if (eps_curr_trace >= 0)
{
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
- // calculation of deviatoric parts
- auto const& P_dev = Invariants::deviatoric_projection;
- KelvinVector const epsd_curr = P_dev * eps;
-
- // Hydrostatic part for the stress and the tangent.
- double const eps_curr_trace = Invariants::trace(eps);
-
- auto const& K =
- LinearElasticIsotropic<DisplacementDim>::_mp.bulk_modulus(t, x);
- auto const& mu = LinearElasticIsotropic<DisplacementDim>::_mp.mu(t, x);
-
- C_tensile = KelvinMatrix::Zero();
- C_compressive = KelvinMatrix::Zero();
-
- if (eps_curr_trace >= 0)
- {
- strain_energy_tensile = K / 2 * eps_curr_trace * eps_curr_trace +
- mu * epsd_curr.transpose() * epsd_curr;
- sigma_tensile.noalias() =
- K * eps_curr_trace * Invariants::identity2 + 2 * mu * epsd_curr;
- sigma_compressive.noalias() = KelvinVector::Zero();
- C_tensile.template topLeftCorner<3, 3>().setConstant(K);
- C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
- elastic_energy = degradation * strain_energy_tensile;
- }
- else
- {
- strain_energy_tensile = mu * epsd_curr.transpose() * epsd_curr;
- sigma_tensile.noalias() = 2 * mu * epsd_curr;
- sigma_compressive.noalias() =
- K * eps_curr_trace * Invariants::identity2;
- C_tensile.noalias() = 2 * mu * P_dev * KelvinMatrix::Identity();
- C_compressive.template topLeftCorner<3, 3>().setConstant(K);
- elastic_energy = K / 2 * eps_curr_trace * eps_curr_trace +
- mu * epsd_curr.transpose() * epsd_curr;
- }
-
- sigma_real.noalias() = degradation * sigma_tensile + sigma_compressive;
- return true;
+ double const strain_energy_tensile =
+ bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
+ mu * epsd_curr.transpose() * epsd_curr;
+ KelvinVector const sigma_tensile =
+ bulk_modulus * eps_curr_trace * Invariants::identity2 +
+ 2 * mu * epsd_curr;
+ C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
+ C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
+ double const elastic_energy = degradation * strain_energy_tensile;
+ KelvinVector const sigma_real =
+ degradation * sigma_tensile; // + sigma_compressive, which is zero;
+ return std::make_tuple(sigma_real, sigma_tensile, C_tensile,
+ C_compressive, strain_energy_tensile,
+ elastic_energy);
}
-};
-
-extern template class LinearElasticIsotropicPhaseField<2>;
-extern template class LinearElasticIsotropicPhaseField<3>;
-
+ double const strain_energy_tensile = mu * epsd_curr.transpose() * epsd_curr;
+ KelvinVector const sigma_tensile = 2 * mu * epsd_curr;
+ KelvinVector const sigma_compressive =
+ bulk_modulus * eps_curr_trace * Invariants::identity2;
+ C_tensile.noalias() = 2 * mu * P_dev * KelvinMatrix::Identity();
+ C_compressive.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
+ double const elastic_energy =
+ bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
+ mu * epsd_curr.transpose() * epsd_curr;
+ KelvinVector const sigma_real =
+ degradation * sigma_tensile + sigma_compressive;
+ return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
+ strain_energy_tensile, elastic_energy);
+}
+} // namespace Phasefield
} // namespace Solids
} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/PhaseFieldExtension.h b/MaterialLib/SolidModels/PhaseFieldExtension.h
deleted file mode 100644
index 425da16cf2b..00000000000
--- a/MaterialLib/SolidModels/PhaseFieldExtension.h
+++ /dev/null
@@ -1,89 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include "MechanicsBase.h"
-
-namespace MaterialLib
-{
-namespace Solids
-{
-template <int DisplacementDim>
-struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
-{
- using KelvinVector =
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
- using KelvinMatrix =
- MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
- virtual bool calculateDegradedStress(double const t,
- ProcessLib::SpatialPosition const& x,
- KelvinVector const& eps,
- double& strain_energy_tensile,
- KelvinVector& sigma_tensile,
- KelvinVector& sigma_compressive,
- KelvinMatrix& C_tensile,
- KelvinMatrix& C_compressive,
- KelvinVector& sigma_real,
- double const degradation,
- double& elastic_energy) const = 0;
-
- /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
- /// constitutive relation compute function.
- bool calculateDegradedStress(
- double const t,
- ProcessLib::SpatialPosition const& x,
- Eigen::Matrix<double, Eigen::Dynamic, 1> const& eps,
- double& strain_energy_tensile,
- Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_tensile,
- Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_compressive,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- C_tensile,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- C_compressive,
- Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_real,
- double const degradation,
- double& elastic_energy) const
- {
- // TODO Avoid copies of data:
- // Using MatrixBase<Derived> not possible because template functions
- // cannot be virtual. Maybe there is a workaround for this. Using
- // Map<Matrix<double, ...>> makes the interface (for the material model
- // implementation) unnecessary difficult.
- KelvinVector const eps_{eps};
- KelvinVector sigma_tensile_{sigma_tensile};
- KelvinVector sigma_compressive_{sigma_compressive};
- KelvinMatrix C_tensile_{C_tensile};
- KelvinMatrix C_compressive_{C_compressive};
- KelvinVector sigma_real_{sigma_real};
-
- bool const result = calculateDegradedStress(t,
- x,
- eps_,
- strain_energy_tensile,
- sigma_tensile_,
- sigma_compressive_,
- C_tensile_,
- C_compressive_,
- sigma_real_,
- degradation,
- elastic_energy);
-
- sigma_tensile = sigma_tensile_;
- sigma_compressive = sigma_compressive_;
- C_tensile = C_tensile_;
- C_compressive = C_compressive_;
- sigma_real = sigma_real_;
- return result;
- }
-};
-
-} // namespace Solids
-} // namespace MaterialLib
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index ec85389e166..32466594ff3 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -70,20 +70,28 @@ struct IntegrationPointData final
template <typename DisplacementVectorType>
void updateConstitutiveRelation(double const t,
- SpatialPosition const& x_position,
+ SpatialPosition const& x,
double const /*dt*/,
DisplacementVectorType const& /*u*/,
double const degradation)
{
- static_cast<
- MaterialLib::Solids::PhaseFieldExtension<DisplacementDim> const&>(
- solid_material)
- .calculateDegradedStress(t, x_position, eps, strain_energy_tensile,
- sigma_tensile, sigma_compressive,
- C_tensile, C_compressive, sigma,
- degradation, elastic_energy);
+ auto linear_elastic_mp =
+ static_cast<MaterialLib::Solids::LinearElasticIsotropic<
+ DisplacementDim> const&>(solid_material)
+ .getMaterialProperties();
+
+ auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
+ auto const mu = linear_elastic_mp.mu(t, x);
+
+ std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
+ strain_energy_tensile, elastic_energy) =
+ MaterialLib::Solids::Phasefield::calculateDegradedStress<
+ DisplacementDim>(degradation, bulk_modulus, mu, eps);
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+
+ static constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
};
/// Used by for extrapolation of the integration point values. It is ordered
--
GitLab