/**
 * \author Norihiro Watanabe
 * \date   2013-08-13
 *
 * \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 <boost/math/constants/constants.hpp>
#include <cassert>

#include "MeshLib/Elements/Element.h"
#include "MeshLib/Node.h"
#include "NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"

namespace NumLib
{
/**
 * \brief Template class for isoparametric elements
 *
 * \tparam ShapeFunctionType_   The shape function type.
 * \tparam ShapeMatrixTypes_    An aggregate of shape matrix types.
 */
template <class ShapeFunctionType_, class ShapeMatrixTypes_>
class TemplateIsoparametric
{
public:
    using ShapeFunctionType = ShapeFunctionType_;

    /// Coordinate mapping matrices type.
    using ShapeMatrices = typename ShapeMatrixTypes_::ShapeMatrices;

    /// Natural coordinates mapping tools specialization for specific
    /// MeshElement and ShapeFunction types.
    using NaturalCoordsMappingType =
        NaturalCoordinatesMapping<ShapeFunctionType, ShapeMatrices>;

    /**
     * Constructor without specifying a mesh element. setMeshElement() must be
     * called afterwards.
     */
    TemplateIsoparametric() : _ele(nullptr) {}

    /**
     * Construct this object for the given mesh element.
     *
     * @param e                      Mesh element object
     */
    explicit TemplateIsoparametric(const MeshLib::Element& e) : _ele(&e) {}

    /// return current mesh element
    const MeshLib::Element* getMeshElement() const { return _ele; }

    /// Sets the mesh element
    void setMeshElement(const MeshLib::Element& e) { this->_ele = &e; }
    /**
     * compute shape functions
     *
     * @param natural_pt            position in natural coordinates
     * @param shape                 evaluated shape function matrices
     * @param global_dim            global dimension
     * @param is_axially_symmetric  if true, the integral measure for cylinder
     *                              coordinates is used, and 1 otherwise.
     */
    void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
                               const unsigned global_dim,
                               bool is_axially_symmetric) const
    {
        NaturalCoordsMappingType::computeShapeMatrices(*_ele, natural_pt, shape,
                                                       global_dim);
        computeIntegralMeasure(is_axially_symmetric, shape);
    }

    /**
     * compute shape functions
     *
     * @tparam T_SHAPE_MATRIX_TYPE  shape matrix types to be calculated
     * @param natural_pt            position in natural coordinates
     * @param shape                 evaluated shape function matrices
     * @param global_dim            global dimension
     * @param is_axially_symmetric  if true, the integral measure for cylinder
     *                              coordinates is used, and 1 otherwise.
     */
    template <ShapeMatrixType T_SHAPE_MATRIX_TYPE>
    void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
                               const unsigned global_dim,
                               bool is_axially_symmetric) const
    {
        NaturalCoordsMappingType::template computeShapeMatrices<
            T_SHAPE_MATRIX_TYPE>(*_ele, natural_pt, shape, global_dim);
        computeIntegralMeasure(is_axially_symmetric, shape);
    }

    /// Interpolates the x coordinate of the element with the given shape
    /// matrix.
    double interpolateZerothCoordinate(
        typename ShapeMatrices::ShapeType const& N) const
    {
        auto* nodes = _ele->getNodes();
        typename ShapeMatrices::ShapeType rs(N.size());
        for (int i = 0; i < rs.size(); ++i)
        {
            rs[i] = (*nodes[i])[0];
        }
        auto const r = N.dot(rs);
        return r;
    }

    /// Interpolates the coordinates of the element with the given shape matrix.
    std::array<double, 3> interpolateCoordinates(
        typename ShapeMatrices::ShapeType const& N) const
    {
        auto* nodes = _ele->getNodes();
        typename ShapeMatrices::ShapeType rs(N.size());
        std::array<double, 3> interpolated_coords;
        for (int d = 0; d < 3; ++d)
        {
            for (int i = 0; i < rs.size(); ++i)
            {
                rs[i] = (*nodes[i])[d];
            }
            interpolated_coords[d] = N.dot(rs);
        }
        return interpolated_coords;
    }

private:
    void computeIntegralMeasure(bool is_axially_symmetric,
                                ShapeMatrices& shape) const
    {
        if (!is_axially_symmetric)
        {
            shape.integralMeasure = 1.0;
            return;
        }

        // Note: If an integration point is located at the rotation axis, r will
        // be zero which might lead to problems with the assembled equation
        // system.
        // E.g., for triangle elements, if an integration point is
        // located at edge of the unit triangle, it is possible that
        // r becomes zero.
        auto const r = interpolateZerothCoordinate(shape.N);
        shape.integralMeasure = boost::math::constants::two_pi<double>() * r;
    }

    const MeshLib::Element* _ele;
};

/// Creates a TemplateIsoparametric element for the given shape functions and
/// the underlying mesh element.
template <typename ShapeFunction, typename ShapeMatricesType>
NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>
createIsoparametricFiniteElement(MeshLib::Element const& e)
{
    using FemType =
        NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;

    return FemType{e};
}
}  // namespace NumLib