//! \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::value. */ template struct Recursivity { //! Type of the current element; might be either GlobalMatrixWithCoefficient or GlobalVectorWithCoefficient. using current_type = typename std::tuple_element::type; // Check here the type of current item is one of those expected. static_assert(std::is_same::value || std::is_same::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 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 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 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 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 struct ZeroSpecialCase { //! Type of the current element; might be either GlobalMatrixWithCoefficient or GlobalVectorWithCoefficient. using current_type = typename std::tuple_element::type; /*! * \brief This static method computes the coefficient required by * \a Recursivity::InjectInGlobalLinearAlgebraImpl. * * See this function for a detailed explanation. */ template static double FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple); }; // ============================ // Stop points of recursive loops. // ============================ //! \cond IGNORE_BLOCK_IN_DOXYGEN template struct Recursivity { static void Assembly(const LinearAlgebraTupleT& linear_algebra_tuple); template static void InjectInGlobalLinearAlgebra(const LinearAlgebraTupleT& linear_algebra_tuple, const LocalFEltSpace& local_felt_space, LocalVariationalOperatorT& local_variational_operator); }; template struct ZeroSpecialCase { template 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_