From 3bb4a25b9227db5e0ba6840936b8b7393ae684a0 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Fri, 14 Feb 2020 16:17:21 +0100
Subject: [PATCH] [PL] RM; Porosity element average output.
---
ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h | 6 ++++++
ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp | 4 ++++
ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h | 1 +
3 files changed, 11 insertions(+)
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index a283e400a3c..0f0126813e3 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -910,6 +910,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
_integration_method.getNumberOfPoints();
double saturation_avg = 0;
+ double porosity_avg = 0;
using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
KV sigma_avg = KV::Zero();
@@ -936,12 +937,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
saturation_avg += S_L;
+ porosity_avg +=
+ _ip_data[ip].porosity; // Note, this is not updated, because needs
+ // xdot and dt to be passed.
sigma_avg += _ip_data[ip].sigma_eff;
}
saturation_avg /= n_integration_points;
+ porosity_avg /= n_integration_points;
sigma_avg /= n_integration_points;
(*_process_data.element_saturation)[_element.getID()] = saturation_avg;
+ (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
KV::RowsAtCompileTime]) =
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index bfbf25f5d34..c05497d0809 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -206,6 +206,10 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
MeshLib::MeshItemType::Cell, 1);
+ _process_data.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
+ MeshLib::MeshItemType::Cell, 1);
+
_process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "stress_avg",
MeshLib::MeshItemType::Cell,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index 97dfcfb101e..93511f2e8fc 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -58,6 +58,7 @@ struct RichardsMechanicsProcessData
bool const apply_mass_lumping;
MeshLib::PropertyVector<double>* element_saturation = nullptr;
+ MeshLib::PropertyVector<double>* element_porosity = nullptr;
MeshLib::PropertyVector<double>* element_stresses = nullptr;
MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
--
GitLab