From 3410895e4005a09528ff08f137aaa31dda9853d1 Mon Sep 17 00:00:00 2001
From: KeitaYoshioka <keita.yoshioka@ufz.de>
Date: Wed, 7 Feb 2018 13:04:36 +0100
Subject: [PATCH] [PL] PF: Add energy computation at postTimestep.
---
.../LinearElasticIsotropicPhaseField.h | 6 +-
MaterialLib/SolidModels/PhaseFieldExtension.h | 9 ++-
.../PhaseField/LocalAssemblerInterface.h | 10 +++
ProcessLib/PhaseField/PhaseFieldFEM-impl.h | 79 +++++++++++++++++++
ProcessLib/PhaseField/PhaseFieldFEM.h | 15 +++-
.../PhaseField/PhaseFieldProcess-impl.h | 27 ++++++-
ProcessLib/PhaseField/PhaseFieldProcessData.h | 3 +
7 files changed, 139 insertions(+), 10 deletions(-)
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
index 7bce42d01c2..7314b1130fd 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 4ef4dbfea62..425da16cf2b 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/PhaseField/LocalAssemblerInterface.h b/ProcessLib/PhaseField/LocalAssemblerInterface.h
index c7f88e48888..e8d6a46ae1b 100644
--- a/ProcessLib/PhaseField/LocalAssemblerInterface.h
+++ b/ProcessLib/PhaseField/LocalAssemblerInterface.h
@@ -102,6 +102,16 @@ struct PhaseFieldLocalAssemblerInterface
GlobalVector const& x, double const t, double& crack_volume,
bool const use_monolithic_scheme,
CoupledSolutionsForStaggeredScheme const* const cpl_xs) = 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
index b41a17021bc..e44074d184c 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -502,5 +502,84 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
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 0fe7ed00c99..f524a65d3a1 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;
};
@@ -161,6 +161,7 @@ public:
_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;
@@ -217,6 +218,16 @@ public:
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
{
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
index 684c1f08a46..3c7cabeb4a4 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
@@ -360,11 +360,30 @@ template <int DisplacementDim>
void PhaseFieldProcess<DisplacementDim>::postTimestepConcreteProcess(
GlobalVector const& x, int const process_id)
{
- DBUG("PostTimestep PhaseFieldProcess %d.", process_id);
+ if (isPhaseFieldProcess(process_id))
+ {
+ DBUG("PostTimestep PhaseFieldProcess.");
- GlobalExecutor::executeMemberOnDereferenced(
- &LocalAssemblerInterface::postTimestep, _local_assemblers,
- getDOFTable(process_id), 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>
diff --git a/ProcessLib/PhaseField/PhaseFieldProcessData.h b/ProcessLib/PhaseField/PhaseFieldProcessData.h
index 26fcde7b5b8..75fdace0717 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcessData.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcessData.h
@@ -98,6 +98,9 @@ struct PhaseFieldProcessData
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;
};
--
GitLab