From 8d2b00fd2325c46f98942c432405df9b1a211e48 Mon Sep 17 00:00:00 2001
From: "Dmitry Yu. Naumov" <github@naumov.de>
Date: Tue, 4 Sep 2018 16:09:57 +0200
Subject: [PATCH] [PL] RM; Interpolated pressure output.
This is copy from LIE/HM and needs to be generalized.
---
.../RichardsMechanicsFEM-impl.h | 42 +++++++++++++++++++
.../RichardsMechanicsProcess.cpp | 5 +++
.../RichardsMechanicsProcessData.h | 1 +
3 files changed, 48 insertions(+)
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index c7d2fb5f04e..6d09a5688bd 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -15,6 +15,7 @@
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/KelvinVector.h"
+#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
#include "NumLib/Function/Interpolation.h"
#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -832,6 +833,47 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
saturation_avg /= n_integration_points;
(*_process_data.element_saturation)[_element.getID()] = saturation_avg;
+
+ // For each higher order node evaluate the shape matrices for the lower
+ // order element (the base nodes)
+ // TODO (naumov) Extract this method to be useful for other processes.
+ auto interpolate_p = [&]() {
+ using FemType =
+ NumLib::TemplateIsoparametric<ShapeFunctionPressure,
+ ShapeMatricesTypePressure>;
+
+ FemType fe(
+ *static_cast<const typename ShapeFunctionPressure::MeshElement*>(
+ &_element));
+ int const number_base_nodes = _element.getNumberOfBaseNodes();
+ int const number_all_nodes = _element.getNumberOfNodes();
+
+ for (int n = 0; n < number_base_nodes; ++n)
+ {
+ std::size_t const global_index = _element.getNodeIndex(n);
+ (*_process_data.pressure_interpolated)[global_index] = p_L[n];
+ }
+
+ for (int n = number_base_nodes; n < number_all_nodes; ++n)
+ {
+ // Evaluated at higher order nodes' coordinates.
+ typename ShapeMatricesTypePressure::ShapeMatrices shape_matrices_p{
+ ShapeFunctionPressure::DIM, DisplacementDim,
+ ShapeFunctionPressure::NPOINTS};
+
+ fe.computeShapeFunctions(
+ NumLib::NaturalCoordinates<typename ShapeFunctionDisplacement::
+ MeshElement>::coordinates[n]
+ .data(),
+ shape_matrices_p, DisplacementDim, _is_axially_symmetric);
+
+ auto const& N_p = shape_matrices_p.N;
+
+ std::size_t const global_index = _element.getNodeIndex(n);
+ (*_process_data.pressure_interpolated)[global_index] = N_p * p_L;
+ }
+ };
+ interpolate_p();
}
} // namespace RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index d1ef1a3cd66..2765c51b04d 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -195,6 +195,11 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
_process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
MeshLib::MeshItemType::Cell, 1);
+
+ _process_data.pressure_interpolated =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
+ MeshLib::MeshItemType::Node, 1);
}
template <int DisplacementDim>
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index f737f7ce28a..c788c645b81 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -104,6 +104,7 @@ struct RichardsMechanicsProcessData
double t = 0.0;
MeshLib::PropertyVector<double>* element_saturation = nullptr;
+ MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
--
GitLab