/**
 * \file
 * \copyright
 * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

#pragma once

#include <array>
#include <cassert>

#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/InitShapeMatrices.h"

namespace NumLib
{

namespace detail
{

//! \see ::NumLib::shapeFunctionInterpolate()
template<unsigned DOFOffset, typename NodalValues, typename ShapeMatrix>
void shapeFunctionInterpolate(
        const NodalValues &/*nodal_values*/,
        const ShapeMatrix &/*shape_matrix_N*/)
{}

//! \see ::NumLib::shapeFunctionInterpolate()
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

/*!
 * Interpolates variables given at element nodes according to the given shape matrix.
 *
 * This function simply does the usual finite-element interpolation, i.e. multiplication
 * of nodal values with the shape function.
 *
 * \param nodal_values vector of nodal values, ordered by component
 * \param shape_matrix_N shape matrix of the point to which will be interpolated
 * \param interpolated_value interpolated value of the first d.o.f. (output parameter)
 * \param interpolated_values interpolated value of further d.o.f. (output parameter)
 *
 * \tparam NodalValues  type of the container where nodal values are stored
 * \tparam ShapeMatrix  type of the shape matrix \f$N\f$.
 * \tparam ScalarValues all of the types in this pack must currently be \c double.
 *
 * \note
 * \c nodal_values have to be ordered by component and it is assumed that all passed d.o.f. are
 * single-component and are interpolated using the same shape function.
 */
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...);
}

/// Interpolates scalar \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)
{
    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 = element.getNodeIndex(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 =
            element.getNodeIndex(number_base_nodes + n);
        interpolated_values_global_vector[global_index] =
            shape_matrices[n].N * node_values;
    }
}

} // namespace NumLib