diff --git a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
index ce26e691ca9e976ffaa14987a6c142ed806310f4..a9d8f3fb07c2f8d81e21634c7eff87b410ffc0d4 100644
--- a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
@@ -17,6 +17,7 @@
#include "ParameterLib/Utils.h"
#include "PhaseFieldProcess.h"
#include "PhaseFieldProcessData.h"
+#include "ProcessLib/Common/HydroMechanics/CreateInitialStress.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h"
@@ -210,6 +211,10 @@ std::unique_ptr<Process> createPhaseFieldProcess(
phasefield_model_string.c_str());
}();
+ // Initial stress conditions
+ auto initial_stress = ProcessLib::createInitialStress<DisplacementDim>(
+ config, parameters, mesh);
+
auto const softening_curve = [&]
{
auto const softening_curve_string =
@@ -284,6 +289,7 @@ std::unique_ptr<Process> createPhaseFieldProcess(
crack_resistance,
crack_length_scale,
solid_density,
+ initial_stress,
specific_body_force,
pressurized_crack,
propagating_pressurized_crack,
diff --git a/ProcessLib/PhaseField/LocalAssemblerInterface.h b/ProcessLib/PhaseField/LocalAssemblerInterface.h
index 90d78743094e7ede4fcd842a1a6eb91d9480ba8f..bcd8cfe76baa2ad07750bfae2d01f86625f10b2d 100644
--- a/ProcessLib/PhaseField/LocalAssemblerInterface.h
+++ b/ProcessLib/PhaseField/LocalAssemblerInterface.h
@@ -23,6 +23,14 @@ struct PhaseFieldLocalAssemblerInterface
: public ProcessLib::LocalAssemblerInterface,
public NumLib::ExtrapolatableElement
{
+ virtual std::size_t setIPDataInitialConditions(
+ std::string_view const name, double const* values,
+ int const integration_order) = 0;
+
+ virtual std::vector<double> getSigma() const = 0;
+
+ virtual std::vector<double> getEpsilon() const = 0;
+
virtual std::vector<double> const& getIntPtSigma(
const double t,
std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 79e6b0a4aa499e58d1c0d9626d070bab3a90cfee..543149ea6b6b03391860b0a53c4edbbf7e0cd6f8 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -415,5 +415,73 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::computeEnergy(
surface_energy += element_surface_energy;
pressure_work += element_pressure_work;
}
+
+template <typename ShapeFunctionDisplacement, int DisplacementDim>
+std::size_t PhaseFieldLocalAssembler<
+ ShapeFunctionDisplacement,
+ DisplacementDim>::setIPDataInitialConditions(std::string_view const name,
+ double const* values,
+ int const integration_order)
+{
+ if (integration_order !=
+ static_cast<int>(_integration_method.getIntegrationOrder()))
+ {
+ OGS_FATAL(
+ "Setting integration point initial conditions; The integration "
+ "order of the local assembler for element {:d} is different from "
+ "the integration order in the initial condition.",
+ _element.getID());
+ }
+
+ if (name == "sigma")
+ {
+ if (_process_data.initial_stress.value != nullptr)
+ {
+ OGS_FATAL(
+ "Setting initial conditions for stress from integration "
+ "point data and from a parameter '{:s}' is not possible "
+ "simultaneously.",
+ _process_data.initial_stress.value->name);
+ }
+
+ return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+ values, _ip_data, &IpData::sigma);
+ }
+
+ return 0;
+}
+
+template <typename ShapeFunctionDisplacement, int DisplacementDim>
+std::vector<double> PhaseFieldLocalAssembler<ShapeFunctionDisplacement,
+ DisplacementDim>::getSigma() const
+{
+ return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
+ _ip_data, &IpData::sigma);
+}
+
+template <typename ShapeFunctionDisplacement, int DisplacementDim>
+std::vector<double> PhaseFieldLocalAssembler<
+ ShapeFunctionDisplacement, DisplacementDim>::getEpsilon() const
+{
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ std::vector<double> ip_epsilon_values;
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
+ ip_epsilon_values, n_integration_points, kelvin_vector_size);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ auto const& eps = _ip_data[ip].eps;
+ cache_mat.row(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
+ }
+
+ return ip_epsilon_values;
+}
+
} // namespace PhaseField
} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index 17bfe2b9f1e8764c671c5c3227814b97890a2185..c95506fe78a1c3bc361dca7579837b5e82b8a321 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -262,6 +262,12 @@ public:
}
}
+ /// Returns number of read integration points.
+ std::size_t setIPDataInitialConditions(
+ std::string_view const name,
+ double const* values,
+ int const integration_order) override;
+
void assemble(double const /*t*/, double const /*dt*/,
std::vector<double> const& /*local_x*/,
std::vector<double> const& /*local_x_prev*/,
@@ -287,7 +293,27 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- _ip_data[ip].pushBackState();
+ auto& ip_data = _ip_data[ip];
+
+ ParameterLib::SpatialPosition const x_position{
+ std::nullopt, _element.getID(), ip,
+ MathLib::Point3d(
+ NumLib::interpolateCoordinates<ShapeFunction,
+ ShapeMatricesType>(
+ _element, ip_data.N))};
+
+ /// Set initial stress from parameter.
+ if (_process_data.initial_stress.value)
+ {
+ ip_data.sigma =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector<
+ DisplacementDim>((*_process_data.initial_stress.value)(
+ std::numeric_limits<
+ double>::quiet_NaN() /* time independent */,
+ x_position));
+ }
+
+ ip_data.pushBackState();
}
}
@@ -327,6 +353,13 @@ public:
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
+ // TODO (naumov) This method is same as getIntPtSigma but for arguments and
+ // the ordering of the cache_mat.
+ // There should be only one.
+ std::vector<double> getSigma() const override;
+
+ std::vector<double> getEpsilon() const override;
+
private:
std::vector<double> const& getIntPtSigma(
const double /*t*/,
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.cpp b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
index 38233cf3db2ae52ae7e4165409d4562329d80eb2..3020a93570299bd68a0b1945c35cebdef71f1edd 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
@@ -17,6 +17,7 @@
#include "PhaseFieldFEM.h"
#include "ProcessLib/Process.h"
#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/SetIPDataInitialConditions.h"
namespace ProcessLib
{
@@ -47,6 +48,13 @@ PhaseFieldProcess<DisplacementDim>::PhaseFieldProcess(
_nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+
+ _integration_point_writer.emplace_back(
+ std::make_unique<MeshLib::IntegrationPointWriter>(
+ "sigma_ip",
+ static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
+ integration_order, _local_assemblers,
+ &LocalAssemblerInterface::getSigma));
}
template <int DisplacementDim>
@@ -159,6 +167,9 @@ void PhaseFieldProcess<DisplacementDim>::initializeConcreteProcess(
getExtrapolator(), _local_assemblers,
&LocalAssemblerInterface::getIntPtEpsilonTensile));
+ setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(),
+ _local_assemblers);
+
// Initialize local assemblers after all variables have been set.
GlobalExecutor::executeMemberOnDereferenced(
&LocalAssemblerInterface::initialize, _local_assemblers,
diff --git a/ProcessLib/PhaseField/PhaseFieldProcessData.h b/ProcessLib/PhaseField/PhaseFieldProcessData.h
index 65fb8ab0a580ad81c833a9220de1d55263bd8d79..21cc7e38a9c5dede4ac21e19b719c95657f5e23d 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcessData.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcessData.h
@@ -16,6 +16,7 @@
#include "MeshLib/PropertyVector.h"
#include "ParameterLib/Parameter.h"
+#include "ProcessLib/Common/HydroMechanics/InitialStress.h"
namespace MaterialLib
{
@@ -209,6 +210,10 @@ struct PhaseFieldProcessData
ParameterLib::Parameter<double> const& crack_resistance;
ParameterLib::Parameter<double> const& crack_length_scale;
ParameterLib::Parameter<double> const& solid_density;
+ /// Optional, initial stress field. A symmetric tensor, short vector
+ /// representation of length 4 or 6, ParameterLib::Parameter<double>.
+ InitialStress const initial_stress;
+
Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
bool pressurized_crack = false;
bool propagating_pressurized_crack = false;