From f982090688df016cbb1d666ce47b7aac31ce7ec8 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Mon, 17 Sep 2018 09:49:11 +0200
Subject: [PATCH] [PL] Move interpolation algorithm to NumLib.
---
NumLib/Function/Interpolation.h | 51 +++++++++++++++++
.../RichardsMechanicsFEM-impl.h | 55 ++-----------------
2 files changed, 55 insertions(+), 51 deletions(-)
diff --git a/NumLib/Function/Interpolation.h b/NumLib/Function/Interpolation.h
index 8999d16d0b4..876ac7d8ee0 100644
--- a/NumLib/Function/Interpolation.h
+++ b/NumLib/Function/Interpolation.h
@@ -12,6 +12,9 @@
#include<array>
#include<cassert>
+#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+
namespace NumLib
{
@@ -85,4 +88,52 @@ void shapeFunctionInterpolate(
interpolated_values...);
}
+/// Interpolates \c node_values given in lower order element nodes (e.g. the
+/// base nodes) to higher order element's nodes (e.g. quadratic nodes) and
+/// writes the result into the global property vector.
+///
+/// The base nodes' values are copied. For each higher order node the shape
+/// matrices are evaluated for the lower order element (the base nodes), and
+/// used for the the scalar quantity interpolation.
+template <typename LowerOrderShapeFunction, typename HigherOrderMeshElementType,
+ int GlobalDim, typename EigenMatrixType>
+void interpolateToHigherOrderNodes(
+ MeshLib::Element const& element, bool const is_axially_symmetric,
+ Eigen::MatrixBase<EigenMatrixType> const& node_values,
+ MeshLib::PropertyVector<double>& interpolated_values_global_vector)
+{
+ using SF = LowerOrderShapeFunction;
+ using ShapeMatricesType = ShapeMatrixPolicyType<SF, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+ using FemType = TemplateIsoparametric<SF, ShapeMatricesType>;
+
+ FemType fe(*static_cast<const typename SF::MeshElement*>(&element));
+ int const number_base_nodes = element.getNumberOfBaseNodes();
+ int const number_all_nodes = element.getNumberOfNodes();
+
+ // Copy the values for linear nodes.
+ for (int n = 0; n < number_base_nodes; ++n)
+ {
+ std::size_t const global_index = element.getNodeIndex(n);
+ interpolated_values_global_vector[global_index] = node_values[n];
+ }
+
+ // Shape matrices storage reused in the interpolation loop.
+ ShapeMatrices shape_matrices{SF::DIM, GlobalDim, SF::NPOINTS};
+
+ // Interpolate values for higher order nodes.
+ for (int n = number_base_nodes; n < number_all_nodes; ++n)
+ {
+ // Evaluated at higher order nodes' coordinates.
+ fe.template computeShapeFunctions<ShapeMatrixType::N>(
+ NaturalCoordinates<HigherOrderMeshElementType>::coordinates[n]
+ .data(),
+ shape_matrices, GlobalDim, is_axially_symmetric);
+
+ std::size_t const global_index = element.getNodeIndex(n);
+ interpolated_values_global_vector[global_index] =
+ shape_matrices.N * node_values;
+ }
+}
+
} // namespace NumLib
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 0442978c6b5..f386dc0ddd5 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -15,7 +15,6 @@
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/KelvinVector.h"
-#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
#include "NumLib/Function/Interpolation.h"
#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -23,53 +22,6 @@ namespace ProcessLib
{
namespace RichardsMechanics
{
-/// For each higher order node evaluate the shape matrices for the lower order
-/// element (the base nodes)
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim, typename EigenMatrixType>
-void interpolate_p(MeshLib::Element const& element,
- bool const is_axially_symmetric,
- Eigen::MatrixBase<EigenMatrixType> const& p_L,
- MeshLib::PropertyVector<double>& pressure_interpolated)
-{
- using ShapeMatricesTypePressure =
- ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
- 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();
-
- // Copy the values for linear nodes.
- for (int n = 0; n < number_base_nodes; ++n)
- {
- std::size_t const global_index = element.getNodeIndex(n);
- pressure_interpolated[global_index] = p_L[n];
- }
-
- // Interpolate values for higher order nodes.
- 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.template computeShapeFunctions<NumLib::ShapeMatrixType::N>(
- 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);
- pressure_interpolated[global_index] = N_p * p_L;
- }
-}
-
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -881,9 +833,10 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
(*_process_data.element_saturation)[_element.getID()] = saturation_avg;
- interpolate_p<ShapeFunctionDisplacement, ShapeFunctionPressure,
- DisplacementDim>(_element, _is_axially_symmetric, p_L,
- *_process_data.pressure_interpolated);
+ NumLib::interpolateToHigherOrderNodes<
+ ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
+ DisplacementDim>(_element, _is_axially_symmetric, p_L,
+ *_process_data.pressure_interpolated);
}
} // namespace RichardsMechanics
--
GitLab