//! \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_
*