OGS
Interpolation.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <array>
14#include <cassert>
15
20
21namespace NumLib
22{
23namespace detail
24{
26template <unsigned DOFOffset, typename NodalValues, typename ShapeMatrix>
27void shapeFunctionInterpolate(const NodalValues& /*nodal_values*/,
28 const ShapeMatrix& /*shape_matrix_N*/)
29{
30}
31
33template <unsigned DOFOffset, typename NodalValues, typename ShapeMatrix,
34 typename... ScalarTypes>
35void shapeFunctionInterpolate(const NodalValues& nodal_values,
36 const ShapeMatrix& shape_matrix_N,
37 double& interpolated_value,
38 ScalarTypes&... interpolated_values)
39{
40 auto const num_nodes = shape_matrix_N.size();
41 double iv = 0.0;
42
43 for (auto n = decltype(num_nodes){0}; n < num_nodes; ++n)
44 {
45 iv += nodal_values[DOFOffset * num_nodes + n] * shape_matrix_N[n];
46 }
47
48 interpolated_value = iv;
49
50 shapeFunctionInterpolate<DOFOffset + 1>(nodal_values, shape_matrix_N,
51 interpolated_values...);
52}
53
54} // namespace detail
55
79template <typename NodalValues, typename ShapeMatrix, typename... ScalarTypes>
80void shapeFunctionInterpolate(const NodalValues& nodal_values,
81 const ShapeMatrix& shape_matrix_N,
82 double& interpolated_value,
83 ScalarTypes&... interpolated_values)
84{
85 auto const num_nodal_dof = sizeof...(interpolated_values) + 1;
86 auto const num_nodes = shape_matrix_N.size();
87
88 assert(num_nodes * num_nodal_dof ==
89 static_cast<std::size_t>(nodal_values.size()));
90 (void)num_nodal_dof;
91 (void)num_nodes; // no warnings when not in debug build
92
93 detail::shapeFunctionInterpolate<0>(nodal_values, shape_matrix_N,
94 interpolated_value,
95 interpolated_values...);
96}
97
105template <typename LowerOrderShapeFunction, typename HigherOrderMeshElementType,
106 int GlobalDim, typename EigenMatrixType>
108 MeshLib::Element const& element, bool const is_axially_symmetric,
109 Eigen::MatrixBase<EigenMatrixType> const& node_values,
110 MeshLib::PropertyVector<double>& interpolated_values_global_vector)
111{
112 assert(dynamic_cast<HigherOrderMeshElementType const*>(&element));
113 assert(node_values.cols() == 1); // Scalar quantity only.
114
115 using SF = LowerOrderShapeFunction;
116 using ShapeMatricesType = ShapeMatrixPolicyType<SF, GlobalDim>;
117
118 int const number_base_nodes = element.getNumberOfBaseNodes();
119 int const number_all_nodes = element.getNumberOfNodes();
120
121 // Copy the values for linear nodes.
122 for (int n = 0; n < number_base_nodes; ++n)
123 {
124 std::size_t const global_index = getNodeIndex(element, n);
125 interpolated_values_global_vector[global_index] = node_values[n];
126 }
127
128 //
129 // Interpolate values for higher order nodes.
130 //
131
132 int const number_higher_order_nodes = number_all_nodes - number_base_nodes;
133 std::vector<MathLib::Point3d> higher_order_nodes;
134 higher_order_nodes.reserve(number_higher_order_nodes);
135 for (int n = 0; n < number_higher_order_nodes; ++n)
136 {
137 higher_order_nodes.emplace_back(
139 [number_base_nodes + n]);
140 }
141
142 // Shape matrices evaluated at higher order nodes' coordinates.
143 auto const shape_matrices =
144 computeShapeMatrices<SF, ShapeMatricesType, GlobalDim,
145 ShapeMatrixType::N>(element, is_axially_symmetric,
146 higher_order_nodes);
147
148 for (int n = 0; n < number_higher_order_nodes; ++n)
149 {
150 std::size_t const global_index =
151 getNodeIndex(element, number_base_nodes + n);
152 interpolated_values_global_vector[global_index] =
153 shape_matrices[n].N * node_values;
154 }
155}
156
157} // namespace NumLib
virtual unsigned getNumberOfNodes() const =0
virtual unsigned getNumberOfBaseNodes() const =0
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)