diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/PhaseFieldIrreversibleDamageOracleBoundaryCondition/c_PhaseFieldIrreversibleDamageOracleBoundaryCondition.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/PhaseFieldIrreversibleDamageOracleBoundaryCondition/c_PhaseFieldIrreversibleDamageOracleBoundaryCondition.md
new file mode 100644
index 0000000000000000000000000000000000000000..24deec36a1cd164d6ca9d648e1dabb1308a20ce8
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/PhaseFieldIrreversibleDamageOracleBoundaryCondition/c_PhaseFieldIrreversibleDamageOracleBoundaryCondition.md
@@ -0,0 +1 @@
+Phase field irreversibility is set. When its value becomes smaller than 0.05, it will be set to 0 in dirichlet boundary condition.
diff --git a/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_hydro_crack_scheme.md b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_hydro_crack_scheme.md
new file mode 100644
index 0000000000000000000000000000000000000000..e938d33adaba54df328d623ca8c11f732bf1902b
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_hydro_crack_scheme.md
@@ -0,0 +1 @@
+In case of static, pressure is set to 1.0. For propagating crack, evolution of pressure is found through enforcing the crack volume balance.
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
index 7bce42d01c2ff387bba59e81a4b69526b6baa00a..7314b1130fd3755cf758b056817d43e975eaa148 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
@@ -94,7 +94,8 @@ public:
KelvinMatrix& C_tensile,
KelvinMatrix& C_compressive,
KelvinVector& sigma_real,
- double const degradation) const override
+ double const degradation,
+ double& elastic_energy) const override
{
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
// calculation of deviatoric parts
@@ -120,6 +121,7 @@ public:
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
{
@@ -129,6 +131,8 @@ public:
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;
diff --git a/MaterialLib/SolidModels/PhaseFieldExtension.h b/MaterialLib/SolidModels/PhaseFieldExtension.h
index 4ef4dbfea6272bfaecdd6b7e6a000adfec33754c..425da16cf2b6ad27b77d878dd100da01705007d6 100644
--- a/MaterialLib/SolidModels/PhaseFieldExtension.h
+++ b/MaterialLib/SolidModels/PhaseFieldExtension.h
@@ -32,7 +32,8 @@ struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
KelvinMatrix& C_tensile,
KelvinMatrix& C_compressive,
KelvinVector& sigma_real,
- double const degradation) const = 0;
+ double const degradation,
+ double& elastic_energy) const = 0;
/// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
/// constitutive relation compute function.
@@ -48,7 +49,8 @@ struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
C_compressive,
Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_real,
- double const degradation) const
+ double const degradation,
+ double& elastic_energy) const
{
// TODO Avoid copies of data:
// Using MatrixBase<Derived> not possible because template functions
@@ -71,7 +73,8 @@ struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
C_tensile_,
C_compressive_,
sigma_real_,
- degradation);
+ degradation,
+ elastic_energy);
sigma_tensile = sigma_tensile_;
sigma_compressive = sigma_compressive_;
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 64a633bfc0d69d172f59c91e8efe11062b80ad0f..ead1f9bbc2bf9820a653d60844c20dd91bc45737 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -18,6 +18,7 @@
#include "NonuniformDirichletBoundaryCondition.h"
#include "NonuniformNeumannBoundaryCondition.h"
#include "NormalTractionBoundaryCondition.h"
+#include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
namespace ProcessLib
@@ -62,7 +63,6 @@ BoundaryConditionBuilder::createBoundaryCondition(
config, dof_table, mesh, variable_id, integration_order,
shapefunction_order);
}
-
//
// Special boundary conditions
//
@@ -72,6 +72,12 @@ BoundaryConditionBuilder::createBoundaryCondition(
config, dof_table, mesh, variable_id, integration_order,
shapefunction_order, parameters);
}
+ if (type == "PhaseFieldIrreversibleDamageOracleBoundaryCondition")
+ {
+ return createPhaseFieldIrreversibleDamageOracleBoundaryCondition(
+ config, dof_table, mesh, variable_id, integration_order,
+ shapefunction_order, parameters);
+ }
OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str());
}
@@ -230,6 +236,21 @@ BoundaryConditionBuilder::createNormalTractionBoundaryCondition(
parameters);
}
+std::unique_ptr<BoundaryCondition> BoundaryConditionBuilder::
+ createPhaseFieldIrreversibleDamageOracleBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table,
+ const MeshLib::Mesh& mesh, const int variable_id,
+ const unsigned /*integration_order*/,
+ const unsigned /*shapefunction_order*/,
+ const std::vector<
+ std::unique_ptr<ProcessLib::ParameterBase>>& /*parameters*/)
+{
+ return ProcessLib::
+ createPhaseFieldIrreversibleDamageOracleBoundaryCondition(
+ config.config, dof_table, mesh, variable_id, *config.component_id);
+}
+
std::vector<MeshLib::Element*> BoundaryConditionBuilder::getClonedElements(
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
GeoLib::GeoObject const& geometry)
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index 6f02c9ae29e159f00418d4e53c72a11febdcae89..a654c57659168f5c3447bae767c4e3c55e3f0cda 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -62,6 +62,11 @@ public:
// Therefore there is nothing to do here.
}
+ virtual void preTimestep(const double /*t*/, GlobalVector const& /*x*/)
+ {
+ // A hook added for solution dependent dirichlet
+ }
+
virtual ~BoundaryCondition() = default;
};
@@ -127,6 +132,16 @@ protected:
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters);
+ virtual std::unique_ptr<BoundaryCondition>
+ createPhaseFieldIrreversibleDamageOracleBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table,
+ const MeshLib::Mesh& mesh, const int variable_id,
+ const unsigned /*integration_order*/,
+ const unsigned /*shapefunction_order*/,
+ const std::vector<
+ std::unique_ptr<ProcessLib::ParameterBase>>& /*parameters*/);
+
static std::vector<MeshLib::Element*> getClonedElements(
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
GeoLib::GeoObject const& geometry);
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
index 42aec640e3fdd9072e9d1216b5a16705586f2d1c..679f76579600208e9c06ea9d56425fb3bfa7ce2d 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
@@ -44,6 +44,14 @@ public:
NumLib::LocalToGlobalIndexMap const& dof_table,
unsigned const integration_order);
+ void preTimestep(const double t, GlobalVector const& x)
+ {
+ for (auto const& bc_ptr : _boundary_conditions)
+ {
+ bc_ptr->preTimestep(t, x);
+ }
+ }
+
private:
mutable std::vector<NumLib::IndexValueVector<GlobalIndexType>> _dirichlet_bcs;
std::vector<std::unique_ptr<BoundaryCondition>> _boundary_conditions;
diff --git a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bc72d5b47334e51860f4ee502c9b6b6797c7312e
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
@@ -0,0 +1,85 @@
+/**
+ * \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 "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
+
+#include <algorithm>
+#include <logog/include/logog.hpp>
+#include <vector>
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+void PhaseFieldIrreversibleDamageOracleBoundaryCondition::getEssentialBCValues(
+ const double /*t*/,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+{
+ SpatialPosition pos;
+
+ bc_values.ids.clear();
+ bc_values.values.clear();
+
+ // convert mesh node ids to global index for the given component
+ bc_values.ids.reserve(bc_values.ids.size() + _bc_values.ids.size());
+ bc_values.values.reserve(bc_values.values.size() +
+ _bc_values.values.size());
+
+ std::copy(_bc_values.ids.begin(), _bc_values.ids.end(),
+ std::back_inserter(bc_values.ids));
+ std::copy(_bc_values.values.begin(), _bc_values.values.end(),
+ std::back_inserter(bc_values.values));
+}
+
+// update new values and corresponding indices.
+void PhaseFieldIrreversibleDamageOracleBoundaryCondition::preTimestep(
+ const double /*t*/, const GlobalVector& x)
+{
+ // phase-field variable is considered irreversible if it loses more than 95%
+ // of the stiffness, which is a widely used threshold.
+ double irreversibleDamage = 0.05;
+
+ _bc_values.ids.clear();
+ _bc_values.values.clear();
+
+ auto const mesh_id = _mesh.getID();
+ auto const& nodes = _mesh.getNodes();
+ for (auto const* n : nodes)
+ {
+ std::size_t node_id = n->getID();
+ MeshLib::Location l(mesh_id, MeshLib::MeshItemType::Node, node_id);
+ const auto g_idx =
+ _dof_table.getGlobalIndex(l, _variable_id, _component_id);
+
+ if (x[node_id] <= irreversibleDamage)
+ {
+ _bc_values.ids.emplace_back(g_idx);
+ _bc_values.values.emplace_back(0.0);
+ }
+ }
+}
+
+std::unique_ptr<PhaseFieldIrreversibleDamageOracleBoundaryCondition>
+createPhaseFieldIrreversibleDamageOracleBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ NumLib::LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
+ int const variable_id, int const component_id)
+{
+ DBUG(
+ "Constructing PhaseFieldIrreversibleDamageOracleBoundaryCondition from "
+ "config.");
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+ config.checkConfigParameter(
+ "type", "PhaseFieldIrreversibleDamageOracleBoundaryCondition");
+
+ return std::make_unique<
+ PhaseFieldIrreversibleDamageOracleBoundaryCondition>(
+ dof_table, mesh, variable_id, component_id);
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..9afa9f2426c3ceeb3016736c3ea7f5538cf87788
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
@@ -0,0 +1,66 @@
+/**
+ * \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 "BoundaryCondition.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/IndexValueVector.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace ProcessLib
+{
+class PhaseFieldIrreversibleDamageOracleBoundaryCondition final
+ : public BoundaryCondition
+{
+public:
+ PhaseFieldIrreversibleDamageOracleBoundaryCondition(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh, int const variable_id,
+ int const component_id)
+ : _dof_table(dof_table),
+ _mesh(mesh),
+ _variable_id(variable_id),
+ _component_id(component_id)
+ {
+ if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
+ component_id >=
+ dof_table.getNumberOfVariableComponents(variable_id))
+ {
+ OGS_FATAL(
+ "Variable id or component id too high. Actual values: (%d, "
+ "%d), "
+ "maximum values: (%d, %d).",
+ variable_id, component_id, dof_table.getNumberOfVariables(),
+ dof_table.getNumberOfVariableComponents(variable_id));
+ }
+ }
+
+ void getEssentialBCValues(
+ const double t,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
+
+ void preTimestep(const double t, const GlobalVector& x) override;
+
+private:
+ NumLib::LocalToGlobalIndexMap const& _dof_table;
+ MeshLib::Mesh const& _mesh;
+ int const _variable_id;
+ int const _component_id;
+
+ NumLib::IndexValueVector<GlobalIndexType> _bc_values;
+};
+
+std::unique_ptr<PhaseFieldIrreversibleDamageOracleBoundaryCondition>
+createPhaseFieldIrreversibleDamageOracleBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ NumLib::LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
+ int const variable_id, int const component_id);
+
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
index 8116905574936e1471daebded5a895bc4545687b..674bb831fa55c0baada25f33375e133cdd087fbc 100644
--- a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
@@ -22,7 +22,6 @@ namespace ProcessLib
{
namespace PhaseField
{
-
template <int DisplacementDim>
std::unique_ptr<Process> createPhaseFieldProcess(
MeshLib::Mesh& mesh,
@@ -56,9 +55,9 @@ std::unique_ptr<Process> createPhaseFieldProcess(
auto per_process_variables = findProcessVariables(
variables, pv_config,
{//! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__process_variables__phasefield}
- "phasefield",
+ "phasefield",
//! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__process_variables__displacement}
- "displacement"});
+ "displacement"});
variable_ph = &per_process_variables[0].get();
variable_u = &per_process_variables[1].get();
process_variables.push_back(std::move(per_process_variables));
@@ -66,14 +65,14 @@ std::unique_ptr<Process> createPhaseFieldProcess(
else // staggered scheme.
{
using namespace std::string_literals;
- for (auto const& variable_name : {"phasefield"s, "displacement"s})
+ for (auto const& variable_name : {"displacement"s, "phasefield"s})
{
auto per_process_variables =
findProcessVariables(variables, pv_config, {variable_name});
process_variables.push_back(std::move(per_process_variables));
}
- variable_ph = &process_variables[0][0].get();
- variable_u = &process_variables[1][0].get();
+ variable_u = &process_variables[0][0].get();
+ variable_ph = &process_variables[1][0].get();
}
DBUG("Associate displacement with process variable \'%s\'.",
@@ -94,7 +93,7 @@ std::unique_ptr<Process> createPhaseFieldProcess(
if (variable_ph->getNumberOfComponents() != 1)
{
OGS_FATAL(
- "Pressure 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());
@@ -118,11 +117,13 @@ std::unique_ptr<Process> createPhaseFieldProcess(
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()));
+ 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
{
@@ -191,10 +192,29 @@ std::unique_ptr<Process> createPhaseFieldProcess(
std::copy_n(b.data(), b.size(), specific_body_force.data());
}
+ auto const crack_scheme =
+ //! \ogs_file_param{prj__processes__process__PHASE_FIELD__hydro_crack_scheme}
+ config.getConfigParameterOptional<std::string>("hydro_crack_scheme");
+ if (crack_scheme &&
+ ((*crack_scheme != "propagating") && (*crack_scheme != "static")))
+ {
+ OGS_FATAL(
+ "hydro_crack_scheme must be \"propagating\" or \"static\" but "
+ "\"%s\" was given",
+ crack_scheme->c_str());
+ }
+
+ const bool propagating_crack =
+ (crack_scheme && (*crack_scheme == "propagating"));
+ const bool crack_pressure =
+ (crack_scheme &&
+ ((*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};
+ history_field, specific_body_force, propagating_crack,
+ crack_pressure};
SecondaryVariableCollection secondary_variables;
@@ -205,10 +225,10 @@ std::unique_ptr<Process> createPhaseFieldProcess(
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),
- use_monolithic_scheme);
+ 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),
+ use_monolithic_scheme);
}
template std::unique_ptr<Process> createPhaseFieldProcess<2>(
diff --git a/ProcessLib/PhaseField/LocalAssemblerInterface.h b/ProcessLib/PhaseField/LocalAssemblerInterface.h
index 72d970218ebd0f1f2a9ba77e4152147b6124c1d9..b0953d7ab36aa04e7d28d31ebee4b9fd76a9598c 100644
--- a/ProcessLib/PhaseField/LocalAssemblerInterface.h
+++ b/ProcessLib/PhaseField/LocalAssemblerInterface.h
@@ -18,82 +18,40 @@ namespace ProcessLib
{
namespace PhaseField
{
-
struct PhaseFieldLocalAssemblerInterface
: public ProcessLib::LocalAssemblerInterface,
public NumLib::ExtrapolatableElement
{
- virtual std::vector<double> const& getIntPtSigmaXX(
+ virtual std::vector<double> const& getIntPtSigma(
const double /*t*/,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
- virtual std::vector<double> const& getIntPtSigmaYY(
+ virtual std::vector<double> const& getIntPtEpsilon(
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 void computeCrackIntegral(
+ std::size_t mesh_item_id,
+ std::vector<
+ std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+ dof_tables,
+ GlobalVector const& x, double const t, double& crack_volume,
+ bool const use_monolithic_scheme,
+ CoupledSolutionsForStaggeredScheme const* const cpl_xs) = 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;
+ virtual void computeEnergy(
+ std::size_t mesh_item_id,
+ std::vector<
+ std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+ dof_tables,
+ GlobalVector const& x, double const t, double& elastic_energy,
+ double& surface_energy, double& pressure_work,
+ bool const use_monolithic_scheme,
+ CoupledSolutionsForStaggeredScheme const* const cpl_xs) = 0;
};
} // namespace PhaseField
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..e44074d184cc413d46ce14689df25e55cc978218
--- /dev/null
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -0,0 +1,585 @@
+/**
+ * \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
+ *
+ * \file PhaseFieldFEM-impl.h
+ * Created on January 8, 2018, 3:00 PM
+ */
+#pragma once
+
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
+
+namespace ProcessLib
+{
+namespace PhaseField
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>::
+ 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)
+{
+ using RhsVector =
+ typename ShapeMatricesType::template VectorType<phasefield_size +
+ displacement_size>;
+ using JacobianMatrix = typename ShapeMatricesType::template MatrixType<
+ phasefield_size + displacement_size,
+ phasefield_size + displacement_size>;
+
+ 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& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ 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
+ //
+
+ // Temporary storage used in the Kud and for potential reuse in the
+ // Kdu matrix.
+ auto const Kud_ip_contribution =
+ (B.transpose() * 2 * d_ip * sigma_tensile * N * w).eval();
+
+ Kud.noalias() += Kud_ip_contribution;
+
+ if (history_variable_prev < strain_energy_tensile)
+ {
+ history_variable = strain_energy_tensile;
+ Kdu.noalias() += Kud_ip_contribution.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];
+ double const d_dot_ip = N.dot(d_dot);
+
+ 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;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>::
+ assembleWithJacobianForStaggeredScheme(
+ double const t, 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_b_data, std::vector<double>& local_Jac_data,
+ LocalCoupledSolutions const& local_coupled_solutions)
+{
+ // For the equations with phase field.
+ if (local_coupled_solutions.process_id == 1)
+ {
+ assembleWithJacobianPhaseFiledEquations(
+ t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+ local_b_data, local_Jac_data, local_coupled_solutions);
+ return;
+ }
+
+ // For the equations with deformation
+ assembleWithJacobianForDeformationEquations(
+ t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+ local_b_data, local_Jac_data, local_coupled_solutions);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>::
+ assembleWithJacobianForDeformationEquations(
+ double const t, 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_b_data, std::vector<double>& local_Jac_data,
+ LocalCoupledSolutions const& local_coupled_solutions)
+{
+ using DeformationVector =
+ typename ShapeMatricesType::template VectorType<displacement_size>;
+ using DeformationMatrix =
+ typename ShapeMatricesType::template MatrixType<displacement_size,
+ displacement_size>;
+ using PhaseFieldVector =
+ typename ShapeMatricesType::template VectorType<phasefield_size>;
+
+ auto const& local_u = local_coupled_solutions.local_coupled_xs[0];
+ auto const& local_d = local_coupled_solutions.local_coupled_xs[1];
+ assert(local_u.size() == displacement_size);
+ assert(local_d.size() == phasefield_size);
+
+ auto const local_matrix_size = local_u.size();
+ auto d =
+ Eigen::Map<PhaseFieldVector const>(local_d.data(), phasefield_size);
+
+ auto u =
+ Eigen::Map<DeformationVector const>(local_u.data(), displacement_size);
+
+ auto local_Jac = MathLib::createZeroedMatrix<DeformationMatrix>(
+ local_Jac_data, local_matrix_size, local_matrix_size);
+
+ auto local_rhs = MathLib::createZeroedVector<DeformationVector>(
+ local_b_data, local_matrix_size);
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+ double const& dt = _process_data.dt;
+
+ auto local_pressure = 0.0;
+ if (_process_data.crack_pressure)
+ local_pressure = _process_data.unity_pressure;
+
+ int const n_integration_points = _integration_method.getNumberOfPoints();
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ 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& eps = _ip_data[ip].eps;
+ eps.noalias() = B * u;
+ 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;
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ degradation);
+
+ auto& sigma_real = _ip_data[ip].sigma_real;
+ auto const& C_tensile = _ip_data[ip].C_tensile;
+ auto const& C_compressive = _ip_data[ip].C_compressive;
+
+ 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.noalias() -=
+ (B.transpose() * sigma_real - N_u.transpose() * rho_sr * b -
+ local_pressure * N_u.transpose() * dNdx * d) *
+ w;
+
+ local_Jac.noalias() +=
+ B.transpose() * (degradation * C_tensile + C_compressive) * B * w;
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>::
+ assembleWithJacobianPhaseFiledEquations(
+ double const t, 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_b_data, std::vector<double>& local_Jac_data,
+ LocalCoupledSolutions const& local_coupled_solutions)
+{
+ using DeformationVector =
+ typename ShapeMatricesType::template VectorType<displacement_size>;
+ using PhaseFieldVector =
+ typename ShapeMatricesType::template VectorType<phasefield_size>;
+ using PhaseFieldMatrix =
+ typename ShapeMatricesType::template MatrixType<phasefield_size,
+ phasefield_size>;
+
+ auto const& local_u = local_coupled_solutions.local_coupled_xs[0];
+ auto const& local_d = local_coupled_solutions.local_coupled_xs[1];
+ assert(local_u.size() == displacement_size);
+ assert(local_d.size() == phasefield_size);
+
+ auto const local_matrix_size = local_d.size();
+ auto d =
+ Eigen::Map<PhaseFieldVector const>(local_d.data(), phasefield_size);
+ auto u =
+ Eigen::Map<DeformationVector const>(local_u.data(), displacement_size);
+
+ auto local_Jac = MathLib::createZeroedMatrix<PhaseFieldMatrix>(
+ local_Jac_data, local_matrix_size, local_matrix_size);
+ auto local_rhs = MathLib::createZeroedVector<PhaseFieldVector>(
+ local_b_data, local_matrix_size);
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+ double const& dt = _process_data.dt;
+
+ auto local_pressure = 0.0;
+ if (_process_data.crack_pressure)
+ local_pressure = _process_data.unity_pressure;
+ else if (_process_data.propagating_crack)
+ local_pressure = _process_data.pressure;
+
+ int const n_integration_points = _integration_method.getNumberOfPoints();
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ double const gc = _process_data.crack_resistance(t, x_position)[0];
+ double const ls = _process_data.crack_length_scale(t, x_position)[0];
+
+ // for propagating crack, u is rescaled.
+ if (_process_data.propagating_crack)
+ {
+ 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;
+ 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& eps = _ip_data[ip].eps;
+ eps.noalias() = B * u;
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ degradation);
+ }
+
+ auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
+
+ auto& ip_data = _ip_data[ip];
+ ip_data.strain_energy_tensile = strain_energy_tensile;
+
+ 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;
+
+ local_Jac.noalias() +=
+ (2 * N.transpose() * N * strain_energy_tensile +
+ gc * (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls)) *
+ w;
+
+ local_rhs.noalias() -=
+ (N.transpose() * N * d * 2 * strain_energy_tensile +
+ gc * ((N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) * d -
+ N.transpose() / ls) -
+ local_pressure * dNdx.transpose() * N_u * u) *
+ w;
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>::
+ computeCrackIntegral(std::size_t mesh_item_id,
+ std::vector<std::reference_wrapper<
+ NumLib::LocalToGlobalIndexMap>> const& dof_tables,
+ GlobalVector const& /*x*/, double const /*t*/,
+ double& crack_volume,
+ bool const /*use_monolithic_scheme*/,
+ CoupledSolutionsForStaggeredScheme const* const cpl_xs)
+{
+ assert(cpl_xs != nullptr);
+
+ std::vector<std::vector<GlobalIndexType>> indices_of_processes;
+ indices_of_processes.reserve(dof_tables.size());
+ std::transform(dof_tables.begin(), dof_tables.end(),
+ std::back_inserter(indices_of_processes),
+ [&](NumLib::LocalToGlobalIndexMap const& dof_table) {
+ return NumLib::getIndices(mesh_item_id, dof_table);
+ });
+
+ auto local_coupled_xs =
+ getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
+ assert(local_coupled_xs.size() == 2);
+
+ auto const& local_u = local_coupled_xs[0];
+ auto const& local_d = local_coupled_xs[1];
+
+ assert(local_u.size() == displacement_size);
+ assert(local_d.size() == phasefield_size);
+
+ auto d = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<phasefield_size> const>(
+ local_d.data(), phasefield_size);
+
+ auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ displacement_size> const>(local_u.data(), displacement_size);
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ int const n_integration_points = _integration_method.getNumberOfPoints();
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ 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;
+
+ crack_volume += (N_u * u).dot(dNdx * d) * w;
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>::
+ computeEnergy(std::size_t mesh_item_id,
+ std::vector<std::reference_wrapper<
+ NumLib::LocalToGlobalIndexMap>> const& dof_tables,
+ GlobalVector const& /*x*/, double const t,
+ double& elastic_energy, double& surface_energy,
+ double& pressure_work, bool const /*use_monolithic_scheme*/,
+ CoupledSolutionsForStaggeredScheme const* const cpl_xs)
+{
+ assert(cpl_xs != nullptr);
+
+ std::vector<std::vector<GlobalIndexType>> indices_of_processes;
+ indices_of_processes.reserve(dof_tables.size());
+ std::transform(dof_tables.begin(), dof_tables.end(),
+ std::back_inserter(indices_of_processes),
+ [&](NumLib::LocalToGlobalIndexMap const& dof_table) {
+ return NumLib::getIndices(mesh_item_id, dof_table);
+ });
+
+ auto local_coupled_xs =
+ getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
+ assert(local_coupled_xs.size() == 2);
+
+ auto const& local_u = local_coupled_xs[0];
+ auto const& local_d = local_coupled_xs[1];
+
+ assert(local_u.size() == displacement_size);
+ assert(local_d.size() == phasefield_size);
+
+ auto d = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<phasefield_size> const>(
+ local_d.data(), phasefield_size);
+
+ auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ displacement_size> const>(local_u.data(), displacement_size);
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ int const n_integration_points = _integration_method.getNumberOfPoints();
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+ double const d_ip = N.dot(d);
+ auto pressure_ip = _process_data.pressure;
+ auto u_corrected = pressure_ip * u;
+ 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::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;
+
+ elastic_energy += _ip_data[ip].elastic_energy * w;
+
+ surface_energy +=
+ 0.5 * gc *
+ ((1 - d_ip) * (1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
+ w;
+
+ if (_process_data.crack_pressure)
+ pressure_work +=
+ pressure_ip * (N_u * u_corrected).dot(dNdx * d) * w;
+ }
+}
+} // namespace PhaseField
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index c9bde13feaba2385cff617d1b59af6ecdcf7bce7..c168e0f8089fe25feed65af7c57f6f056d6fbc65 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -46,7 +46,7 @@ struct IntegrationPointData final
typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
sigma_real_prev, sigma_real;
- double strain_energy_tensile;
+ double strain_energy_tensile, elastic_energy;
MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
@@ -81,7 +81,7 @@ struct IntegrationPointData final
.calculateDegradedStress(t, x_position, eps, strain_energy_tensile,
sigma_tensile, sigma_compressive,
C_tensile, C_compressive, sigma_real,
- degradation);
+ degradation, elastic_energy);
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
@@ -101,18 +101,12 @@ 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;
@@ -166,6 +160,8 @@ public:
ip_data.history_variable_prev =
_process_data.history_field(0, x_position)[0];
ip_data.sigma_real.setZero(kelvin_vector_size);
+ ip_data.strain_energy_tensile = 0.0;
+ ip_data.elastic_energy = 0.0;
ip_data.N = shape_matrices[ip].N;
ip_data.dNdx = shape_matrices[ip].dNdx;
@@ -191,180 +187,14 @@ public:
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);
+ std::vector<double>& local_Jac_data) override;
- 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& N = _ip_data[ip].N;
- auto const& dNdx = _ip_data[ip].dNdx;
-
- 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
- //
-
- // Temporary storage used in the Kud and for potential reuse in the
- // Kdu matrix.
- auto const Kud_ip_contribution =
- (B.transpose() * 2 * d_ip * sigma_tensile * N * w).eval();
-
- Kud.noalias() += Kud_ip_contribution;
-
- if (history_variable_prev < strain_energy_tensile)
- {
- history_variable = strain_energy_tensile;
- Kdu.noalias() += Kud_ip_contribution.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];
- double const d_dot_ip = N.dot(d_dot);
-
- 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 assembleWithJacobianForStaggeredScheme(
+ double const t, 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_b_data, std::vector<double>& local_Jac_data,
+ LocalCoupledSolutions const& local_coupled_solutions) override;
void preTimestepConcrete(std::vector<double> const& /*local_x*/,
double const /*t*/,
@@ -379,6 +209,25 @@ public:
}
}
+ void computeCrackIntegral(
+ std::size_t mesh_item_id,
+ std::vector<
+ std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+ dof_tables,
+ GlobalVector const& x, double const t, double& crack_volume,
+ bool const use_monolithic_scheme,
+ CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
+
+ void computeEnergy(
+ std::size_t mesh_item_id,
+ std::vector<
+ std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+ dof_tables,
+ GlobalVector const& x, double const t, double& elastic_energy,
+ double& surface_energy, double& pressure_work,
+ bool const use_monolithic_scheme,
+ CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
+
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -388,152 +237,72 @@ public:
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(
+private:
+ std::vector<double> const& getIntPtSigma(
const double /*t*/,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
- return getIntPtEpsilon(cache, 1);
- }
+ static const int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ auto const num_intpts = _ip_data.size();
- 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);
- }
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, num_intpts);
- 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);
- }
+ for (unsigned ip = 0; ip < num_intpts; ++ip)
+ {
+ auto const& sigma = _ip_data[ip].sigma_real;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma);
+ }
- 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);
+ return cache;
}
- std::vector<double> const& getIntPtEpsilonXZ(
+ virtual std::vector<double> const& getIntPtEpsilon(
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);
- }
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ auto const num_intpts = _ip_data.size();
-private:
- std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
- std::size_t const component) const
- {
cache.clear();
- cache.reserve(_ip_data.size());
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, num_intpts);
- for (auto const& ip_data : _ip_data)
+ for (unsigned ip = 0; ip < num_intpts; ++ip)
{
- 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));
+ auto const& eps = _ip_data[ip].eps;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
}
return cache;
}
- std::vector<double> const& getIntPtEpsilon(
- std::vector<double>& cache, std::size_t const component) const
- {
- cache.clear();
- cache.reserve(_ip_data.size());
+ void assembleWithJacobianPhaseFiledEquations(
+ double const t, 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_b_data, std::vector<double>& local_Jac_data,
+ LocalCoupledSolutions const& local_coupled_solutions);
- 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;
- }
+ void assembleWithJacobianForDeformationEquations(
+ double const t, 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_b_data, std::vector<double>& local_Jac_data,
+ LocalCoupledSolutions const& local_coupled_solutions);
PhaseFieldProcessData<DisplacementDim>& _process_data;
@@ -557,3 +326,5 @@ private:
} // namespace PhaseField
} // namespace ProcessLib
+
+#include "PhaseFieldFEM-impl.h"
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
index eedd3abc08af808b1f93005b477cc099cb09354a..1f49952126b504ec14a73f1c2fe86bf13be79838 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
@@ -11,6 +11,8 @@
#include <cassert>
+#include "NumLib/DOF/ComputeSparsityPattern.h"
+
#include "ProcessLib/Process.h"
#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
@@ -38,6 +40,12 @@ PhaseFieldProcess<DisplacementDim>::PhaseFieldProcess(
use_monolithic_scheme),
_process_data(std::move(process_data))
{
+ // Disable nodal forces for monolithic scheme.
+ if (!_use_monolithic_scheme)
+ {
+ _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+ }
}
template <int DisplacementDim>
@@ -47,18 +55,49 @@ bool PhaseFieldProcess<DisplacementDim>::isLinear() const
}
template <int DisplacementDim>
-void PhaseFieldProcess<DisplacementDim>::initializeConcreteProcess(
- NumLib::LocalToGlobalIndexMap const& dof_table,
- MeshLib::Mesh const& mesh,
- unsigned const integration_order)
+MathLib::MatrixSpecifications
+PhaseFieldProcess<DisplacementDim>::getMatrixSpecifications(
+ const int process_id) const
{
- ProcessLib::SmallDeformation::createLocalAssemblers<
- DisplacementDim, PhaseFieldLocalAssembler>(
- mesh.getElements(), dof_table, _local_assemblers,
- mesh.isAxiallySymmetric(), integration_order, _process_data);
+ // For the monolithic scheme or the M process (deformation) in the staggered
+ // scheme.
+ if (_use_monolithic_scheme || process_id == 0)
+ {
+ auto const& l = *_local_to_global_index_map;
+ return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+ &l.getGhostIndices(), &this->_sparsity_pattern};
+ }
+
+ // For staggered scheme and phase field process.
+ auto const& l = *_local_to_global_index_map_single_component;
+ return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+ &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
+}
+
+template <int DisplacementDim>
+NumLib::LocalToGlobalIndexMap const&
+PhaseFieldProcess<DisplacementDim>::getDOFTable(const int process_id) const
+{
+ // If monolithic scheme is used or the equation of deformation is solved in
+ // the staggered scheme.
+ if (_use_monolithic_scheme || process_id == 0)
+ {
+ return *_local_to_global_index_map;
+ }
+
+ // For the equation of phasefield
+ return *_local_to_global_index_map_single_component;
+}
+
+template <int DisplacementDim>
+void PhaseFieldProcess<DisplacementDim>::constructDofTable()
+{
+ // Create single component dof in every of the mesh's nodes.
+ _mesh_subset_all_nodes =
+ std::make_unique<MeshLib::MeshSubset>(_mesh, &_mesh.getNodes());
// TODO move the two data members somewhere else.
- // for extrapolation of secondary variables
+ // for extrapolation of secondary variables of stress or strain
std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
all_mesh_subsets_single_component.emplace_back(
_mesh_subset_all_nodes.get());
@@ -68,71 +107,94 @@ void PhaseFieldProcess<DisplacementDim>::initializeConcreteProcess(
// by location order is needed for output
NumLib::ComponentOrder::BY_LOCATION);
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_xx",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXX));
+ assert(_local_to_global_index_map_single_component);
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_yy",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaYY));
+ if (_use_monolithic_scheme)
+ {
+ const int process_id = 0; // Only one process in the monolithic scheme.
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+ std::generate_n(
+ std::back_inserter(all_mesh_subsets),
+ getProcessVariables(process_id)[1].get().getNumberOfComponents() +
+ 1,
+ [&]() {
+ return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
+ });
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_zz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaZZ));
+ std::vector<int> const vec_n_components{1, DisplacementDim};
+ _local_to_global_index_map =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), vec_n_components,
+ NumLib::ComponentOrder::BY_LOCATION);
+ assert(_local_to_global_index_map);
+ }
+ else
+ {
+ // For displacement equation.
+ const int process_id = 0;
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+ std::generate_n(
+ std::back_inserter(all_mesh_subsets),
+ getProcessVariables(process_id)[0].get().getNumberOfComponents(),
+ [&]() {
+ return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
+ });
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_xy",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXY));
+ std::vector<int> const vec_n_components{DisplacementDim};
+ _local_to_global_index_map =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), vec_n_components,
+ NumLib::ComponentOrder::BY_LOCATION);
- if (DisplacementDim == 3)
- {
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_xz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXZ));
-
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_yz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaYZ));
+ // For phase field equation.
+ _sparsity_pattern_with_single_component =
+ NumLib::computeSparsityPattern(
+ *_local_to_global_index_map_single_component, _mesh);
}
+}
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xx",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonXX));
-
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yy",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonYY));
+template <int DisplacementDim>
+void PhaseFieldProcess<DisplacementDim>::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ ProcessLib::SmallDeformation::createLocalAssemblers<
+ DisplacementDim, PhaseFieldLocalAssembler>(
+ mesh.getElements(), dof_table, _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_zz",
+ "sigma",
makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonZZ));
+ &LocalAssemblerInterface::getIntPtSigma));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xy",
+ "epsilon",
makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonXY));
+ &LocalAssemblerInterface::getIntPtEpsilon));
+}
- if (DisplacementDim == 3)
+template <int DisplacementDim>
+void PhaseFieldProcess<DisplacementDim>::initializeBoundaryConditions()
+{
+ if (_use_monolithic_scheme)
{
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonYZ));
-
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonXZ));
+ const int process_id_of_phasefieldmechanics = 0;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map, process_id_of_phasefieldmechanics);
+ return;
}
+
+ // Staggered scheme:
+ // for the equations of deformation.
+ const int mechanical_process_id = 0;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map, mechanical_process_id);
+ // for the phase field
+ const int phasefield_process_id = 1;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map_single_component, phasefield_process_id);
}
template <int DisplacementDim>
@@ -143,7 +205,33 @@ void PhaseFieldProcess<DisplacementDim>::assembleConcreteProcess(
DBUG("Assemble PhaseFieldProcess.");
std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
- dof_tables = {std::ref(*_local_to_global_index_map)};
+ dof_tables;
+
+ // For the monolithic scheme
+ if (_use_monolithic_scheme)
+ {
+ DBUG("Assemble PhaseFieldProcess for the monolithic scheme.");
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ }
+ else
+ {
+ // For the staggered scheme
+ if (_coupled_solutions->process_id == 1)
+ {
+ DBUG(
+ "Assemble the equations of phase field in "
+ "PhaseFieldProcess for the staggered scheme.");
+ }
+ else
+ {
+ DBUG(
+ "Assemble the equations of deformation in "
+ "PhaseFieldProcess for the staggered scheme.");
+ }
+ dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ }
+
// Call global assembler for each local assembly item.
GlobalExecutor::executeMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
@@ -156,42 +244,144 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
GlobalVector& b, GlobalMatrix& Jac)
{
- // DBUG("AssembleJacobian PhaseFieldProcess.");
-
std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
- dof_tables = {std::ref(*_local_to_global_index_map)};
+ dof_tables;
+
+ // For the monolithic scheme
+ if (_use_monolithic_scheme)
+ {
+ DBUG("AssembleJacobian PhaseFieldProcess for the monolithic scheme.");
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ }
+ else
+ {
+ // For the staggered scheme
+ if (_coupled_solutions->process_id == 1)
+ {
+ DBUG(
+ "Assemble the Jacobian equations of phase field in "
+ "PhaseFieldProcess for the staggered scheme.");
+ }
+ else
+ {
+ DBUG(
+ "Assemble the Jacobian equations of deformation in "
+ "PhaseFieldProcess for the staggered scheme.");
+ }
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+ }
// Call global assembler for each local assembly item.
+
GlobalExecutor::executeMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
- _local_assemblers, dof_tables, t, x, xdot, dxdot_dx,
- dx_dx, M, K, b, Jac, _coupled_solutions);
+ _local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+ Jac, _coupled_solutions);
+
+ if (!_use_monolithic_scheme && (_coupled_solutions->process_id == 0))
+ {
+ b.copyValues(*_nodal_forces);
+ std::transform(_nodal_forces->begin(), _nodal_forces->end(),
+ _nodal_forces->begin(), [](double val) { return -val; });
+ }
}
template <int DisplacementDim>
void PhaseFieldProcess<DisplacementDim>::preTimestepConcreteProcess(
GlobalVector const& x, double const t, double const dt,
- const int /*process_id*/)
+ const int process_id)
{
- DBUG("PreTimestep PhaseFieldProcess.");
+ DBUG("PreTimestep PhaseFieldProcess %d.", process_id);
_process_data.dt = dt;
_process_data.t = t;
+ _process_data.injected_volume = _process_data.t;
GlobalExecutor::executeMemberOnDereferenced(
&LocalAssemblerInterface::preTimestep, _local_assemblers,
- *_local_to_global_index_map, x, t, dt);
+ getDOFTable(process_id), x, t, dt);
}
template <int DisplacementDim>
void PhaseFieldProcess<DisplacementDim>::postTimestepConcreteProcess(
- GlobalVector const& x, int const /*process_id*/)
+ GlobalVector const& x, int const process_id)
{
- DBUG("PostTimestep PhaseFieldProcess.");
+ if (isPhaseFieldProcess(process_id))
+ {
+ DBUG("PostTimestep PhaseFieldProcess.");
- GlobalExecutor::executeMemberOnDereferenced(
- &LocalAssemblerInterface::postTimestep, _local_assemblers,
- *_local_to_global_index_map, x);
+ _process_data.elastic_energy = 0.0;
+ _process_data.surface_energy = 0.0;
+ _process_data.pressure_work = 0.0;
+
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_tables;
+
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::computeEnergy, _local_assemblers,
+ dof_tables, x, _process_data.t, _process_data.elastic_energy,
+ _process_data.surface_energy, _process_data.pressure_work,
+ _use_monolithic_scheme, _coupled_solutions);
+
+ INFO("Elastic energy: %g Surface energy: %g Pressure work: %g ",
+ _process_data.elastic_energy, _process_data.surface_energy,
+ _process_data.pressure_work);
+ }
}
+template <int DisplacementDim>
+void PhaseFieldProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
+ GlobalVector const& x, const double t, const int process_id)
+{
+ if (_use_monolithic_scheme)
+ return;
+
+ _process_data.crack_volume = 0.0;
+
+ if (!isPhaseFieldProcess(process_id))
+ {
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_tables;
+
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+
+ DBUG("PostNonLinearSolver crack volume computation.");
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::computeCrackIntegral, _local_assemblers,
+ dof_tables, x, t, _process_data.crack_volume,
+ _use_monolithic_scheme, _coupled_solutions);
+
+ INFO("Integral of crack: %g", _process_data.crack_volume);
+
+ if (_process_data.propagating_crack)
+ {
+ _process_data.pressure_old = _process_data.pressure;
+ _process_data.pressure =
+ _process_data.injected_volume / _process_data.crack_volume;
+ _process_data.pressure_error =
+ abs(_process_data.pressure_old - _process_data.pressure) /
+ _process_data.pressure;
+ INFO("Internal pressure: %g and Pressure error: %.4e",
+ _process_data.pressure, _process_data.pressure_error);
+ auto& u = _coupled_solutions->coupled_xs[0].get();
+ MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
+ _process_data.pressure);
+ }
+ }
+ else
+ {
+ if (_process_data.propagating_crack)
+ {
+ auto& u = _coupled_solutions->coupled_xs[0].get();
+ MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
+ 1 / _process_data.pressure);
+ }
+ }
+}
} // namespace PhaseField
} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.h b/ProcessLib/PhaseField/PhaseFieldProcess.h
index efdc9d2b5954ae0b7a5b8be02edc726979a1c798..b0bdfaf9995d947ffa1a577c784c777dbd97202f 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.h
@@ -18,7 +18,6 @@ namespace ProcessLib
{
namespace PhaseField
{
-
/**
* \brief A class to simulate mechanical fracturing process
* using phase-field approach in solids described by
@@ -69,30 +68,44 @@ public:
bool isLinear() const override;
//! @}
+ MathLib::MatrixSpecifications getMatrixSpecifications(
+ const int process_id) const override;
+
+ NumLib::LocalToGlobalIndexMap const& getDOFTable(
+ const int process_id) const override;
+
private:
using LocalAssemblerInterface = PhaseFieldLocalAssemblerInterface;
+ void constructDofTable() override;
+
+ void initializeBoundaryConditions() override;
+
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh,
unsigned const integration_order) override;
- void assembleConcreteProcess(
- const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b) override;
+ void assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) override;
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) override;
- void preTimestepConcreteProcess(
- GlobalVector const& x, double const t, double const dt,
- const int process_id) override;
+ void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+ double const dt,
+ const int process_id) override;
void postTimestepConcreteProcess(GlobalVector const& x,
int const process_id) override;
+ void postNonLinearSolverConcreteProcess(GlobalVector const& x,
+ const double t,
+ int const process_id) override;
+
private:
PhaseFieldProcessData<DisplacementDim> _process_data;
@@ -100,6 +113,20 @@ private:
std::unique_ptr<NumLib::LocalToGlobalIndexMap>
_local_to_global_index_map_single_component;
+
+ MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
+
+ /// Sparsity pattern for the phase field equation, and it is initialized
+ /// only if the staggered scheme is used.
+ GlobalSparsityPattern _sparsity_pattern_with_single_component;
+
+ /// Check whether the process represented by \c process_id is/has
+ /// mechanical process. In the present implementation, the mechanical
+ /// process has process_id == 0 in the staggered scheme.
+ bool isPhaseFieldProcess(int const process_id) const
+ {
+ return !_use_monolithic_scheme && process_id == 1;
+ }
};
extern template class PhaseFieldProcess<2>;
diff --git a/ProcessLib/PhaseField/PhaseFieldProcessData.h b/ProcessLib/PhaseField/PhaseFieldProcessData.h
index 7184f64298b4c091e92e005d3fc42b423484360e..75fdace07173ad99b6c67235b9c3a0bae855f470 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcessData.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcessData.h
@@ -21,7 +21,7 @@ namespace Solids
template <int DisplacementDim>
struct MechanicsBase;
}
-}
+} // namespace MaterialLib
namespace ProcessLib
{
template <typename T>
@@ -41,7 +41,8 @@ struct PhaseFieldProcessData
Parameter<double> const& kinetic_coefficient_,
Parameter<double> const& solid_density_,
Parameter<double>& history_field_,
- Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_)
+ Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_,
+ bool propagating_crack_, bool crack_pressure_)
: material{std::move(material_)},
residual_stiffness(residual_stiffness_),
crack_resistance(crack_resistance_),
@@ -49,7 +50,9 @@ struct PhaseFieldProcessData
kinetic_coefficient(kinetic_coefficient_),
solid_density(solid_density_),
history_field(history_field_),
- specific_body_force(specific_body_force_)
+ specific_body_force(specific_body_force_),
+ propagating_crack(propagating_crack_),
+ crack_pressure(crack_pressure_)
{
}
@@ -63,7 +66,9 @@ struct PhaseFieldProcessData
history_field(other.history_field),
specific_body_force(other.specific_body_force),
dt(other.dt),
- t(other.t)
+ t(other.t),
+ propagating_crack(other.propagating_crack),
+ crack_pressure(other.crack_pressure)
{
}
@@ -87,6 +92,17 @@ struct PhaseFieldProcessData
Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
double dt = 0.0;
double t = 0.0;
+ double const unity_pressure = 1.0;
+ double pressure = 0.0;
+ double pressure_old = 0.0;
+ double pressure_error = 0.0;
+ double injected_volume = 0.0;
+ double crack_volume = 0.0;
+ double elastic_energy = 0.0;
+ double surface_energy = 0.0;
+ double pressure_work = 0.0;
+ bool propagating_crack = false;
+ bool crack_pressure = false;
};
} // namespace PhaseField
diff --git a/ProcessLib/PhaseField/Tests.cmake b/ProcessLib/PhaseField/Tests.cmake
index 29c3a780e48a70d7e00e759b590cc699177a4f0d..b7a8a39025b91538779cdcb582d6be52065d72f8 100644
--- a/ProcessLib/PhaseField/Tests.cmake
+++ b/ProcessLib/PhaseField/Tests.cmake
@@ -10,3 +10,56 @@ AddTest(
expected_cube_1e0_pcs_0_ts_10000_t_1.000000.vtu cube_1e0_pcs_0_ts_10000_t_1.000000.vtu displacement displacement 1e-6 1e-6
expected_cube_1e0_pcs_0_ts_10000_t_1.000000.vtu cube_1e0_pcs_0_ts_10000_t_1.000000.vtu phasefield phasefield 1e-6 1e-6
)
+
+AddTest(
+ NAME PhaseField_2D_StaticHydraulicFracture
+ PATH PhaseField
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS 2D_static.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ DIFF_DATA
+ expected_2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu 2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu displacement displacement 1e-15 1e-15
+ expected_2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu 2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu phasefield phasefield 1e-15 1e-15
+ )
+
+AddTest(
+ NAME PhaseField_3D_beam_LARGE
+ PATH PhaseField
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS beam3d_stag_1pcs.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ DIFF_DATA
+ expected_beam3_stag1pcsAT2_pcs_1_ts_10_t_1.000000.vtu beam3_stag1pcsAT2_pcs_1_ts_10_t_1.000000.vtu displacement displacement 1e-5 0
+ expected_beam3_stag1pcsAT2_pcs_1_ts_10_t_1.000000.vtu beam3_stag1pcsAT2_pcs_1_ts_10_t_1.000000.vtu phasefield phasefield 1e-6 0
+ )
+
+
+AddTest(
+ NAME LARGE_PhaseField_Staggered_square_line
+ PATH PhaseField
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_line_h_400.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ DIFF_DATA
+ square_line_h_400_pcs_0_ts_400000_t_400.000000.vtu square_line_h_400_pcs_0_ts_400000_t_400.000000.vtu displacement displacement 1e-16 0
+ square_line_h_400_pcs_0_ts_400000_t_400.000000.vtu square_line_h_400_pcs_0_ts_400000_t_400.000000.vtu phasefield phasefield 1e-16 0
+ )
+
+AddTest(
+ NAME LARGE_PhaseField_Staggered_square_shear
+ PATH PhaseField
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_shear_h_400.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ DIFF_DATA
+ square_shear_h_400_pcs_0_ts_30000_t_300.000000.vtu square_shear_h_400_pcs_0_ts_30000_t_300.000000.vtu displacement displacement 1e-16 0
+ square_shear_h_400_pcs_0_ts_30000_t_300.000000.vtu square_shear_h_400_pcs_0_ts_30000_t_300.000000.vtu phasefield phasefield 1e-16 0
+ )
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 4647714e169eeeb86e119c127b63cbc16e60408a..9e9d9ab2e7e2b2c6edfcfcc10f40f6fadc56c780 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -346,6 +346,8 @@ void Process::preTimestep(GlobalVector const& x, const double t,
MathLib::LinAlg::setLocalAccessibleVector(x);
preTimestepConcreteProcess(x, t, delta_t, process_id);
+
+ _boundary_conditions[process_id].preTimestep(t, x);
}
void Process::postTimestep(GlobalVector const& x, int const process_id)
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index 93474a224f8c1238941233ec4bb69fb3ee696933..48d45cd23714ee5b1fe49743e6610847ae032155 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -165,11 +165,15 @@ ProcessVariable::createBoundaryConditions(
std::vector<std::unique_ptr<ParameterBase>> const& parameters)
{
std::vector<std::unique_ptr<BoundaryCondition>> bcs;
+ bcs.reserve(_bc_configs.size());
for (auto& config : _bc_configs)
- bcs.emplace_back(_bc_builder->createBoundaryCondition(
+ {
+ auto bc = _bc_builder->createBoundaryCondition(
config, dof_table, _mesh, variable_id, integration_order,
- _shapefunction_order, parameters));
+ _shapefunction_order, parameters);
+ bcs.push_back(std::move(bc));
+ }
return bcs;
}
diff --git a/Tests/Data/PhaseField/2D_static.prj b/Tests/Data/PhaseField/2D_static.prj
new file mode 100644
index 0000000000000000000000000000000000000000..d848c5d108aa41c773cb4bc6bed248f52c6a6429
--- /dev/null
+++ b/Tests/Data/PhaseField/2D_static.prj
@@ -0,0 +1,307 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>cracked_medium_ic.vtu</mesh>
+ <geometry>bar.gml</geometry>
+ <processes>
+ <process>
+ <name>PhaseField</name>
+ <type>PHASE_FIELD</type>
+ <coupling_scheme>staggered</coupling_scheme>
+ <hydro_crack_scheme>static</hydro_crack_scheme>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>LinearElasticIsotropic</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ </constitutive_relation>
+ <phasefield_parameters>
+ <residual_stiffness>k</residual_stiffness>
+ <crack_resistance>gc</crack_resistance>
+ <crack_length_scale>ls</crack_length_scale>
+ <kinetic_coefficient>M</kinetic_coefficient>
+ <history_field>H</history_field>
+ </phasefield_parameters>
+ <solid_density>rho_sr</solid_density>
+ <process_variables>
+ <phasefield>phasefield</phasefield>
+ <displacement>displacement</displacement>
+ </process_variables>
+ <specific_body_force>0 0 </specific_body_force>
+ </process>
+ </processes>
+ <time_loop>
+ <global_process_coupling>
+ <max_iter> 1000 </max_iter>
+ <convergence_criteria>
+ <!-- convergence criterion for the first process -->
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>INFINITY_N</norm_type>
+ <abstol>1.e-1</abstol>
+ <reltol>1.e-1</reltol>
+ </convergence_criterion>
+ <!-- convergence criterion for the second process -->
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>INFINITY_N</norm_type>
+ <abstol>1.e-4</abstol>
+ <reltol>1.e-10</reltol>
+ </convergence_criterion>
+ </convergence_criteria>
+ </global_process_coupling>
+ <processes>
+ <process ref="PhaseField">
+ <nonlinear_solver>basic_newton_u</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1.e-14</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable>displacement</variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>1</t_end>
+ <timesteps>
+ <pair>
+ <repeat>2</repeat>
+ <delta_t>1.0</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ <process ref="PhaseField">
+ <nonlinear_solver>basic_newton_d</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1.e-14</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable>phasefield</variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>1</t_end>
+ <timesteps>
+ <pair>
+ <repeat>2</repeat>
+ <delta_t>1.0</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>2D_StaticCrack</prefix>
+ <timesteps>
+ <pair>
+ <repeat>10</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ <pair>
+ <repeat>10</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ <!--<pair><repeat>50</repeat><each_steps>100</each_steps></pair>-->
+ </timesteps>
+ </output>
+ </time_loop>
+ <parameters>
+ <!-- Mechanics -->
+ <parameter>
+ <name>E</name>
+ <type>Constant</type>
+ <value>1.0</value>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>Constant</type>
+ <value>0.15</value>
+ </parameter>
+ <parameter>
+ <name>k</name>
+ <type>Constant</type>
+ <value>1e-8</value>
+ </parameter>
+ <parameter>
+ <name>gc</name>
+ <type>Constant</type>
+ <value>1.0</value>
+ </parameter>
+ <parameter>
+ <name>ls</name>
+ <type>Constant</type>
+ <value>0.02</value>
+ </parameter>
+ <parameter>
+ <name>H</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>phasefield_ic</name>
+ <type>MeshNode</type>
+ <field_name>phase-field</field_name>
+ <!--type>Constant</type>
+ <value>1</value>-->
+ </parameter>
+ <parameter>
+ <name>phasefield_bc</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>dirichlet0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>M</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>Dirichlet_load</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>Dirichlet_top_y</name>
+ <type>CurveScaled</type>
+ <curve>Dirichlet_top_time</curve>
+ <parameter>Dirichlet_load</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>Dirichlet_top_time</name>
+ <coords>0 1</coords>
+ <values>0 1.0</values>
+ </curve>
+ </curves>
+ <process_variables>
+ <process_variable>
+ <name>displacement</name>
+ <components>2</components>
+ <order>1</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>phasefield</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>phasefield_ic</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>phasefield_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>phasefield_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>right</geometry>
+ <type>PhaseFieldIrreversibleDamageOracleBoundaryCondition</type>
+ <component>0</component>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton_d</name>
+ <type>Newton</type>
+ <max_iter>50</max_iter>
+ <linear_solver>linear_solver_d</linear_solver>
+ </nonlinear_solver>
+ <nonlinear_solver>
+ <name>basic_newton_u</name>
+ <type>Newton</type>
+ <max_iter>50</max_iter>
+ <linear_solver>linear_solver_u</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>linear_solver_d</name>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ </linear_solver>
+ <linear_solver>
+ <name>linear_solver_u</name>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/PhaseField/bar.gml b/Tests/Data/PhaseField/bar.gml
new file mode 100644
index 0000000000000000000000000000000000000000..d226fe62c21fdde3fde76e6d94cd2b29f34016c2
--- /dev/null
+++ b/Tests/Data/PhaseField/bar.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:19c58dced67204b13be71364487e0f0591191151b6162d19bc29db839471c1fa
+size 910
diff --git a/Tests/Data/PhaseField/beam3d.prj b/Tests/Data/PhaseField/beam3d.prj
index ee6df112f292188cb8fa905f4c443e614e444169..e41182aa564c32769a806bf5f5ee3a3864508a36 100644
--- a/Tests/Data/PhaseField/beam3d.prj
+++ b/Tests/Data/PhaseField/beam3d.prj
@@ -25,14 +25,8 @@
<phasefield>phasefield</phasefield>
</process_variables>
<secondary_variables>
- <secondary_variable type="static" internal_name="sigma_xx" output_name="sigma_xx"/>
- <secondary_variable type="static" internal_name="sigma_yy" output_name="sigma_yy"/>
- <secondary_variable type="static" internal_name="sigma_zz" output_name="sigma_zz"/>
- <secondary_variable type="static" internal_name="sigma_xy" output_name="sigma_xy"/>
- <secondary_variable type="static" internal_name="epsilon_xx" output_name="epsilon_xx"/>
- <secondary_variable type="static" internal_name="epsilon_yy" output_name="epsilon_yy"/>
- <secondary_variable type="static" internal_name="epsilon_zz" output_name="epsilon_zz"/>
- <secondary_variable type="static" internal_name="epsilon_xy" output_name="epsilon_xy"/>
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
<specific_body_force>0 0 0</specific_body_force>
</process>
@@ -53,14 +47,8 @@
<variables>
<variable>displacement</variable>
<variable>phasefield</variable>
- <variable>sigma_xx</variable>
- <variable>sigma_yy</variable>
- <variable>sigma_zz</variable>
- <variable>sigma_xy</variable>
- <variable>epsilon_xx</variable>
- <variable>epsilon_yy</variable>
- <variable>epsilon_zz</variable>
- <variable>epsilon_xy</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
</variables>
</output>
<time_stepping>
diff --git a/Tests/Data/PhaseField/beam3d_stag_1pcs.prj b/Tests/Data/PhaseField/beam3d_stag_1pcs.prj
new file mode 100644
index 0000000000000000000000000000000000000000..2d71e978074adec4655b7c873e065d8573e157c2
--- /dev/null
+++ b/Tests/Data/PhaseField/beam3d_stag_1pcs.prj
@@ -0,0 +1,322 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>beam3d.vtu</mesh>
+ <geometry>beam3d.gml</geometry>
+ <processes>
+ <process>
+ <name>PhaseField</name>
+ <type>PHASE_FIELD</type>
+ <coupling_scheme>staggered</coupling_scheme>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>LinearElasticIsotropic</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ </constitutive_relation>
+ <phasefield_parameters>
+ <residual_stiffness>k</residual_stiffness>
+ <crack_resistance>gc</crack_resistance>
+ <crack_length_scale>ls</crack_length_scale>
+ <kinetic_coefficient>M</kinetic_coefficient>
+ <history_field>H</history_field>
+ </phasefield_parameters>
+ <solid_density>rho_sr</solid_density>
+ <process_variables>
+ <phasefield>phasefield</phasefield>
+ <displacement>displacement</displacement>
+ </process_variables>
+ <specific_body_force>0 0 0</specific_body_force>
+ </process>
+ </processes>
+ <time_loop>
+ <global_process_coupling>
+ <max_iter> 1000 </max_iter>
+ <convergence_criteria>
+ <!-- convergence criterion for the first process -->
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>INFINITY_N</norm_type>
+ <abstol>1.e-1</abstol>
+ <reltol>1.e-1</reltol>
+ </convergence_criterion>
+ <!-- convergence criterion for the second process -->
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>INFINITY_N</norm_type>
+ <abstol>1.e-4</abstol>
+ <reltol>1.e-10</reltol>
+ </convergence_criterion>
+ </convergence_criteria>
+ </global_process_coupling>
+ <processes>
+ <process ref="PhaseField">
+ <nonlinear_solver>basic_newton_u</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1.e-14</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable>displacement</variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>1</t_end>
+ <timesteps>
+ <pair>
+ <repeat>10</repeat>
+ <delta_t>1e-1</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ <process ref="PhaseField">
+ <nonlinear_solver>basic_newton_d</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1.e-14</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable>phasefield</variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>1</t_end>
+ <timesteps>
+ <pair>
+ <repeat>10</repeat>
+ <delta_t>1e-1</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>beam3_stag1pcsAT2</prefix>
+ <timesteps>
+ <pair>
+ <repeat>10</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ <pair>
+ <repeat>10</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ </timesteps>
+ </output>
+ </time_loop>
+ <parameters>
+ <!-- Mechanics -->
+ <parameter>
+ <name>E</name>
+ <type>Constant</type>
+ <value>2.1e5</value>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>Constant</type>
+ <value>0.3</value>
+ </parameter>
+ <parameter>
+ <name>k</name>
+ <type>Constant</type>
+ <value>1e-8</value>
+ </parameter>
+ <parameter>
+ <name>gc</name>
+ <type>Constant</type>
+ <value>2.7</value>
+ </parameter>
+ <parameter>
+ <name>ls</name>
+ <type>Constant</type>
+ <value>0.02</value>
+ </parameter>
+ <parameter>
+ <name>H</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0 0</values>
+ </parameter>
+ <parameter>
+ <name>phasefield_ic</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>phasefield_bc</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>dirichlet0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>Dirichlet_spatial</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>dirichlet1</name>
+ <type>CurveScaled</type>
+ <curve>Dirichlet_temporal</curve>
+ <parameter>Dirichlet_spatial</parameter>
+ </parameter>
+ <parameter>
+ <name>dirichlet2</name>
+ <type>CurveScaled</type>
+ <curve>Dirichlet_temporal_left</curve>
+ <parameter>Dirichlet_spatial</parameter>
+ </parameter>
+ <parameter>
+ <name>M</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>Dirichlet_temporal</name>
+ <coords>0 1</coords>
+ <values>0 0.013</values>
+ </curve>
+ <curve>
+ <name>Dirichlet_temporal_left</name>
+ <coords>0 1</coords>
+ <values>0 -0.013</values>
+ </curve>
+ </curves>
+ <process_variables>
+ <process_variable>
+ <name>phasefield</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>phasefield_ic</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>phasefield_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>phasefield_bc</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>displacement</name>
+ <components>3</components>
+ <order>1</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>back</geometry>
+ <type>Dirichlet</type>
+ <component>2</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>front</geometry>
+ <type>Dirichlet</type>
+ <component>2</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet2</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>beam3d_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet1</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton_d</name>
+ <type>Newton</type>
+ <max_iter>50</max_iter>
+ <linear_solver>linear_solver_d</linear_solver>
+ </nonlinear_solver>
+ <nonlinear_solver>
+ <name>basic_newton_u</name>
+ <type>Newton</type>
+ <max_iter>50</max_iter>
+ <linear_solver>linear_solver_u</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>linear_solver_d</name>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ </linear_solver>
+ <linear_solver>
+ <name>linear_solver_u</name>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/PhaseField/cracked_medium_ic.vtu b/Tests/Data/PhaseField/cracked_medium_ic.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..18ffcbda95817f483abb18ee1e691c407c898bc6
--- /dev/null
+++ b/Tests/Data/PhaseField/cracked_medium_ic.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f22a75eb26ee0fcd85f30edb7536f8fd3bb41a20fc3a6465674fbbeec760b4e4
+size 608625
diff --git a/Tests/Data/PhaseField/cube_1e0.prj b/Tests/Data/PhaseField/cube_1e0.prj
index 9c65f5f238e0286e0967c07f1481a6e13634ced8..015d7d44df0c9276317cdd3260440e04b38e555c 100644
--- a/Tests/Data/PhaseField/cube_1e0.prj
+++ b/Tests/Data/PhaseField/cube_1e0.prj
@@ -25,14 +25,8 @@
<phasefield>phasefield</phasefield>
</process_variables>
<secondary_variables>
- <secondary_variable type="static" internal_name="sigma_xx" output_name="sigma_xx"/>
- <secondary_variable type="static" internal_name="sigma_yy" output_name="sigma_yy"/>
- <secondary_variable type="static" internal_name="sigma_zz" output_name="sigma_zz"/>
- <secondary_variable type="static" internal_name="sigma_xy" output_name="sigma_xy"/>
- <secondary_variable type="static" internal_name="epsilon_xx" output_name="epsilon_xx"/>
- <secondary_variable type="static" internal_name="epsilon_yy" output_name="epsilon_yy"/>
- <secondary_variable type="static" internal_name="epsilon_zz" output_name="epsilon_zz"/>
- <secondary_variable type="static" internal_name="epsilon_xy" output_name="epsilon_xy"/>
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
<specific_body_force>0 0 0</specific_body_force>
</process>
@@ -53,14 +47,8 @@
<variables>
<variable>displacement</variable>
<variable>phasefield</variable>
- <variable>sigma_xx</variable>
- <variable>sigma_yy</variable>
- <variable>sigma_zz</variable>
- <variable>sigma_xy</variable>
- <variable>epsilon_xx</variable>
- <variable>epsilon_yy</variable>
- <variable>epsilon_zz</variable>
- <variable>epsilon_xy</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
</variables>
</output>
<time_stepping>
diff --git a/Tests/Data/PhaseField/expected_2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu b/Tests/Data/PhaseField/expected_2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..4c094b9b19beb756be46e249969967ed79a4955b
--- /dev/null
+++ b/Tests/Data/PhaseField/expected_2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7242128dbc6adc0dc37c92f9e962c53d5218b32d5ce51fb46aac600ca32fec82
+size 350563
diff --git a/Tests/Data/PhaseField/expected_beam3_stag1pcsAT2_pcs_1_ts_10_t_1.000000.vtu b/Tests/Data/PhaseField/expected_beam3_stag1pcsAT2_pcs_1_ts_10_t_1.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..8a9eff31f2ccb1c087d6c187cd4f4317a2ae12e5
--- /dev/null
+++ b/Tests/Data/PhaseField/expected_beam3_stag1pcsAT2_pcs_1_ts_10_t_1.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6ec0d3b7e9c65734abcc0f7bcfb6014c04ace62f4a0658d13f3704d5a8c31528
+size 42664
diff --git a/Tests/Data/PhaseField/square_line_h_400.prj b/Tests/Data/PhaseField/square_line_h_400.prj
index d6b2a4d726404d5f4f0923615dedb8aaa52a17d2..ce3f93e05d02278d891165085513f19dd654eb4c 100644
--- a/Tests/Data/PhaseField/square_line_h_400.prj
+++ b/Tests/Data/PhaseField/square_line_h_400.prj
@@ -25,14 +25,8 @@
<phasefield>phasefield</phasefield>
</process_variables>
<secondary_variables>
- <secondary_variable type="static" internal_name="sigma_xx" output_name="sigma_xx"/>
- <secondary_variable type="static" internal_name="sigma_yy" output_name="sigma_yy"/>
- <secondary_variable type="static" internal_name="sigma_zz" output_name="sigma_zz"/>
- <secondary_variable type="static" internal_name="sigma_xy" output_name="sigma_xy"/>
- <secondary_variable type="static" internal_name="epsilon_xx" output_name="epsilon_xx"/>
- <secondary_variable type="static" internal_name="epsilon_yy" output_name="epsilon_yy"/>
- <secondary_variable type="static" internal_name="epsilon_zz" output_name="epsilon_zz"/>
- <secondary_variable type="static" internal_name="epsilon_xy" output_name="epsilon_xy"/>
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
<specific_body_force>0 0</specific_body_force>
</process>
@@ -53,14 +47,8 @@
<variables>
<variable>displacement</variable>
<variable>phasefield</variable>
- <variable>sigma_xx</variable>
- <variable>sigma_yy</variable>
- <variable>sigma_zz</variable>
- <variable>sigma_xy</variable>
- <variable>epsilon_xx</variable>
- <variable>epsilon_yy</variable>
- <variable>epsilon_zz</variable>
- <variable>epsilon_xy</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
</variables>
</output>
<time_stepping>
diff --git a/Tests/Data/PhaseField/square_shear_h_400.prj b/Tests/Data/PhaseField/square_shear_h_400.prj
index 6eaabc34a16a38ca55662b41aeb7b680d2663b6f..5b569530e298f3fe6128fa8400032ff4f8aa7dbb 100644
--- a/Tests/Data/PhaseField/square_shear_h_400.prj
+++ b/Tests/Data/PhaseField/square_shear_h_400.prj
@@ -25,14 +25,8 @@
<phasefield>phasefield</phasefield>
</process_variables>
<secondary_variables>
- <secondary_variable type="static" internal_name="sigma_xx" output_name="sigma_xx"/>
- <secondary_variable type="static" internal_name="sigma_yy" output_name="sigma_yy"/>
- <secondary_variable type="static" internal_name="sigma_zz" output_name="sigma_zz"/>
- <secondary_variable type="static" internal_name="sigma_xy" output_name="sigma_xy"/>
- <secondary_variable type="static" internal_name="epsilon_xx" output_name="epsilon_xx"/>
- <secondary_variable type="static" internal_name="epsilon_yy" output_name="epsilon_yy"/>
- <secondary_variable type="static" internal_name="epsilon_zz" output_name="epsilon_zz"/>
- <secondary_variable type="static" internal_name="epsilon_xy" output_name="epsilon_xy"/>
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
<specific_body_force>0 0</specific_body_force>
</process>
@@ -53,14 +47,8 @@
<variables>
<variable>displacement</variable>
<variable>phasefield</variable>
- <variable>sigma_xx</variable>
- <variable>sigma_yy</variable>
- <variable>sigma_zz</variable>
- <variable>sigma_xy</variable>
- <variable>epsilon_xx</variable>
- <variable>epsilon_yy</variable>
- <variable>epsilon_zz</variable>
- <variable>epsilon_xy</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
</variables>
</output>
<time_stepping>