diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.h b/MaterialLib/SolidModels/LinearElasticIsotropic.h
index 2bd88582428d34a70668280d0816076ccb72ea17..5e9271dc0192cb4a000497de952495ebbd183ad6 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.h
@@ -124,8 +124,8 @@ public:
ProcessLib::SpatialPosition const& x,
double const T) const;
+ MaterialProperties getMaterialProperties() const { return _mp; }
- MaterialProperties getMaterialProperties() { return _mp; }
protected:
MaterialProperties _mp;
};
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp
deleted file mode 100644
index 2cfed161f5dfe6b4d610b81eb6a3b19358afa4a8..0000000000000000000000000000000000000000
--- 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 003c68e70ffafcbb3a7413cf94d727e95ad5e825..ae497050da9dd8e685e08df0381f82d4ef214ee2 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 425da16cf2b6ad27b77d878dd100da01705007d6..0000000000000000000000000000000000000000
--- 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/CreatePhaseFieldProcess.cpp b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
index 1d0195881275ddbd8724db66b7b14b7e06906249..63dad638ad686b30a365c4ba8e0fac62bfbd3846 100644
--- a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
@@ -11,9 +11,10 @@
#include <cassert>
-#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
-#include "MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
#include "PhaseFieldProcess.h"
#include "PhaseFieldProcessData.h"
@@ -99,40 +100,14 @@ std::unique_ptr<Process> createPhaseFieldProcess(
variable_ph->getNumberOfComponents());
}
- // Constitutive relation.
- // read type;
- auto const constitutive_relation_config =
- //! \ogs_file_param{prj__processes__process__PHASE_FIELD__constitutive_relation}
- config.getConfigSubtree("constitutive_relation");
+ auto solid_constitutive_relations =
+ MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
+ parameters, config);
auto const phasefield_parameters_config =
//! \ogs_file_param{prj__processes__process__PHASE_FIELD__phasefield_parameters}
config.getConfigSubtree("phasefield_parameters");
- auto const type =
- //! \ogs_file_param{prj__processes__process__PHASE_FIELD__constitutive_relation__type}
- constitutive_relation_config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<MaterialLib::Solids::PhaseFieldExtension<DisplacementDim>>
- material = nullptr;
- if (type == "LinearElasticIsotropic")
- {
- const bool skip_type_checking = false;
- auto elastic_model =
- MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
- parameters, constitutive_relation_config, skip_type_checking);
- material = std::make_unique<
- MaterialLib::Solids::LinearElasticIsotropicPhaseField<
- DisplacementDim>>(
- std::move(elastic_model->getMaterialProperties()));
- }
- else
- {
- OGS_FATAL(
- "Cannot construct constitutive relation of given type \'%s\'.",
- type.c_str());
- }
-
// Residual stiffness
auto& residual_stiffness = findParameter<double>(
phasefield_parameters_config,
@@ -212,9 +187,11 @@ std::unique_ptr<Process> createPhaseFieldProcess(
((*crack_scheme == "propagating") || (*crack_scheme == "static")));
PhaseFieldProcessData<DisplacementDim> process_data{
- std::move(material), residual_stiffness, crack_resistance,
- crack_length_scale, kinetic_coefficient, solid_density,
- history_field, specific_body_force, propagating_crack,
+ materialIDs(mesh), std::move(solid_constitutive_relations),
+ residual_stiffness, crack_resistance,
+ crack_length_scale, kinetic_coefficient,
+ solid_density, history_field,
+ specific_body_force, propagating_crack,
crack_pressure};
SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index d8baae9231332189f29424d86ceb330afc52385d..32466594ff33e7be0aa24aeb0e7597186bfd45b5 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -12,7 +12,9 @@
#include <memory>
#include <vector>
-#include "MaterialLib/SolidModels/PhaseFieldExtension.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -68,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
@@ -126,6 +136,19 @@ public:
_ip_data.reserve(n_integration_points);
_secondary_data.N.resize(n_integration_points);
+ auto& solid_material =
+ MaterialLib::Solids::selectSolidConstitutiveRelation(
+ _process_data.solid_materials,
+ _process_data.material_ids,
+ e.getID());
+ if (!dynamic_cast<MaterialLib::Solids::LinearElasticIsotropic<
+ DisplacementDim> const*>(&solid_material))
+ {
+ OGS_FATAL(
+ "For phasefield process only linear elastic solid material "
+ "support is implemented.");
+ }
+
auto const shape_matrices =
initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, DisplacementDim>(
@@ -136,8 +159,7 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- // displacement (subscript u)
- _ip_data.emplace_back(*_process_data.material);
+ _ip_data.emplace_back(solid_material);
auto& ip_data = _ip_data[ip];
ip_data.integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
diff --git a/ProcessLib/PhaseField/PhaseFieldProcessData.h b/ProcessLib/PhaseField/PhaseFieldProcessData.h
index 75fdace07173ad99b6c67235b9c3a0bae855f470..069c0dc3254b850e690a3eebcc7e197e3bc4c74c 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcessData.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcessData.h
@@ -33,8 +33,11 @@ template <int DisplacementDim>
struct PhaseFieldProcessData
{
PhaseFieldProcessData(
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
- material_,
+ MeshLib::PropertyVector<int> const* const material_ids_,
+ std::map<int,
+ std::unique_ptr<
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+ solid_materials_,
Parameter<double> const& residual_stiffness_,
Parameter<double> const& crack_resistance_,
Parameter<double> const& crack_length_scale_,
@@ -42,8 +45,10 @@ struct PhaseFieldProcessData
Parameter<double> const& solid_density_,
Parameter<double>& history_field_,
Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_,
- bool propagating_crack_, bool crack_pressure_)
- : material{std::move(material_)},
+ bool const propagating_crack_,
+ bool const crack_pressure_)
+ : material_ids(material_ids_),
+ solid_materials{std::move(solid_materials_)},
residual_stiffness(residual_stiffness_),
crack_resistance(crack_resistance_),
crack_length_scale(crack_length_scale_),
@@ -56,21 +61,7 @@ struct PhaseFieldProcessData
{
}
- PhaseFieldProcessData(PhaseFieldProcessData&& other)
- : material{std::move(other.material)},
- residual_stiffness(other.residual_stiffness),
- crack_resistance(other.crack_resistance),
- crack_length_scale(other.crack_length_scale),
- kinetic_coefficient(other.kinetic_coefficient),
- solid_density(other.solid_density),
- history_field(other.history_field),
- specific_body_force(other.specific_body_force),
- dt(other.dt),
- t(other.t),
- propagating_crack(other.propagating_crack),
- crack_pressure(other.crack_pressure)
- {
- }
+ PhaseFieldProcessData(PhaseFieldProcessData&& other) = default;
//! Copies are forbidden.
PhaseFieldProcessData(PhaseFieldProcessData const&) = delete;
@@ -81,8 +72,11 @@ struct PhaseFieldProcessData
//! Assignments are not needed.
void operator=(PhaseFieldProcessData&&) = delete;
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
- material;
+ MeshLib::PropertyVector<int> const* const material_ids;
+
+ std::map<int, std::unique_ptr<
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+ solid_materials;
Parameter<double> const& residual_stiffness;
Parameter<double> const& crack_resistance;
Parameter<double> const& crack_length_scale;
diff --git a/ProcessLib/ThermoMechanicalPhaseField/CreateThermoMechanicalPhaseFieldProcess.cpp b/ProcessLib/ThermoMechanicalPhaseField/CreateThermoMechanicalPhaseFieldProcess.cpp
index be679ddd19e2982b0d8e3b3161d88f898ae37648..2deccc667b4c8c298298261b3168ba70814161d1 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/CreateThermoMechanicalPhaseFieldProcess.cpp
+++ b/ProcessLib/ThermoMechanicalPhaseField/CreateThermoMechanicalPhaseFieldProcess.cpp
@@ -9,9 +9,12 @@
#include "CreateThermoMechanicalPhaseFieldProcess.h"
-#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
-#include "MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h"
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
#include "ThermoMechanicalPhaseFieldProcess.h"
#include "ThermoMechanicalPhaseFieldProcessData.h"
@@ -88,8 +91,7 @@ std::unique_ptr<Process> createThermoMechanicalPhaseFieldProcess(
if (variable_ph->getNumberOfComponents() != 1)
{
OGS_FATAL(
- "Phase field process variable '%s' is not a scalar variable but "
- "has "
+ "Phasefield process variable '%s' is not a scalar variable but has "
"%d components.",
variable_ph->getName().c_str(),
variable_ph->getNumberOfComponents());
@@ -107,11 +109,9 @@ std::unique_ptr<Process> createThermoMechanicalPhaseFieldProcess(
variable_T->getNumberOfComponents());
}
- // Constitutive relation.
- // read type;
- auto const constitutive_relation_config =
- //! \ogs_file_param{prj__processes__process__THERMO_MECHANICAL_PHASE_FIELD__constitutive_relation}
- config.getConfigSubtree("constitutive_relation");
+ auto solid_constitutive_relations =
+ MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
+ parameters, config);
auto const phasefield_parameters_config =
//! \ogs_file_param{prj__processes__process__THERMO_MECHANICAL_PHASE_FIELD__phasefield_parameters}
@@ -121,30 +121,6 @@ std::unique_ptr<Process> createThermoMechanicalPhaseFieldProcess(
//! \ogs_file_param{prj__processes__process__THERMO_MECHANICAL_PHASE_FIELD__thermal_parameters}
config.getConfigSubtree("thermal_parameters");
- auto const type =
- //! \ogs_file_param{prj__processes__process__THERMO_MECHANICAL_PHASE_FIELD__constitutive_relation__type}
- constitutive_relation_config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<MaterialLib::Solids::PhaseFieldExtension<DisplacementDim>>
- material = nullptr;
- if (type == "LinearElasticIsotropic")
- {
- const bool skip_type_checking = false;
- auto elastic_model =
- MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
- parameters, constitutive_relation_config, skip_type_checking);
- material = std::make_unique<
- MaterialLib::Solids::LinearElasticIsotropicPhaseField<
- DisplacementDim>>(
- std::move(elastic_model->getMaterialProperties()));
- }
- else
- {
- OGS_FATAL(
- "Cannot construct constitutive relation of given type \'%s\'.",
- type.c_str());
- }
-
// Residual stiffness
auto& residual_stiffness = findParameter<double>(
phasefield_parameters_config,
@@ -234,7 +210,8 @@ std::unique_ptr<Process> createThermoMechanicalPhaseFieldProcess(
}
ThermoMechanicalPhaseFieldProcessData<DisplacementDim> process_data{
- std::move(material),
+ materialIDs(mesh),
+ std::move(solid_constitutive_relations),
residual_stiffness,
crack_resistance,
crack_length_scale,
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
index d7741e0e8abcde82671d2b55392df10fe9cf0ec7..774b7fc7fe79c9856c398de875050b98fcb981e1 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
@@ -12,7 +12,9 @@
#include <memory>
#include <vector>
-#include "MaterialLib/SolidModels/PhaseFieldExtension.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -70,7 +72,7 @@ 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 alpha,
@@ -79,13 +81,18 @@ struct IntegrationPointData final
{
eps_m.noalias() = eps - alpha * delta_T * Invariants::identity2;
- static_cast<
- MaterialLib::Solids::PhaseFieldExtension<DisplacementDim> const&>(
- solid_material)
- .calculateDegradedStress(
- t, x_position, eps_m, 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_m);
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
@@ -143,6 +150,20 @@ public:
_ip_data.reserve(n_integration_points);
_secondary_data.N.resize(n_integration_points);
+ auto& solid_material =
+ MaterialLib::Solids::selectSolidConstitutiveRelation(
+ _process_data.solid_materials,
+ _process_data.material_ids,
+ e.getID());
+ if (!dynamic_cast<MaterialLib::Solids::LinearElasticIsotropic<
+ DisplacementDim> const*>(&solid_material))
+ {
+ OGS_FATAL(
+ "For phasefield process only linear elastic solid material "
+ "support is implemented.");
+ }
+
+
auto const shape_matrices =
initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, DisplacementDim>(
@@ -153,8 +174,7 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- // displacement (subscript u)
- _ip_data.emplace_back(*_process_data.material);
+ _ip_data.emplace_back(solid_material);
auto& ip_data = _ip_data[ip];
ip_data.integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
@@ -336,8 +356,8 @@ private:
IntegrationMethod _integration_method;
MeshLib::Element const& _element;
- bool const _is_axially_symmetric;
SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+ bool const _is_axially_symmetric;
static const int temperature_index = 0;
static const int temperature_size = ShapeFunction::NPOINTS;
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcessData.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcessData.h
index cc667114bb69c59d0bc00f1f91f620077c23aedf..46dbe556f1ddfdee0d165e8fcc37d50f0650a72f 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcessData.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcessData.h
@@ -33,8 +33,11 @@ template <int DisplacementDim>
struct ThermoMechanicalPhaseFieldProcessData
{
ThermoMechanicalPhaseFieldProcessData(
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
- material_,
+ MeshLib::PropertyVector<int> const* const material_ids_,
+ std::map<int,
+ std::unique_ptr<
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+ solid_materials_,
Parameter<double> const& residual_stiffness_,
Parameter<double> const& crack_resistance_,
Parameter<double> const& crack_length_scale_,
@@ -46,7 +49,8 @@ struct ThermoMechanicalPhaseFieldProcessData
Parameter<double> const& residual_thermal_conductivity_,
double const reference_temperature_,
Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_)
- : material{std::move(material_)},
+ : material_ids(material_ids_),
+ solid_materials{std::move(solid_materials_)},
residual_stiffness(residual_stiffness_),
crack_resistance(crack_resistance_),
crack_length_scale(crack_length_scale_),
@@ -63,24 +67,7 @@ struct ThermoMechanicalPhaseFieldProcessData
}
ThermoMechanicalPhaseFieldProcessData(
- ThermoMechanicalPhaseFieldProcessData&& other)
- : material{std::move(other.material)},
- residual_stiffness(other.residual_stiffness),
- crack_resistance(other.crack_resistance),
- crack_length_scale(other.crack_length_scale),
- kinetic_coefficient(other.kinetic_coefficient),
- solid_density(other.solid_density),
- linear_thermal_expansion_coefficient(
- other.linear_thermal_expansion_coefficient),
- specific_heat_capacity(other.specific_heat_capacity),
- thermal_conductivity(other.thermal_conductivity),
- residual_thermal_conductivity(other.residual_thermal_conductivity),
- reference_temperature(other.reference_temperature),
- specific_body_force(other.specific_body_force),
- dt(other.dt),
- t(other.t)
- {
- }
+ ThermoMechanicalPhaseFieldProcessData&& other) = default;
//! Copies are forbidden.
ThermoMechanicalPhaseFieldProcessData(
@@ -92,8 +79,11 @@ struct ThermoMechanicalPhaseFieldProcessData
//! Assignments are not needed.
void operator=(ThermoMechanicalPhaseFieldProcessData&&) = delete;
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
- material;
+ MeshLib::PropertyVector<int> const* const material_ids;
+
+ std::map<int, std::unique_ptr<
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+ solid_materials;
Parameter<double> const& residual_stiffness;
Parameter<double> const& crack_resistance;
Parameter<double> const& crack_length_scale;