#531 Fewer scale of local matrices; result remains the same.

Sources/FiniteElement/Operators/GlobalVariationalOperator/Private/Impl/Recursivity.hpp

@@ -38,7 +38,7 @@ namespace HappyHeart
                 static void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
                                                     const LocalFEltSpace& local_felt_space,
                                                     LocalVariationalOperatorT& local_variational_operator);
-
+                
                 
             };
             
@@ -61,6 +61,29 @@ namespace HappyHeart
 
             
             
+            template<class LinearAlgebraTupleT, std::size_t I>
+            struct ZeroSpecialCase
+            {
+                
+                
+                static double FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple);
+                
+                
+            };
+            
+            
+            
+            template<class LinearAlgebraTupleT>
+            struct ZeroSpecialCase<LinearAlgebraTupleT, 0>
+            {
+                
+                
+                static double FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple);
+                
+                
+            };
+            
+            
             
         } //namespace Impl
         

Sources/FiniteElement/Operators/GlobalVariationalOperator/Private/Impl/Recursivity.hxx

@@ -34,7 +34,7 @@ namespace HappyHeart
                 Recursivity<LinearAlgebraTupleT, I + 1, TupleSizeT>::Assembly531(linear_algebra_tuple);
             }
             
-            
+
             
             template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
             template<class LocalVariationalOperatorT>
@@ -54,8 +54,19 @@ namespace HappyHeart
                 
                 const double coefficient = std::get<I>(linear_algebra_tuple).second;
                 
-                if (!Utilities::AreEqual(coefficient, 1.))
-                    local_matrix *= coefficient;
+                const double previous_coefficient =
+                    ZeroSpecialCase<LinearAlgebraTupleT, I>::FetchPreviousCoefficient(linear_algebra_tuple);
+                
+                assert(!Utilities::IsZero(previous_coefficient));
+                
+                
+                
+                // Trick here: local matrix was multiplied by a coefficient in previous iteration; we here both apply
+                // the new one and cancel the previous one.
+                const auto factor = coefficient / previous_coefficient;
+                
+                if (!Utilities::AreEqual(factor, 1.))
+                    local_matrix *= factor;
                 
                 const auto& row_numbering_subset = global_matrix.GetRowNumberingSubset();
                 const auto& col_numbering_subset = global_matrix.GetColNumberingSubset();
@@ -96,11 +107,6 @@ namespace HappyHeart
                     
                 }
                 
-
-                 if (!Utilities::AreEqual(coefficient, 1.)) // \todo #531 Be clever and get rid of it (use value from previous iteration in the first product).
-                    local_matrix *= 1. / coefficient;
-
-                
                 Recursivity<LinearAlgebraTupleT, I + 1, TupleSizeT>::InjectInGlobalMatrix531(linear_algebra_tuple,
                                                                                              local_felt_space,
                                                                                              local_variational_operator);
@@ -130,6 +136,25 @@ namespace HappyHeart
             }
 
            
+            template<class LinearAlgebraTupleT>
+            double ZeroSpecialCase<LinearAlgebraTupleT, 0>
+            ::FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple)
+            {
+                static_cast<void>(linear_algebra_tuple);
+                return 1.;
+            }
+            
+            
+            
+            template<class LinearAlgebraTupleT, std::size_t I>
+            double ZeroSpecialCase<LinearAlgebraTupleT, I>
+            ::FetchPreviousCoefficient(const LinearAlgebraTupleT& linear_algebra_tuple)
+            {
+                return std::get<I-1>(linear_algebra_tuple).second;
+            }
+
+            
+            
             
         } //namespace Impl