//! \file 
//
//
//  Recursivity.hpp
//  HappyHeart
//
//  Created by Sebastien Gilles on 29/04/15.
//  Copyright (c) 2015 Inria. All rights reserved.
//

#ifndef HAPPY_HEART_x_OPERATORS_x_GLOBAL_VARIATIONAL_OPERATOR_x_PRIVATE_x_IMPL_x_RECURSIVITY_HPP_
# define HAPPY_HEART_x_OPERATORS_x_GLOBAL_VARIATIONAL_OPERATOR_x_PRIVATE_x_IMPL_x_RECURSIVITY_HPP_

# include "Utilities/Numeric/Numeric.hpp"

# include "Core/LinearAlgebra/GlobalMatrix.hpp"
# include "Core/LinearAlgebra/GlobalVector.hpp"

# include "FiniteElement/FiniteElement/LocalFEltSpace.hpp"


namespace HappyHeart
{
    
    
    namespace Private
    {
        
        
        namespace Impl
        {
            
            
            /*!
             * \brief Helper class used for metaprogramming iterations.
             *
             * The main purpose here is to apply the same operation to each item of \a LinearAlgebraTupleT.
             *
             * \tparam LinearAlgebraTupleT A tuple that may include either \a GlobalMatrixWithCoefficient (for bilinear
             * operators) or \a GlobalVectorWithCoefficient (for linear operators) objects. Some non linear operators
             * may include both types of objects; the ordering doesn't matter.
             * \tparam I Current element of the tuple.
             * \tparam TupleSizeT The result of std::tuple_size<LinearAlgebraTupleT>::value.
             */
            template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
            struct Recursivity
            {
                
                //! Type of the current element; might be either GlobalMatrixWithCoefficient or GlobalVectorWithCoefficient.
                using current_type = typename std::tuple_element<I, LinearAlgebraTupleT>::type;
                
                // Check here the type of current item is one of those expected.
                static_assert(std::is_same<current_type, GlobalMatrixWithCoefficient&>::value
                              || std::is_same<current_type, GlobalVectorWithCoefficient&>::value
                              , "Tuple is expected to include only one of those types.");

                
                /*!
                 * \brief Call Assembly() on current element and then apply recursion.
                 */
                static void Assembly(const LinearAlgebraTupleT& linear_algebra_tuple);
                
                /*!
                 * \brief Report the result of elementary computation into the matrix pointed by current item and
                 * then apply recursion.
                 *
                 * \param[in] local_felt_space Local finite element space being assembled. Remind a local finite element
                 * space is the object linked to a geometric element that hold finite elements information.
                 * \param[in,out] local_variational_operator Local variational operator in charge of the computation of 
                 * the local linear algebra. It also holds the results of these computations, hence its presence here.
                 */
                template<class LocalVariationalOperatorT>
                static void InjectInGlobalLinearAlgebra(const LinearAlgebraTupleT& linear_algebra_tuple,
                                                        const LocalFEltSpace& local_felt_space,
                                                        LocalVariationalOperatorT& local_variational_operator);
                
            private:

                
                /*!
                 * \brief Private function called by \a InjectInGlobalLinearAlgebra when current element refers to a 
                 * matrix.
                 *
                 * \param[in] previous_coefficient The elementary matrix has been computed once but during each assembling
                 * the coefficient applied is not the same. Unfortunately Petsc doesn't provide a function that 
                 * set the values multiplied by a factor, so I have to multiply the local matrix by a factor.
                 *
                 * Let's say for instance the operator is assembled in two matrices (this is pseudo code, not actual
                 * HappyHeart call):
                 *
                 * \code
                 * operator.Assemble({ Matrix1, factor1 }, {Matrix2, factor2});
                 * \endcode
                 *
                 * One way to do injection is (still pseudo code):
                 * \code
                 * local_matrix *= factor1;
                 * Inject local_matrix into  Matrix1;
                 * local_matrix /= factor1; // of course with check no 0 here...
                 * local_matrix *= factor2;
                 * Inject local_matrix into  Matrix2;
                 * local_matrix /= factor2; // of course with check no 0 here...
                 * \endcode
                 *
                 * However we see here that on one hand the last division isn't useful, on the other hand two of the 
                 * operations could be done in one step.
                 *
                 * The latter case is exactly the point of \a previous_coefficient: when assembling into Matrix2
                 * the program looks the previous matrix being assembled to determine the correct factor. So in
                 * pseudo-code what is done is:
                 *
                 * \code
                 * local_matrix *= factor1;
                 * Inject local_matrix into  Matrix1;
                 * local_matrix *= factor2 / factor1; of course with check no 0 here...
                 * Inject local_matrix into  Matrix2;
                 * \endcode
                 *
                 * If there are both matrices and vectors in \a LinearAlgebraTupleT, it still works whatever the
                 * ordering: when computing a matrix vectors objects are skipped in the determination of the 
                 * previous factor.
                 *
                 */
                template<class LocalVariationalOperatorT>
                static void InjectInGlobalLinearAlgebraImpl(const LocalFEltSpace& local_felt_space,
                                                               LocalVariationalOperatorT& local_variational_operator,
                                                               const GlobalMatrixWithCoefficient& global_matrix,
                                                               double previous_coefficient);

