diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 34f2d5a6588cee3a23762c7fa58100865f0ed970..d494bdec19bd28494d9f2c071b1cb6c4410c9763 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -42,6 +42,7 @@
#include "ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.h"
#include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h"
+#include "ProcessLib/PhaseField/CreatePhaseFieldProcess.h"
#include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/TES/CreateTESProcess.h"
@@ -399,6 +400,26 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
_process_variables, _parameters, integration_order,
process_config);
}
+ else if (type == "PHASE_FIELD")
+ {
+ switch (_mesh_vec[0]->getDimension())
+ {
+ case 2:
+ process =
+ ProcessLib::PhaseField::createPhaseFieldProcess<
+ 2>(*_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters,
+ integration_order, process_config);
+ break;
+ case 3:
+ process =
+ ProcessLib::PhaseField::createPhaseFieldProcess<
+ 3>(*_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters,
+ integration_order, process_config);
+ break;
+ }
+ }
else if (type == "SMALL_DEFORMATION")
{
switch (_mesh_vec[0]->getDimension())
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/c_PHASE_FIELD.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/c_PHASE_FIELD.md
new file mode 100644
index 0000000000000000000000000000000000000000..47f149e88d320fea48a74ea181e7c64b468290be
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/c_PHASE_FIELD.md
@@ -0,0 +1 @@
+Phase-field modelling for brittle fracture. It is implemented monolithically.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/i_constitutive_relation.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/i_constitutive_relation.md
new file mode 100644
index 0000000000000000000000000000000000000000..4a2457d3f4c6929e8670842215d3ab145424d306
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/i_constitutive_relation.md
@@ -0,0 +1 @@
+The constitutive relation for the brittle fracture.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/t_type.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/t_type.md
new file mode 100644
index 0000000000000000000000000000000000000000..c66714fc4af7a6a2e11c293f1bc98e75a4c7c09c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/constitutive_relation/t_type.md
@@ -0,0 +1,2 @@
+The type of constitutive relation available.
+1. Linear elastic
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/i_phasefield_parameters.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/i_phasefield_parameters.md
new file mode 100644
index 0000000000000000000000000000000000000000..fa401a8b3df36ff3689832b1cc7abf33053f5c53
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/i_phasefield_parameters.md
@@ -0,0 +1 @@
+A tag for phase-field parameters.
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_length_scale.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_length_scale.md
new file mode 100644
index 0000000000000000000000000000000000000000..2d624ad71c95f7e22d0fdda55333d60442394328
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_length_scale.md
@@ -0,0 +1 @@
+A length-scale parameter which controls the regularisation.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_resistance.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_resistance.md
new file mode 100644
index 0000000000000000000000000000000000000000..29068c95b59e62700e3d5634c3c05e80ecd9717a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_crack_resistance.md
@@ -0,0 +1 @@
+The critical Griffith-type fracture energy.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_history_field.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_history_field.md
new file mode 100644
index 0000000000000000000000000000000000000000..f90cb95d6a2d87cafe3c8fab042027fc716ad5eb
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_history_field.md
@@ -0,0 +1 @@
+A damage-driving history field which is associated with the maximum local tensile strain energy.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_kinetic_coefficient.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_kinetic_coefficient.md
new file mode 100644
index 0000000000000000000000000000000000000000..ad80c785d41859db77470c3594581fd9f8a370bf
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_kinetic_coefficient.md
@@ -0,0 +1 @@
+A kinetic coefficient which is defined to address rate-dependent cases.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_residual_stiffness.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_residual_stiffness.md
new file mode 100644
index 0000000000000000000000000000000000000000..33a233ac922b2cd95aeed1c9212997e2e63e6ebd
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/phasefield_parameters/t_residual_stiffness.md
@@ -0,0 +1 @@
+A residual stiffness which is used in numerical calculations.
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/i_process_variables.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/i_process_variables.md
new file mode 100644
index 0000000000000000000000000000000000000000..949045cf3b383e6fb962472086e6c261e08383df
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/i_process_variables.md
@@ -0,0 +1 @@
+The process variables for displacement and phasefield.
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_displacement.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_displacement.md
new file mode 100644
index 0000000000000000000000000000000000000000..b6e0f53b6b056cb721efd2ffa7ddfed60a5c5731
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_displacement.md
@@ -0,0 +1 @@
+Process variable name of displacement.
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_phasefield.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_phasefield.md
new file mode 100644
index 0000000000000000000000000000000000000000..cc5910921a5d5735b3fba88c729016ec67738aaa
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_phasefield.md
@@ -0,0 +1 @@
+Process variable name of phasefield.
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_solid_density.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_solid_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..54bbc62fe6cba8cce1195247a1c7339cde00233d
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_solid_density.md
@@ -0,0 +1 @@
+Reference solid density.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_specific_body_force.md
new file mode 100644
index 0000000000000000000000000000000000000000..dd77e37995b3f3628c82ccd1f99f6315167028f8
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_specific_body_force.md
@@ -0,0 +1 @@
+Specific body forces that are used to apply gravitational forces. A vector of displacement dimension's length.
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.h b/MaterialLib/SolidModels/LinearElasticIsotropic.h
index afa40e8a627d4172b060c6e5d61aa43ced82e16d..23cadbab63990e2c82b3f30ec229a87f1ce891f2 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.h
@@ -17,7 +17,7 @@ namespace MaterialLib
namespace Solids
{
template <int DisplacementDim>
-class LinearElasticIsotropic final : public MechanicsBase<DisplacementDim>
+class LinearElasticIsotropic : public MechanicsBase<DisplacementDim>
{
public:
/// Variables specific to the material model
@@ -47,6 +47,13 @@ public:
(2 * (1 + _poissons_ratio(t, x)[0]));
}
+ /// the bulk modulus.
+ double bulk_modulus(double const t, X const& x) const
+ {
+ return _youngs_modulus(t, x)[0] /
+ (3 * (1 - 2 * _poissons_ratio(t, x)[0]));
+ }
+
private:
P const& _youngs_modulus;
P const& _poissons_ratio;
@@ -101,7 +108,9 @@ public:
typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
material_state_variables) override;
-private:
+ MaterialProperties getMaterialProperties() {return _mp;}
+
+protected:
MaterialProperties _mp;
};
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e46fa45017a46097940740d8f3040700d73bb44a
--- /dev/null
+++ b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp
@@ -0,0 +1,21 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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
new file mode 100644
index 0000000000000000000000000000000000000000..1a5eeaf141d53b7c445321b2cfe2abbf6c314a76
--- /dev/null
+++ b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
@@ -0,0 +1,135 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "PhaseFieldExtension.h"
+#include "LinearElasticIsotropic.h"
+#include "KelvinVector.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+class LinearElasticIsotropicPhaseField final
+ : public LinearElasticIsotropic<DisplacementDim>,
+ public PhaseFieldExtension<DisplacementDim>
+{
+public:
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+ explicit LinearElasticIsotropicPhaseField(
+ typename LinearElasticIsotropic<
+ DisplacementDim>::MaterialProperties&& material_properties)
+ : LinearElasticIsotropic<DisplacementDim>(std::move(material_properties))
+ {
+ }
+
+ std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+ createMaterialStateVariables() 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) override
+ {
+ return LinearElasticIsotropic<DisplacementDim>::
+ integrateStress(t,
+ x,
+ dt,
+ eps_prev,
+ eps,
+ sigma_prev,
+ material_state_variables);
+ }
+
+ /** 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) const override
+ {
+ using Invariants =
+ MaterialLib::SolidModels::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();
+ }
+ 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);
+ }
+
+ sigma_real.noalias() = degradation * sigma_tensile + sigma_compressive;
+ return true;
+ }
+};
+
+extern template class LinearElasticIsotropicPhaseField<2>;
+extern template class LinearElasticIsotropicPhaseField<3>;
+
+} // namespace Solids
+} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/PhaseFieldExtension.h b/MaterialLib/SolidModels/PhaseFieldExtension.h
new file mode 100644
index 0000000000000000000000000000000000000000..77bb5a40027cce13779cbe2b1d985bf125916288
--- /dev/null
+++ b/MaterialLib/SolidModels/PhaseFieldExtension.h
@@ -0,0 +1,89 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+#include "MechanicsBase.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+
+template <int DisplacementDim>
+struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
+{
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix = ProcessLib::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) 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) 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);
+
+ sigma_tensile = sigma_tensile_;
+ sigma_compressive = sigma_compressive_;
+ C_tensile = C_tensile_;
+ C_compressive = C_compressive_;
+ sigma_real = sigma_real_;
+ return result;
+ }
+
+ virtual ~PhaseFieldExtension() = default;
+};
+
+} // namespace Solids
+} // namespace MaterialLib
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index e64cf8f67828ab3664c09f4704f28bfba4935a4d..f73380f5ba850aac7e4c504efe9028c539be0aae 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -17,6 +17,7 @@ APPEND_SOURCE_FILES(SOURCES LIE/SmallDeformation)
APPEND_SOURCE_FILES(SOURCES LIE/SmallDeformation/LocalAssembler)
APPEND_SOURCE_FILES(SOURCES LiquidFlow)
APPEND_SOURCE_FILES(SOURCES Parameter)
+APPEND_SOURCE_FILES(SOURCES PhaseField)
APPEND_SOURCE_FILES(SOURCES RichardsFlow)
APPEND_SOURCE_FILES(SOURCES SmallDeformation)
APPEND_SOURCE_FILES(SOURCES TES)
@@ -58,6 +59,7 @@ include(HydroMechanics/Tests.cmake)
include(LIE/HydroMechanics/Tests.cmake)
include(LIE/SmallDeformation/Tests.cmake)
include(LiquidFlow/Tests.cmake)
+include(PhaseField/Tests.cmake)
include(RichardsFlow/Tests.cmake)
include(SmallDeformation/Tests.cmake)
include(TES/Tests.cmake)
diff --git a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0dcadcb9de6e2bca3a7db6a1a405352e634432c5
--- /dev/null
+++ b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
@@ -0,0 +1,208 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "CreatePhaseFieldProcess.h"
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+
+#include "PhaseFieldProcess.h"
+#include "PhaseFieldProcessData.h"
+
+namespace ProcessLib
+{
+namespace PhaseField
+{
+template <int DisplacementDim>
+class PhaseFieldProcess;
+
+extern template class PhaseFieldProcess<2>;
+extern template class PhaseFieldProcess<3>;
+
+template <int DisplacementDim>
+std::unique_ptr<Process> createPhaseFieldProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{prj__processes__process__type}
+ config.checkConfigParameter("type", "PHASE_FIELD");
+ DBUG("Create PhaseFieldProcess.");
+
+ // Process variable.
+
+ //! \ogs_file_param{prj__processes__process__PHASE_FIELD__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+
+ auto process_variables = findProcessVariables(
+ variables, pv_config,
+ {//! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__process_variables__phasefield}
+ "phasefield",
+ //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__process_variables__displacement}
+ "displacement"});
+
+ DBUG("Associate displacement with process variable \'%s\'.",
+ process_variables[1].get().getName().c_str());
+
+ if (process_variables[1].get().getNumberOfComponents() != DisplacementDim)
+ {
+ OGS_FATAL(
+ "Number of components of the process variable '%s' is different "
+ "from the displacement dimension: got %d, expected %d",
+ process_variables[1].get().getName().c_str(),
+ process_variables[1].get().getNumberOfComponents(),
+ DisplacementDim);
+ }
+
+ DBUG("Associate phase field with process variable \'%s\'.",
+ process_variables[0].get().getName().c_str());
+ if (process_variables[0].get().getNumberOfComponents() != 1)
+ {
+ OGS_FATAL(
+ "Phase field process variable '%s' is not a scalar variable but has "
+ "%d components.",
+ process_variables[0].get().getName().c_str(),
+ process_variables[0].get().getNumberOfComponents(),
+ DisplacementDim);
+ }
+
+ // Constitutive relation.
+ // read type;
+ auto const constitutive_relation_config =
+ //! \ogs_file_param{prj__processes__process__PHASE_FIELD__constitutive_relation}
+ config.getConfigSubtree("constitutive_relation");
+
+ 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")
+ {
+ auto elastic_model = MaterialLib::Solids::createLinearElasticIsotropic<
+ DisplacementDim>(parameters, constitutive_relation_config);
+ 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,
+ //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__phasefield_parameters__residual_stiffness}
+ "residual_stiffness", parameters, 1);
+ DBUG("Use \'%s\' as residual stiffness.", residual_stiffness.name.c_str());
+
+ // Crack resistance
+ auto& crack_resistance = findParameter<double>(
+ phasefield_parameters_config,
+ //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__phasefield_parameters__crack_resistance}
+ "crack_resistance", parameters, 1);
+ DBUG("Use \'%s\' as crack resistance.", crack_resistance.name.c_str());
+
+ // Crack length scale
+ auto& crack_length_scale = findParameter<double>(
+ phasefield_parameters_config,
+ //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__phasefield_parameters__crack_length_scale}
+ "crack_length_scale", parameters, 1);
+ DBUG("Use \'%s\' as crack length scale.", crack_length_scale.name.c_str());
+
+ // Kinetic coefficient
+ auto& kinetic_coefficient = findParameter<double>(
+ phasefield_parameters_config,
+ //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__phasefield_parameters__kinetic_coefficient}
+ "kinetic_coefficient", parameters, 1);
+ DBUG("Use \'%s\' as kinetic coefficient.",
+ kinetic_coefficient.name.c_str());
+
+ // Solid density
+ auto& solid_density = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__solid_density}
+ "solid_density", parameters, 1);
+ DBUG("Use \'%s\' as solid density parameter.", solid_density.name.c_str());
+
+ // History field
+ auto& history_field = findParameter<double>(
+ phasefield_parameters_config,
+ //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__phasefield_parameters__history_field}
+ "history_field", parameters, 1);
+ DBUG("Use \'%s\' as history field.", history_field.name.c_str());
+
+ // Specific body force
+ Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
+ {
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__PHASE_FIELD__specific_body_force}
+ config.getConfigParameter<std::vector<double>>(
+ "specific_body_force");
+ if (specific_body_force.size() != DisplacementDim)
+ OGS_FATAL(
+ "The size of the specific body force vector does not match the "
+ "displacement dimension. Vector size is %d, displacement "
+ "dimension is %d",
+ specific_body_force.size(), DisplacementDim);
+
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+ }
+
+ PhaseFieldProcessData<DisplacementDim> process_data{
+ std::move(material), residual_stiffness, crack_resistance,
+ crack_length_scale, kinetic_coefficient, solid_density,
+ history_field, specific_body_force};
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"PhaseField_displacement"});
+
+ ProcessLib::parseSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+
+ return std::make_unique<PhaseFieldProcess<DisplacementDim>>(
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(process_data),
+ std::move(secondary_variables), std::move(named_function_caller));
+}
+
+template std::unique_ptr<Process> createPhaseFieldProcess<2>(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<Process> createPhaseFieldProcess<3>(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+
+} // namespace PhaseField
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/CreatePhaseFieldProcess.h b/ProcessLib/PhaseField/CreatePhaseFieldProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..770b8e7bd22f79f938d240c69151347b1d99fb93
--- /dev/null
+++ b/ProcessLib/PhaseField/CreatePhaseFieldProcess.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace PhaseField
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createPhaseFieldProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+
+} // namespace PhaseField
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
new file mode 100644
index 0000000000000000000000000000000000000000..f3be2299be43927ac8a26721c41bad5d6be40f80
--- /dev/null
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -0,0 +1,653 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 <memory>
+#include <vector>
+
+#include "MaterialLib/SolidModels/KelvinVector.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "PhaseFieldProcessData.h"
+
+namespace ProcessLib
+{
+namespace PhaseField
+{
+template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
+struct IntegrationPointData final
+{
+ explicit IntegrationPointData(
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+ : solid_material(solid_material),
+ material_state_variables(
+ solid_material.createMaterialStateVariables())
+ {
+ }
+
+ typename ShapeMatrixType::NodalRowVectorType N;
+ typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
+
+ typename BMatricesType::KelvinVectorType eps, eps_prev;
+
+ typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
+ sigma_real_prev, sigma_real;
+ double strain_energy_tensile;
+
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+ std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables> material_state_variables;
+
+ typename BMatricesType::KelvinMatrixType C_tensile, C_compressive;
+ double integration_weight;
+ double history_variable;
+ double history_variable_prev;
+
+ void pushBackState()
+ {
+ if (history_variable_prev <
+ history_variable)
+ {
+ history_variable_prev =
+ history_variable;
+ }
+ eps_prev = eps;
+ sigma_real_prev = sigma_real;
+ material_state_variables->pushBackState();
+ }
+
+ template <typename DisplacementVectorType>
+ void updateConstitutiveRelation(
+ double const t,
+ SpatialPosition const& x_position,
+ double const dt,
+ DisplacementVectorType const& /*u*/,
+ double const degradation)
+ {
+ static_cast<MaterialLib::Solids::PhaseFieldExtension<DisplacementDim>&>(
+ solid_material)
+ .calculateDegradedStress(t, x_position, eps, strain_energy_tensile,
+ sigma_tensile, sigma_compressive, C_tensile,
+ C_compressive, sigma_real, degradation);
+ }
+
+};
+
+/// Used by for extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+ std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
+};
+
+struct PhaseFieldLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+ virtual std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+class PhaseFieldLocalAssembler : public PhaseFieldLocalAssemblerInterface
+{
+public:
+ using ShapeMatricesType =
+ ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
+
+ // Types for displacement.
+ // (Higher order elements = ShapeFunction).
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+ using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>;
+
+ using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
+ using RhsVector = typename ShapeMatricesType::template VectorType<
+ ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim>;
+ using JacobianMatrix = typename ShapeMatricesType::template MatrixType<
+ ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim,
+ ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim>;
+
+ PhaseFieldLocalAssembler(PhaseFieldLocalAssembler const&) = delete;
+ PhaseFieldLocalAssembler(PhaseFieldLocalAssembler&&) = delete;
+
+ PhaseFieldLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ PhaseFieldProcessData<DisplacementDim>& process_data)
+ : _process_data(process_data),
+ _integration_method(integration_order),
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N.resize(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, DisplacementDim>(
+ e, is_axially_symmetric, _integration_method);
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ // displacement (subscript u)
+ _ip_data.emplace_back(*_process_data.material);
+ auto& ip_data = _ip_data[ip];
+ ip_data.integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
+
+ ip_data.eps.setZero(kelvin_vector_size);
+ ip_data.eps_prev.resize(kelvin_vector_size);
+ ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
+ ip_data.C_compressive.setZero(kelvin_vector_size,
+ kelvin_vector_size);
+ ip_data.sigma_tensile.setZero(kelvin_vector_size);
+ ip_data.sigma_compressive.setZero(kelvin_vector_size);
+ ip_data.history_variable =
+ process_data.history_field(0, x_position)[0];
+ ip_data.history_variable_prev =
+ process_data.history_field(0, x_position)[0];
+ ip_data.sigma_real.setZero(kelvin_vector_size);
+
+ ip_data.N = shape_matrices[ip].N;
+ ip_data.dNdx = shape_matrices[ip].dNdx;
+
+ _secondary_data.N[ip] = shape_matrices[ip].N;
+ }
+ }
+
+ void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& /*local_rhs_data*/) override
+ {
+ OGS_FATAL(
+ "PhaseFieldLocalAssembler: assembly without jacobian is not "
+ "implemented.");
+ }
+
+ void assembleWithJacobian(double const t,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ const double /*dxdot_dx*/, const double /*dx_dx*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_rhs_data,
+ std::vector<double>& local_Jac_data) override
+ {
+ auto const local_matrix_size = local_x.size();
+ assert(local_matrix_size == phasefield_size + displacement_size);
+
+ auto d = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ phasefield_size> const>(local_x.data() + phasefield_index,
+ phasefield_size);
+
+ auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ displacement_size> const>(local_x.data() + displacement_index,
+ displacement_size);
+
+ auto d_dot = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ phasefield_size> const>(local_xdot.data() + phasefield_index,
+ phasefield_size);
+
+ auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
+ local_Jac_data, local_matrix_size, local_matrix_size);
+
+ auto local_rhs = MathLib::createZeroedVector<RhsVector>(
+ local_rhs_data, local_matrix_size);
+
+ typename ShapeMatricesType::template MatrixType<displacement_size,
+ phasefield_size> Kud;
+ Kud.setZero(displacement_size, phasefield_size);
+
+ typename ShapeMatricesType::template MatrixType<phasefield_size,
+ displacement_size> Kdu;
+ Kdu.setZero(phasefield_size, displacement_size);
+
+ typename ShapeMatricesType::NodalMatrixType Kdd;
+ Kdd.setZero(phasefield_size, phasefield_size);
+
+ typename ShapeMatricesType::NodalMatrixType Ddd;
+ Ddd.setZero(phasefield_size, phasefield_size);
+
+ double const& dt = _process_data.dt;
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& dNdx = _ip_data[ip].dNdx;
+ auto const& N = _ip_data[ip].N;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ _element, N);
+ auto const& B = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
+ _is_axially_symmetric);
+
+ auto const& C_tensile = _ip_data[ip].C_tensile;
+ auto const& C_compressive = _ip_data[ip].C_compressive;
+
+ auto& eps = _ip_data[ip].eps;
+ auto const& strain_energy_tensile =
+ _ip_data[ip].strain_energy_tensile;
+ auto const& sigma_tensile = _ip_data[ip].sigma_tensile;
+
+ auto& history_variable = _ip_data[ip].history_variable;
+ auto& history_variable_prev = _ip_data[ip].history_variable_prev;
+
+ auto const& sigma_real = _ip_data[ip].sigma_real;
+
+ // Kdd_1 defines one term which both used in Kdd and local_rhs for
+ // phase field
+ double const gc = _process_data.crack_resistance(t, x_position)[0];
+ double const ls =
+ _process_data.crack_length_scale(t, x_position)[0];
+ typename ShapeMatricesType::NodalMatrixType const Kdd_1 =
+ dNdx.transpose() * 2 * gc * ls * dNdx;
+
+ //
+ // displacement equation, displacement part
+ //
+ double const k = _process_data.residual_stiffness(t, x_position)[0];
+ double const d_ip = N.dot(d);
+ double const degradation = d_ip * d_ip * (1 - k) + k;
+ eps.noalias() = B * u;
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ degradation);
+
+ local_Jac.template block<displacement_size, displacement_size>(
+ displacement_index, displacement_index)
+ .noalias() += B.transpose() *
+ (degradation * C_tensile + C_compressive) * B * w;
+
+ typename ShapeMatricesType::template MatrixType<DisplacementDim,
+ displacement_size>
+ N_u = ShapeMatricesType::template MatrixType<
+ DisplacementDim,
+ displacement_size>::Zero(DisplacementDim,
+ displacement_size);
+
+ for (int i = 0; i < DisplacementDim; ++i)
+ N_u.template block<1, displacement_size / DisplacementDim>(
+ i, i * displacement_size / DisplacementDim)
+ .noalias() = N;
+
+ auto const rho_sr = _process_data.solid_density(t, x_position)[0];
+ auto const& b = _process_data.specific_body_force;
+ local_rhs.template block<displacement_size, 1>(displacement_index, 0)
+ .noalias() -=
+ (B.transpose() * sigma_real - N_u.transpose() * rho_sr * b) * w;
+
+ //
+ // displacement equation, phasefield part
+ //
+ Kud.noalias() += B.transpose() * 2 * d_ip * sigma_tensile * N * w;
+
+ double const d_dot_ip = N.dot(d_dot);
+
+ if (history_variable_prev < strain_energy_tensile)
+ {
+ history_variable = strain_energy_tensile;
+ Kdu.noalias() = Kud.transpose();
+ }
+ else
+ {
+ history_variable = history_variable_prev;
+ }
+
+ //
+ // phasefield equation, phasefield part.
+ //
+ Kdd.noalias() += (Kdd_1 + N.transpose() * 2 * history_variable * N +
+ N.transpose() * 0.5 * gc / ls * N) *
+ w;
+ double const M =
+ _process_data.kinetic_coefficient(t, x_position)[0];
+ local_rhs.template segment<phasefield_size>(phasefield_index)
+ .noalias() -= (N.transpose() * d_dot_ip / M + Kdd_1 * d +
+ N.transpose() * d_ip * 2 * history_variable -
+ N.transpose() * 0.5 * gc / ls * (1 - d_ip)) *
+ w;
+
+ Ddd.noalias() += N.transpose() / M * N * w;
+ }
+ // displacement equation, phasefield part
+ local_Jac.template block<displacement_size, phasefield_size>(
+ displacement_index, phasefield_index)
+ .noalias() += Kud;
+
+ // phasefield equation, phasefield part.
+ local_Jac.template block<phasefield_size, phasefield_size>(
+ phasefield_index, phasefield_index)
+ .noalias() += Kdd + Ddd / dt;
+
+ // phasefield equation, displacement part.
+ local_Jac.template block<phasefield_size, displacement_size>(
+ phasefield_index, displacement_index)
+ .noalias() += Kdu;
+
+ }
+
+ void preTimestepConcrete(std::vector<double> const& /*local_x*/,
+ double const /*t*/,
+ double const /*delta_t*/) override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data[ip].pushBackState();
+ }
+ }
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _secondary_data.N[integration_point];
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+ std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 0);
+ }
+
+ std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 1);
+ }
+
+ std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 2);
+ }
+
+ std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 3);
+ }
+
+ std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtSigma(cache, 4);
+ }
+
+ std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtSigma(cache, 5);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 0);
+ }
+
+ std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 1);
+ }
+
+ std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 2);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 3);
+ }
+
+ std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtEpsilon(cache, 4);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtEpsilon(cache, 5);
+ }
+
+private:
+ std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
+ std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ if (component < 3) // xx, yy, zz components
+ cache.push_back(ip_data.sigma_real[component]);
+ else // mixed xy, yz, xz components
+ cache.push_back(ip_data.sigma_real[component] / std::sqrt(2));
+ }
+
+ return cache;
+ }
+
+ std::vector<double> const& getIntPtEpsilon(
+ std::vector<double>& cache, std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ if (component < 3) // xx, yy, zz components
+ cache.push_back(ip_data.eps[component]);
+ else // mixed xy, yz, xz components
+ cache.push_back(ip_data.eps[component] / std::sqrt(2));
+ }
+
+ return cache;
+ }
+
+ PhaseFieldProcessData<DisplacementDim>& _process_data;
+
+ std::vector<
+ IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>,
+ Eigen::aligned_allocator<IntegrationPointData<
+ BMatricesType, ShapeMatricesType, DisplacementDim>>>
+ _ip_data;
+
+ IntegrationMethod _integration_method;
+ MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
+ SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+
+ static const int phasefield_index = 0;
+ static const int phasefield_size = ShapeFunction::NPOINTS;
+ static const int displacement_index = ShapeFunction::NPOINTS;
+ static const int displacement_size =
+ ShapeFunction::NPOINTS * DisplacementDim;
+ static const int kelvin_vector_size =
+ KelvinVectorDimensions<DisplacementDim>::value;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim, int DisplacementDim>
+class LocalAssemblerData final
+ : public PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>
+{
+public:
+ LocalAssemblerData(LocalAssemblerData const&) = delete;
+ LocalAssemblerData(LocalAssemblerData&&) = delete;
+
+ LocalAssemblerData(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ PhaseFieldProcessData<DisplacementDim>& process_data)
+ : PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>(
+ e, local_matrix_size, is_axially_symmetric, integration_order,
+ process_data)
+ {
+ }
+};
+
+} // namespace PhaseField
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.cpp b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e3243e76ceb7dba633fd6c560456ef52594e8ee
--- /dev/null
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
@@ -0,0 +1,20 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "PhaseFieldProcess.h"
+
+namespace ProcessLib
+{
+namespace PhaseField
+{
+template class PhaseFieldProcess<2>;
+template class PhaseFieldProcess<3>;
+
+} // namespace PhaseField
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.h b/ProcessLib/PhaseField/PhaseFieldProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..df885f24062e0c3bd56f18358dfe47a599ebf675
--- /dev/null
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.h
@@ -0,0 +1,244 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 <cassert>
+
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Process.h"
+#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
+
+#include "PhaseFieldFEM.h"
+#include "PhaseFieldProcessData.h"
+
+namespace ProcessLib
+{
+namespace PhaseField
+{
+
+/**
+ * \brief A class to simulate mechanical fracturing process
+ * using phase-field approach in solids described by
+ *
+ * \f[
+ * \mathrm{div} \left[ \left(d^2 + k \right) \boldsymbol{\sigma}_0^+
+ * + \boldsymbol{\sigma}_0^- \right] +
+ * \varrho \boldsymbol{b} = \boldsymbol{0}
+ * \f]
+ * \f[
+ * 2d H(\boldsymbol{\epsilon}_\mathrm{el})
+ * - \frac{1 - d}{2 \varepsilon} g_\mathrm{c}
+ * - 2 \varepsilon g_\mathrm{c} \mathrm{div}(\mathrm{grad} d) = 0
+ * \f]
+ * where
+ * \f{eqnarray*}{
+ * &d:& \mbox{order parameter,}\\
+ * &\varrho:& \mbox{density,}\\
+ * &H:& \mbox{history field,}\\
+ * &g_\mathrm{c}:& \mbox{fracture energy,}\\
+ * &\varepsilon:& \mbox{length scale}\\
+ * \f}
+ *
+ * Detailed model description can refer
+ * <a href="Miao_Biot2017.pdff" target="_blank"><b>Phase field method</b></a>
+ */
+template <int DisplacementDim>
+class PhaseFieldProcess final : public Process
+{
+ using Base = Process;
+
+public:
+ PhaseFieldProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+ jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::reference_wrapper<ProcessVariable>>&&
+ process_variables,
+ PhaseFieldProcessData<DisplacementDim>&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller)
+ : Process(mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables),
+ std::move(named_function_caller)),
+ _process_data(std::move(process_data))
+ {
+ }
+
+ //! \name ODESystem interface
+ //! @{
+
+ bool isLinear() const override { return false; }
+ //! @}
+
+private:
+ using LocalAssemblerInterface = PhaseFieldLocalAssemblerInterface;
+
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override
+ {
+ ProcessLib::SmallDeformation::createLocalAssemblers<
+ DisplacementDim, PhaseFieldLocalAssembler>(
+ mesh.getElements(), dof_table, _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
+
+ // TODO move the two data members somewhere else.
+ // for extrapolation of secondary variables
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
+ all_mesh_subsets_single_component.emplace_back(
+ _mesh_subset_all_nodes.get());
+ _local_to_global_index_map_single_component.reset(
+ new NumLib::LocalToGlobalIndexMap(
+ std::move(all_mesh_subsets_single_component),
+ // by location order is needed for output
+ NumLib::ComponentOrder::BY_LOCATION));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xx",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtSigmaXX));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_yy",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtSigmaYY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_zz",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtSigmaZZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xy",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtSigmaXY));
+
+ if (DisplacementDim == 3)
+ {
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xz",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtSigmaXZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_yz",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtSigmaYZ));
+ }
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xx",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtEpsilonXX));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_yy",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtEpsilonYY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_zz",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtEpsilonZZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xy",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::getIntPtEpsilonXY));
+ if (DisplacementDim == 3)
+ {
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_yz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::
+ getIntPtEpsilonYZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &PhaseFieldLocalAssemblerInterface::
+ getIntPtEpsilonXZ));
+ }
+ }
+
+ void assembleConcreteProcess(
+ const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, StaggeredCouplingTerm const& coupling_term) override
+ {
+ DBUG("Assemble PhaseFieldProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble,
+ _local_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+ coupling_term);
+ }
+
+ void assembleWithJacobianConcreteProcess(
+ const double t, GlobalVector const& x, GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+ StaggeredCouplingTerm const& coupling_term) override
+ {
+ // DBUG("AssembleJacobian PhaseFieldProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, *_local_to_global_index_map, t, x, xdot,
+ dxdot_dx, dx_dx, M, K, b, Jac, coupling_term);
+ }
+
+ void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+ double const dt) override
+ {
+ DBUG("PreTimestep PhaseFieldProcess.");
+
+ _process_data.dt = dt;
+ _process_data.t = t;
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &PhaseFieldLocalAssemblerInterface::preTimestep, _local_assemblers,
+ *_local_to_global_index_map, x, t, dt);
+ }
+
+ void postTimestepConcreteProcess(GlobalVector const& x) override
+ {
+ DBUG("PostTimestep PhaseFieldProcess.");
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &PhaseFieldLocalAssemblerInterface::postTimestep, _local_assemblers,
+ *_local_to_global_index_map, x);
+ }
+
+private:
+ PhaseFieldProcessData<DisplacementDim> _process_data;
+
+ std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+ _local_to_global_index_map_single_component;
+};
+
+} // namespace PhaseField
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldProcessData.h b/ProcessLib/PhaseField/PhaseFieldProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..bea1a24b92ecc4cb4677f482043bf87b41a251fc
--- /dev/null
+++ b/ProcessLib/PhaseField/PhaseFieldProcessData.h
@@ -0,0 +1,82 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+namespace PhaseField
+{
+template <int DisplacementDim>
+struct PhaseFieldProcessData
+{
+ PhaseFieldProcessData(
+ std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
+ material_,
+ Parameter<double> const& residual_stiffness_,
+ Parameter<double> const& crack_resistance_,
+ Parameter<double> const& crack_length_scale_,
+ Parameter<double> const& kinetic_coefficient_,
+ Parameter<double> const& solid_density_,
+ Parameter<double>& history_field_,
+ Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_)
+ : material{std::move(material_)},
+ residual_stiffness(residual_stiffness_),
+ crack_resistance(crack_resistance_),
+ crack_length_scale(crack_length_scale_),
+ kinetic_coefficient(kinetic_coefficient_),
+ solid_density(solid_density_),
+ history_field(history_field_),
+ specific_body_force(specific_body_force_)
+ {
+ }
+
+ 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)
+ {
+ }
+
+ //! Copies are forbidden.
+ PhaseFieldProcessData(PhaseFieldProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(PhaseFieldProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(PhaseFieldProcessData&&) = delete;
+
+ std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+ material;
+ Parameter<double> const& residual_stiffness;
+ Parameter<double> const& crack_resistance;
+ Parameter<double> const& crack_length_scale;
+ Parameter<double> const& kinetic_coefficient;
+ Parameter<double> const& solid_density;
+ Parameter<double>& history_field;
+ Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force;
+ double dt;
+ double t;
+};
+
+} // namespace PhaseField
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/Tests.cmake b/ProcessLib/PhaseField/Tests.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..8a267e8658df4503c2a6029e07c0438d47635fc5
--- /dev/null
+++ b/ProcessLib/PhaseField/Tests.cmake
@@ -0,0 +1,13 @@
+AddTest(
+ NAME PhaseField_3D_Unconfined_Compression
+ PATH PhaseField
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e0.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ ABSTOL 1e-6 RELTOL 1e-6
+ DIFF_DATA
+ expected_cube_1e0_pcs_0_ts_10000_t_1.000000.vtu cube_1e0_pcs_0_ts_10000_t_1.000000.vtu displacement displacement
+ expected_cube_1e0_pcs_0_ts_10000_t_1.000000.vtu cube_1e0_pcs_0_ts_10000_t_1.000000.vtu phasefield phasefield
+ )
diff --git a/web/content/docs/benchmarks/phase-field/Miao_Biot2017.pdf b/web/content/docs/benchmarks/phase-field/Miao_Biot2017.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..053152c7ecbae216aa541594d7b28e48e73454cf
Binary files /dev/null and b/web/content/docs/benchmarks/phase-field/Miao_Biot2017.pdf differ
diff --git a/web/content/docs/benchmarks/phase-field/beam.png b/web/content/docs/benchmarks/phase-field/beam.png
new file mode 100644
index 0000000000000000000000000000000000000000..1a9bf92fef409308ac9bb976eaeff5fb4e268019
Binary files /dev/null and b/web/content/docs/benchmarks/phase-field/beam.png differ
diff --git a/web/content/docs/benchmarks/phase-field/beam_d.png b/web/content/docs/benchmarks/phase-field/beam_d.png
new file mode 100644
index 0000000000000000000000000000000000000000..8a06309b0455bb898b5558b809489b4c7c5b261b
Binary files /dev/null and b/web/content/docs/benchmarks/phase-field/beam_d.png differ
diff --git a/web/content/docs/benchmarks/phase-field/beam_u.png b/web/content/docs/benchmarks/phase-field/beam_u.png
new file mode 100644
index 0000000000000000000000000000000000000000..1fcd4d021010aed55c5c225b8d8aa3c837c0294b
Binary files /dev/null and b/web/content/docs/benchmarks/phase-field/beam_u.png differ
diff --git a/web/content/docs/benchmarks/phase-field/error_u.png b/web/content/docs/benchmarks/phase-field/error_u.png
new file mode 100644
index 0000000000000000000000000000000000000000..edaa0b050dbb798caabb64f1bbb35be39aeaa476
Binary files /dev/null and b/web/content/docs/benchmarks/phase-field/error_u.png differ
diff --git a/web/content/docs/benchmarks/phase-field/phasefield.md b/web/content/docs/benchmarks/phase-field/phasefield.md
new file mode 100644
index 0000000000000000000000000000000000000000..af1d34b6b75c70f58b3b3f1e61d820acd6c4adc9
--- /dev/null
+++ b/web/content/docs/benchmarks/phase-field/phasefield.md
@@ -0,0 +1,50 @@
++++
+project = "https://github.com/ufz/ogs-data/blob/master/PhaseField/beam3d.prj"
+author = "Xing-Yuan Miao"
+date = "2017-05-19T09:10:33+01:00"
+title = "Crack beam under tension"
+weight = 158
+
+[menu]
+ [menu.benchmarks]
+ parent = "phase-field"
+
++++
+
+{{< project-link >}}
+
+## Problem description
+
+We solve a homogeneous beam model under a given displacement loading. The length of the beam is 2\,mm. Detailed model description can refer [this PDF](../Miao_Biot2017.pdf).
+## Results and evaluation
+
+Results show crack phase-field and displacement field distributions through the length of the beam:
+
+{{< img src="../beam.png" >}}
+{{< img src="../beam_d.png" >}}
+{{< img src="../beam_u.png" >}}
+
+For highlighting the deviation between the analytical and numerical solution, we provide the absolute error of the analytical solution and numerical simulation as follows:
+{{< img src="../error_u.png" >}}
+
+The analytical solution is:
+$$
+\begin{equation}
+d (x) = 1 - {\mathrm{e}}^{\frac{- |x|}{2 \varepsilon}}
+\end{equation}
+$$
+$$
+\begin{equation}
+u (x) = \dfrac{\sigma}{E} \int_0^x \dfrac{1}{d (x)^2 + k} \mathrm{d}x
+\end{equation}
+$$
+with
+$$
+\begin{equation}
+\sigma = \dfrac{E u_0}{I (\varepsilon, k)}
+\end{equation}
+$$
+\begin{equation}
+I (\varepsilon, k) = \int_0^1 \dfrac{1}{d (x)^2 + k} \mathrm{d}x
+\end{equation}
+$$
\ No newline at end of file