Interpolation.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <array>
#include <cassert>
 
#include "CoordinatesMapping/NaturalNodeCoordinates.h"
#include "FiniteElement/TemplateIsoparametric.h"
#include "InitShapeMatrices.h"
#include "ShapeMatrixPolicy.h"
 
namespace NumLib
{
namespace detail
{
template <unsigned DOFOffset, typename NodalValues, typename ShapeMatrix>
void shapeFunctionInterpolate(const NodalValues& /*nodal_values*/,
                              const ShapeMatrix& /*shape_matrix_N*/)
{
}

template <unsigned DOFOffset, typename NodalValues, typename ShapeMatrix,
          typename... ScalarTypes>
void shapeFunctionInterpolate(const NodalValues& nodal_values,
                              const ShapeMatrix& shape_matrix_N,
                              double& interpolated_value,
                              ScalarTypes&... interpolated_values)
{
    auto const num_nodes = shape_matrix_N.size();
    double iv = 0.0;
 
    for (auto n = decltype(num_nodes){0}; n < num_nodes; ++n)
    {
        iv += nodal_values[DOFOffset * num_nodes + n] * shape_matrix_N[n];
    }
 
    interpolated_value = iv;
 
    shapeFunctionInterpolate<DOFOffset + 1>(nodal_values, shape_matrix_N,
                                            interpolated_values...);
}
 
}  // namespace detail

template <typename NodalValues, typename ShapeMatrix, typename... ScalarTypes>
void shapeFunctionInterpolate(const NodalValues& nodal_values,
                              const ShapeMatrix& shape_matrix_N,
                              double& interpolated_value,
                              ScalarTypes&... interpolated_values)
{
    auto const num_nodal_dof = sizeof...(interpolated_values) + 1;
    auto const num_nodes = shape_matrix_N.size();
 
    assert(num_nodes * num_nodal_dof ==
           static_cast<std::size_t>(nodal_values.size()));
    (void)num_nodal_dof;
    (void)num_nodes;  // no warnings when not in debug build
 
    detail::shapeFunctionInterpolate<0>(nodal_values, shape_matrix_N,
                                        interpolated_value,
                                        interpolated_values...);
}

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)
{
    assert(dynamic_cast<HigherOrderMeshElementType const*>(&element));
    assert(node_values.cols() == 1);  // Scalar quantity only.
 
    using SF = LowerOrderShapeFunction;
    using ShapeMatricesType = ShapeMatrixPolicyType<SF, GlobalDim>;
 
    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 = getNodeIndex(element, n);
        interpolated_values_global_vector[global_index] = node_values[n];
    }
 
    //
    // Interpolate values for higher order nodes.
    //
 
    int const number_higher_order_nodes = number_all_nodes - number_base_nodes;
    std::vector<MathLib::Point3d> higher_order_nodes;
    higher_order_nodes.reserve(number_higher_order_nodes);
    for (int n = 0; n < number_higher_order_nodes; ++n)
    {
        higher_order_nodes.emplace_back(
            NaturalCoordinates<HigherOrderMeshElementType>::coordinates
                [number_base_nodes + n]);
    }
 
    // Shape matrices evaluated at higher order nodes' coordinates.
    auto const shape_matrices =
        computeShapeMatrices<SF, ShapeMatricesType, GlobalDim,
                             ShapeMatrixType::N>(element, is_axially_symmetric,
                                                 higher_order_nodes);
 
    for (int n = 0; n < number_higher_order_nodes; ++n)
    {
        std::size_t const global_index =
            getNodeIndex(element, number_base_nodes + n);
        interpolated_values_global_vector[global_index] =
            shape_matrices[n].N * node_values;
    }
}
 
}  // namespace NumLib