                //! Same as previously for vectors.
                template<class LocalVariationalOperatorT>
                static void InjectInGlobalLinearAlgebraImpl(const LocalFEltSpace& local_felt_space,
                                                               LocalVariationalOperatorT& local_variational_operator,
                                                               const GlobalVectorWithCoefficient& global_vector,
                                                               double previous_coefficient);
                
                
                /*!
                 * \brief Reduce the elementary matrix to the subset that must be reported in the global matrix.
                 *
                 * Let's consider for instance a Stokes operator in a non monolithic way: same elementary matrix
                 * must be reported in two blocks of different sizes (said size is deduced from numbering subset 
                 * arguments).
                 * Role of current function is to generate a smaller matrix that includes only the elements to report
                 * into the bigger one. Reason for this is the prototype of Petsc function: it expects the elements
                 * to be set to be contiguous in memory.
                 */
                template<class LocalVariationalOperatorT>
                static LocalMatrix& ComputeReducedLocalMatrix(const NumberingSubset& row_numbering_subset,
                                                              const NumberingSubset& col_numbering_subset,
                                                              LocalVariationalOperatorT& local_variational_operator);
                
            };
            
            
            /*!
             * \brief Another struct for recursion that works in the opposite direction: the stopping condition
             * is when index is 0.
             *
             *
             * \tparam LinearAlgebraTupleT A tuple that may include either \a GlobalMatrixWithCoefficient (for bilinear
             * operators) or \a GlobalVectorWithCoefficient (for linear operators) objects. Some non linear operators
             * may include both types of objects; the ordering doesn't matter.
             * \tparam I Current element of the tuple.
             */
            template<class LinearAlgebraTupleT, std::size_t I>
            struct ZeroSpecialCase
            {
                //! Type of the current element; might be either GlobalMatrixWithCoefficient or GlobalVectorWithCoefficient.
                using current_type = typename std::tuple_element<I, LinearAlgebraTupleT>::type;
          

                /*!
                 * \brief This static method computes the coefficient required by 
                 * \a Recursivity::InjectInGlobalLinearAlgebraImpl.
                 *
                 * See this function for a detailed explanation.
                 */
                template<class CurrentTypeT>
                static double FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple);
                
                
            };
            
            
            
            
            // ============================
            // Stop points of recursive loops.
            // ============================
            //! \cond IGNORE_BLOCK_IN_DOXYGEN
            
            
            template<class LinearAlgebraTupleT, std::size_t TupleSizeT>
            struct Recursivity<LinearAlgebraTupleT, TupleSizeT, TupleSizeT>
            {
                
                
                static void Assembly(const LinearAlgebraTupleT& linear_algebra_tuple);
                
                template<class LocalVariationalOperatorT>
                static void InjectInGlobalLinearAlgebra(const LinearAlgebraTupleT& linear_algebra_tuple,
                                                        const LocalFEltSpace& local_felt_space,
                                                        LocalVariationalOperatorT& local_variational_operator);
                
                
            };
            
            
            template<class LinearAlgebraTupleT>
            struct ZeroSpecialCase<LinearAlgebraTupleT, 0>
            {
                
                template<class CurrentTypeT>
                static double FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple);
                
                
            };
            
            //! \endcond IGNORE_BLOCK_IN_DOXYGEN
            // ============================
            // End of stop points of recursive loops.
            // ============================
            
            
            
            
            
        } // namespace Impl
        
        
    } // namespace Private
    
    
} // namespace HappyHeart


# include "Operators/GlobalVariationalOperator/Private/Impl/Recursivity.hxx"


#endif // HAPPY_HEART_x_OPERATORS_x_GLOBAL_VARIATIONAL_OPERATOR_x_PRIVATE_x_IMPL_x_RECURSIVITY_HPP_