#ifndef MEDUSA_BITS_DOMAINS_NURBSPATCH_FWD_HPP_
#define MEDUSA_BITS_DOMAINS_NURBSPATCH_FWD_HPP_

/**
 * @file
 * Declaration of the NURBS patch class.
 *
 * @example test/domains/NURBSPatch_test.cpp
 */

#include <medusa/Config.hpp>
#include <medusa/bits/types/Vec.hpp>
#include <medusa/bits/types/Range.hpp>
#include <medusa/bits/domains/BoxShape.hpp>
#include <memory>

namespace mm
{

    /**
     * Implementation details of NURBSPatch, contains structures for partial class specialization.
     */
    namespace nurbs_patch_internal
    {
        /**
         * Internal structure of NURBSPatch that helps with partial class specialization.
         * All of its functions have the same signature (apart from accepting a reference to NURBSPatch)
         * as the corresponding functions of NURBSPatch. NURBSPatch holds an instance of this class
         * and calls its functions when specialized functions are requested.
         *
         * Note: partial class specialization can also be done by having a base class with functions,
         * that are not specialized (NURBSPatchBase) and a fully specialized class (NURBSPatch) that
         * inherits from it. This approach was chosen to keep Medusa's class hierarchy cleaner.
         */
        template <typename vec_t, typename param_vec_t>
        struct NURBSPatchHelper
        {
        };
        /**
         * Struct for declaring a nested Range.
         * @tparam dim Number of dimensions (nestings).
         * @tparam T Type.
         */
        template <size_t dim, typename T>
        struct NestedRange
        {
            /// Type of the multidimensional Range.
            typedef Range<typename NestedRange<dim - 1, T>::type> type;
        };

        /// @cond
        template <typename T>
        struct NestedRange<0, T>
        {
            typedef T type;
        };
        /// @endcond
    } // namespace nurbs_patch_internal

    /**
     * Class representing a single NURBS patch in an arbitrary dimensional space,
     * defined on an arbitrary dimensional domain and generated by
     * a tensor product of NURBS curves.
     *
     * @warning Currently only supports NURBS surfaces (2D domain) and NURBS
     * curves (1D domain).
     *
     * Usage example:
     * @snippet domains/NURBSPatch_test.cpp NURBSPatch usage example
     *
     * @tparam param_vec_t Vector type of the parametric space.
     * @tparam vec_t Vector type of ambient space.
     *
     * @ingroup domains
     */
    template <typename vec_t, typename param_vec_t>
    class NURBSPatch
    {
        friend struct nurbs_patch_internal::NURBSPatchHelper<vec_t, param_vec_t>;

    public:
        typedef typename vec_t::scalar_t scalar_t; ///< Scalar type.
        /// Store dimensionality.
        enum
        { /** Dimensionality of space after projection. */ dim = vec_t::dim,
          /** Dimensionality of space before projection. */ proj_dim = dim + 1,
          /** Dimensionality of parametric space. */ param_dim = param_vec_t::dim };
        typedef Vec<scalar_t, proj_dim> proj_vec_t; ///< Vector type before projection.

        template <typename T>
        /// Range of param_dim dimensions.
        using NdRange = typename nurbs_patch_internal::NestedRange<param_dim, T>::type;

    private:
        std::array<int, param_dim> p;                 ///< Underlying B-spline orders.
        NdRange<proj_vec_t> weighted_control_points;  ///< Multidimensional Range of control points.
        std::array<Range<scalar_t>, param_dim> knots; ///< Array of Ranges of knots (knot vector).

        /// NURBS surfaces needed for derivative calculation.
        std::unique_ptr<std::array<NURBSPatch<vec_t, param_vec_t>, param_dim>> der_structure;

    public:
        /**
         * Construct NURBS patch with weighted control points.
         * @param wcp Multidimensional Range of weighted control points,
         * e. g. `Vec<double, 4>(w * x, w * y, w * z, w)`.
         * @param ks Array of Ranges of knots (pair of knot vectors).
         * @param in_p Array of underlying B-spline orders (starting from order 0).
         */
        NURBSPatch(const NdRange<proj_vec_t> &wcp,
                   const std::array<Range<scalar_t>, param_dim> &ks,
                   const std::array<int, param_dim> &in_p);

        /**
         * Construct NURBS patch with control points and weights.
         * @param cp Multidimensional Range of control points.
         * @param w Multidimensional Range of weights.
         * @param ks Array of Ranges of knots (pair of knot vectors).
         * @param in_p Array of underlying B-spline orders (starting from order 0).
         */
        NURBSPatch(const NdRange<vec_t> &cp, const NdRange<double> &w,
                   const std::array<Range<scalar_t>, param_dim> &ks,
                   const std::array<int, param_dim> &in_p);

