/**
 * \file
 * \copyright
 * Copyright (c) 2012-2023, 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 <vector>

#include "FiniteElement/TemplateIsoparametric.h"
#include "MathLib/WeightedPoint.h"
#include "MeshLib/Elements/Element.h"

namespace NumLib
{
template <typename ShapeFunction, typename ShapeMatricesType, int GlobalDim,
          ShapeMatrixType SelectedShapeMatrixType = ShapeMatrixType::ALL,
          typename PointContainer>
std::vector<typename ShapeMatricesType::ShapeMatrices,
            Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
computeShapeMatrices(MeshLib::Element const& e, bool const is_axially_symmetric,
                     PointContainer const& points)
{
    std::vector<
        typename ShapeMatricesType::ShapeMatrices,
        Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
        shape_matrices;

    auto const fe =
        createIsoparametricFiniteElement<ShapeFunction, ShapeMatricesType>(e);

    shape_matrices.reserve(points.size());
    for (auto const& p : points)
    {
        shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
                                    ShapeFunction::NPOINTS);
        fe.template computeShapeFunctions<SelectedShapeMatrixType>(
            p.data(), shape_matrices.back(), GlobalDim, is_axially_symmetric);
    }

    return shape_matrices;
}

template <typename ShapeFunction, typename ShapeMatricesType, int GlobalDim,
          ShapeMatrixType SelectedShapeMatrixType = ShapeMatrixType::ALL,
          typename IntegrationMethod>
std::vector<typename ShapeMatricesType::ShapeMatrices,
            Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
initShapeMatrices(MeshLib::Element const& e, bool const is_axially_symmetric,
                  IntegrationMethod const& integration_method)
{
    int const n_integration_points = integration_method.getNumberOfPoints();

    std::vector<MathLib::WeightedPoint> points;
    points.reserve(n_integration_points);
    for (int ip = 0; ip < n_integration_points; ++ip)
    {
        points.push_back(integration_method.getWeightedPoint(ip));
    }

    return computeShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim,
                                SelectedShapeMatrixType>(
        e, is_axially_symmetric, points);
}

template <typename ShapeFunction, typename ShapeMatricesType>
double interpolateXCoordinate(
    MeshLib::Element const& e,
    typename ShapeMatricesType::ShapeMatrices::ShapeType const& N)
{
    auto const fe =
        createIsoparametricFiniteElement<ShapeFunction, ShapeMatricesType>(e);

    return fe.interpolateZerothCoordinate(N);
}

template <typename ShapeFunction, typename ShapeMatricesType>
std::array<double, 3> interpolateCoordinates(
    MeshLib::Element const& e,
    typename ShapeMatricesType::ShapeMatrices::ShapeType const& N)
{
    auto const fe =
        createIsoparametricFiniteElement<ShapeFunction, ShapeMatricesType>(e);

    return fe.interpolateCoordinates(N);
}

}  // namespace NumLib