Recursivity.hpp 11.2 KB
Newer Older
1
2
//! \file 
//
3
4
5
6
7
8
9
10
//
//  Recursivity.hpp
//  HappyHeart
//
//  Created by Sebastien Gilles on 29/04/15.
//  Copyright (c) 2015 Inria. All rights reserved.
//

11
12
#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_
13

14
15
# include "Utilities/Numeric/Numeric.hpp"

16
17
18
# include "Core/LinearAlgebra/GlobalMatrix.hpp"
# include "Core/LinearAlgebra/GlobalVector.hpp"

19
# include "FiniteElement/FiniteElement/LocalFEltSpace.hpp"
20
21
22
23
24
25
26
27
28
29
30
31
32
33


namespace HappyHeart
{
    
    
    namespace Private
    {
        
        
        namespace Impl
        {
            
            
34
35
36
37
38
            /*!
             * \brief Helper class used for metaprogramming iterations.
             *
             * The main purpose here is to apply the same operation to each item of \a LinearAlgebraTupleT.
             *
39
40
             * \tparam LinearAlgebraTupleT A tuple that may include either \a GlobalMatrixWithCoefficient (for bilinear
             * operators) or \a GlobalVectorWithCoefficient (for linear operators) objects. Some non linear operators
41
42
43
44
             * 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.
             */
45
            template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
GILLES Sebastien's avatar
GILLES Sebastien committed
46
            struct Recursivity
47
48
            {
                
49
                //! Type of the current element; might be either GlobalMatrixWithCoefficient or GlobalVectorWithCoefficient.
50
                using current_type = typename std::tuple_element<I, LinearAlgebraTupleT>::type;
51
52
                
                // Check here the type of current item is one of those expected.
53
54
                static_assert(std::is_same<current_type, GlobalMatrixWithCoefficient&>::value
                              || std::is_same<current_type, GlobalVectorWithCoefficient&>::value
55
                              , "Tuple is expected to include only one of those types.");
56

57
                
58
59
60
                /*!
                 * \brief Call Assembly() on current element and then apply recursion.
                 */
61
                static void Assembly(const LinearAlgebraTupleT& linear_algebra_tuple);
GILLES Sebastien's avatar
GILLES Sebastien committed
62
                
63
64
65
66
67
68
                /*!
                 * \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.
69
70
                 * \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.
71
                 */
72
                template<class LocalVariationalOperatorT>
73
                static void InjectInGlobalLinearAlgebra(const LinearAlgebraTupleT& linear_algebra_tuple,
74
75
                                                        const LocalFEltSpace& local_felt_space,
                                                        LocalVariationalOperatorT& local_variational_operator);
76
                
77
78
79
            private:

                
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
                /*!
                 * \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.
                 *
                 */
124
                template<class LocalVariationalOperatorT>
125
                static void InjectInGlobalLinearAlgebraImpl(const LocalFEltSpace& local_felt_space,
126
127
128
129
                                                               LocalVariationalOperatorT& local_variational_operator,
                                                               const GlobalMatrixWithCoefficient& global_matrix,
                                                               double previous_coefficient);

130
                //! Same as previously for vectors.
131
                template<class LocalVariationalOperatorT>
132
                static void InjectInGlobalLinearAlgebraImpl(const LocalFEltSpace& local_felt_space,
133
134
135
                                                               LocalVariationalOperatorT& local_variational_operator,
                                                               const GlobalVectorWithCoefficient& global_vector,
                                                               double previous_coefficient);
136
137
                
                
138
139
140
                /*!
                 * \brief Reduce the elementary matrix to the subset that must be reported in the global matrix.
                 *
141
                 * Let's consider for instance a Stokes operator in a non monolithic way: same elementary matrix
142
143
144
145
146
147
                 * 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.
                 */
148
                template<class LocalVariationalOperatorT>
149
                static LocalMatrix& ComputeReducedLocalMatrix(const NumberingSubset& row_numbering_subset,
150
151
                                                              const NumberingSubset& col_numbering_subset,
                                                              LocalVariationalOperatorT& local_variational_operator);
152
153
154
155
                
            };
            
            
156
157
158
159
160
            /*!
             * \brief Another struct for recursion that works in the opposite direction: the stopping condition
             * is when index is 0.
             *
             *
161
162
             * \tparam LinearAlgebraTupleT A tuple that may include either \a GlobalMatrixWithCoefficient (for bilinear
             * operators) or \a GlobalVectorWithCoefficient (for linear operators) objects. Some non linear operators
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
             * 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
            
193
194
            
            template<class LinearAlgebraTupleT, std::size_t TupleSizeT>
GILLES Sebastien's avatar
GILLES Sebastien committed
195
            struct Recursivity<LinearAlgebraTupleT, TupleSizeT, TupleSizeT>
196
197
198
            {
                
                
199
                static void Assembly(const LinearAlgebraTupleT& linear_algebra_tuple);
200
                
201
                template<class LocalVariationalOperatorT>
202
                static void InjectInGlobalLinearAlgebra(const LinearAlgebraTupleT& linear_algebra_tuple,
203
204
                                                        const LocalFEltSpace& local_felt_space,
                                                        LocalVariationalOperatorT& local_variational_operator);
205
                
206
207
208
209
                
            };
            
            
210
211
            template<class LinearAlgebraTupleT>
            struct ZeroSpecialCase<LinearAlgebraTupleT, 0>
212
213
            {
                
214
                template<class CurrentTypeT>
215
216
217
218
219
                static double FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple);
                
                
            };
            
220
221
222
223
            //! \endcond IGNORE_BLOCK_IN_DOXYGEN
            // ============================
            // End of stop points of recursive loops.
            // ============================
224
225
226
227
            
            
            
            
228
            
229
        } // namespace Impl
230
231
        
        
232
    } // namespace Private
233
234
235
236
237
    
    
} // namespace HappyHeart


238
# include "Operators/GlobalVariationalOperator/Private/Impl/Recursivity.hxx"
239
240


241
#endif // HAPPY_HEART_x_OPERATORS_x_GLOBAL_VARIATIONAL_OPERATOR_x_PRIVATE_x_IMPL_x_RECURSIVITY_HPP_