        /// Copy constructor.
        NURBSPatch(const NURBSPatch<vec_t, param_vec_t> &patch) noexcept;

        /// Copy assignment.
        NURBSPatch<vec_t, param_vec_t> &operator=(const NURBSPatch<vec_t, param_vec_t> &) noexcept;

        /// Move constructor.
        NURBSPatch(NURBSPatch<vec_t, param_vec_t> &&) noexcept = default;
        /// Move assignment.
        NURBSPatch<vec_t, param_vec_t> &operator=(NURBSPatch<vec_t, param_vec_t> &&) noexcept = default;

        /// Calculate data needed for derivative evaluation.
        void computeDerivativeStructure();

        /**
         * Evaluate NURBS patch in one point along the second dimension.
         * @param t Point of evaluation.
         * @param w Variable in which the weight is stored.
         * @return Vector to the point on the NURBS surface corresponding to parameter `t`.
         */
        vec_t evaluate(const param_vec_t &t, scalar_t *w) const;

        /**
         * Overload if weight is not required.
         * @param t Point of evaluation.
         * @return Vector to the point on the NURBS surface corresponding to parameter `t`.
         */
        vec_t evaluate(const param_vec_t &t) const;

        /**
         * Evaluate NURBS patch in one point.
         * @param t Point of evaluation.
         * @return Vector to the weighted point on the NURBS patch corresponding to parameter `t`.
         */
        proj_vec_t evaluateWeighted(const param_vec_t &t) const;

        /**
         * Evaluate NURBS jacobian matrix in one point.
         * @param t Point of evaluation.
         * @param pt Point on the NURBS patch corresponding to `t`.
         * @param w Weight of `pt`.
         * @return Jacobian matrix of the NURBS patch at point `t`.
         * @warning `computeDerivativeStructure()` must be called on the object before calling
         * this function.
         */
        Eigen::Matrix<scalar_t, dim, param_dim> jacobian(const param_vec_t &t, const vec_t &pt,
                                                         const scalar_t &w) const;

        /**
         * Overload that calculates the point automatically.
         * @warning Less efficient.
         * @param t Point of evaluation.
         * @return Jacobian matrix of the NURBS patch at point `t`.
         * @warning `computeDerivativeStructure()` must be called on the object before calling
         * this function.
         */
        Eigen::Matrix<scalar_t, dim, param_dim> jacobian(const param_vec_t &t) const;

        /**
         * Evaluate NURBS and its jacobian matrix in one parameter.
         * @param t Parameter of evaluation.
         * @return Pair of point on the NURBS and jacobian matrix of the NURBS at parameter `t`.
         */
        std::pair<vec_t, Eigen::Matrix<scalar_t, dim, param_dim>> evaluatePointAndJacobian(
            const param_vec_t &t) const;

        /**
         * Get domain of the NURBS patch.
         * @return Domain of the NURBS in the form of a BoxShape.
         */
        BoxShape<param_vec_t> getDomain() const;

        /**
         * Get all boundaries of the NURBS patch as a Range of NURBS patches.
         * @return An array of NURBS curves representing the boundaries.
         */
        Range<NURBSPatch<vec_t, Vec<scalar_t, param_dim - 1>>> getBoundaries() const;

        /**
         * Get vector from the parametric domain of the NURBS patch from the
         * parameter of the NURBS patch representing one of the boundaries of the NURBS.
         * Patch parameter can be generated `epsilon` times the length of the domain away
         * from the boundary. This is useful when normals are not defined at the boundary of the patch.
         * @param t Parameter of the NURBS patch representing the boundary.
         * @param i Index of the boundary.
         * @param epsilon Epsilon.
         * @return Vector from the parametric domain of the NURBS patch corresponding
         * to the input parameter.
         */
        param_vec_t getPatchParameterFromBoundaryParameter(const Vec<scalar_t, param_dim - 1> &t, int i,
                                                           scalar_t epsilon = 0) const;

        /**
         * Get parameter from the parametric domain of the NURBS patch representing the
         * boundary of the NURBS patch corresponding to a parameter from the parametric domain
         * of the NURBS surface.
         * @param t Parameter of the NURBS patch.
         * @param i Index of the boundary.
         * @return Parameter from the parametric domain of the NURBS boundary curve corresponding
         * to the input parameter.
         */
        Vec<scalar_t, param_dim - 1> getBoundaryParameterFromPatchParameter(const param_vec_t &t,
                                                                            int i) const;
    };

} // namespace mm

#endif // MEDUSA_BITS_DOMAINS_NURBSPATCH_FWD_HPP_