#531 Utilities: add AreEqual (that at the moment uses up Yuni's implementation).

Sources/FiniteElement/Operators/GlobalVariationalOperator/GlobalVariationalOperator.hxx

@@ -556,7 +556,7 @@ namespace HappyHeart
 
         auto& local_matrix = elementary_data.GetNonCstMatrixResult();
         
-      //  Private::ReportValuesInGlobalMatrix531(linear_algebra_tuple);
+        Private::InjectInGlobalMatrix531(linear_algebra_tuple, local_operator);
         
         
         

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

@@ -9,8 +9,9 @@
 #ifndef _HAPPY_HEART__FINITE_ELEMENT__OPERATORS__GLOBAL_VARIATIONAL_OPERATOR__PRIVATE__IMPL__RECURSIVITY_HPP_
 # define _HAPPY_HEART__FINITE_ELEMENT__OPERATORS__GLOBAL_VARIATIONAL_OPERATOR__PRIVATE__IMPL__RECURSIVITY_HPP_
 
-# include <memory>
-# include <vector>
+# include "Utilities/Numeric/Numeric.hpp"
+
+# include "FiniteElement/FiniteElement/LocalFEltSpace.hpp"
 
 
 namespace HappyHeart
@@ -33,7 +34,10 @@ namespace HappyHeart
                 
                 static void Assembly531(const LinearAlgebraTupleT& linear_algebra_tuple);
                 
-                
+                template<class LocalVariationalOperatorT>
+                static void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
+                                                    const LocalFEltSpace& local_felt_space,
+                                                    LocalVariationalOperatorT& local_variational_operator);
 
                 
             };
@@ -47,6 +51,11 @@ namespace HappyHeart
                 
                 static void Assembly531(const LinearAlgebraTupleT& linear_algebra_tuple);
                 
+                template<class LocalVariationalOperatorT>
+                static void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
+                                                    const LocalFEltSpace& local_felt_space,
+                                                    LocalVariationalOperatorT& local_variational_operator);
+                
                 
             };
 

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

@@ -35,12 +35,70 @@ namespace HappyHeart
             }
             
             
+            
+            template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
+            template<class LocalVariationalOperatorT>
+            void Recursivity<LinearAlgebraTupleT, I, TupleSizeT>
+            ::InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
+                                      const LocalFEltSpace& local_felt_space,
+                                      LocalVariationalOperatorT& local_variational_operator)
+            {
+                static_assert(I < TupleSizeT, "Otherwise endless recursion!");
+                
+                auto& global_matrix = std::get<I>(linear_algebra_tuple).first;
+                
+                auto& elementary_data = local_variational_operator.GetNonCstElementaryData();
+                auto& local_matrix = elementary_data.GetNonCstMatrixResult();
+                
+                const double coefficient = std::get<I>(linear_algebra_tuple).second;
+                
+                if (!Utilities::AreEqual(coefficient, 1.))
+                    local_matrix *= coefficient;
+                
+                const auto& unknown = local_variational_operator.GetUnknown();
+                const auto& unknown_and_numbering_subset = elementary_data.GetRefFElt(unknown).GetUnknownAndNumberingSubset();
+                const auto& local2global = local_felt_space.GetFiniteElt(unknown_and_numbering_subset).template GetLocal2Global<MpiScale::program_wise>();
+                
+//                #ifndef NDEBUG
+//                assert(static_cast<int>(row_local_2_global.size()) == local_matrix.GetM());
+//                assert(static_cast<int>(col_local_2_global.size()) == local_matrix.GetN());
+//                
+//                assert(global_matrix.GetRowNumberingSubset() == row_unknown_and_numbering_subset.GetNumberingSubset());
+//                assert(global_matrix.GetColNumberingSubset() == col_unknown_and_numbering_subset.GetNumberingSubset());
+//                #endif // NDEBUG
+//
+//                
+//                global_matrix.SetValues(row_local_2_global,
+//                                        column_local_2_global,
+//                                        local_matrix.GetData(),
+//                                        ADD_VALUES,
+//                                        __FILE__, __LINE__);
+
+                 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);
+            }
+            
+            
             template<class LinearAlgebraTupleT, std::size_t TupleSizeT>
             void Recursivity<LinearAlgebraTupleT, TupleSizeT, TupleSizeT>
             ::Assembly531(const LinearAlgebraTupleT& linear_algebra_tuple)
             {
                 // End recursion.
             }
+            
+            
+            template<class LinearAlgebraTupleT, std::size_t TupleSizeT>
+            template<class LocalVariationalOperatorT>
+            void Recursivity<LinearAlgebraTupleT, TupleSizeT, TupleSizeT>
+            ::InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
+                                      const LocalFEltSpace& local_felt_space,
+                                      LocalVariationalOperatorT& local_variational_operator)
+            {
+                // End recursion.
+            }
 
            
             

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

@@ -27,8 +27,10 @@ namespace HappyHeart
         
         
         
-        template<class LinearAlgebraTupleT>
-        void InjectInGlobalMatrix(const LinearAlgebraTupleT& linear_algebra_tuple);
+        template<class LinearAlgebraTupleT, class LocalVariationalOperatorT>
+        void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
+                                     const LocalFEltSpace& local_felt_space,
+                                     LocalVariationalOperatorT& local_operator);
         
                 
         

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

@@ -26,6 +26,17 @@ namespace HappyHeart
             ::Assembly531(linear_algebra_tuple);
            
         }
+        
+        
+        template<class LinearAlgebraTupleT, class LocalVariationalOperatorT>
+        void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
+                                     LocalVariationalOperatorT& local_variational_operator)
+        {
+            
+            Impl::Recursivity<LinearAlgebraTupleT, 0, std::tuple_size<LinearAlgebraTupleT>::value>
+            ::InjectInGlobalMatrix531(linear_algebra_tuple, local_variational_operator);
+            
+        }
 
         
         

Sources/Utilities/Numeric/Numeric.hpp

@@ -17,6 +17,8 @@
 # include <cmath>
 # include <type_traits>
 
+# include "ThirdParty/IncludeWithoutWarning/Yuni/Core/Math.hpp"
+
 
 namespace HappyHeart
 {
@@ -72,6 +74,24 @@ namespace HappyHeart
         IsZero(T value, T epsilon = DefaultEpsilon<T>()) noexcept;
         
         
+        /*!
+         * \brief Check whether a value is close enough to another.
+         *
+         * \tparam T Floating point type considered (float, double or long double).
+         *
+         * \param[in] lhs Lhs value.
+         * \param[in] rhs Rhs value.
+         *
+         * \internal Currently Yuni's function is called; in due time an HappyHeart one should be given (to be able to 
+         * sever Yuni dependancy).
+         *
+         */
+        template<class T>
+        constexpr std::enable_if_t<std::is_floating_point<T>::value, bool>
+        AreEqual(T lhs, T rhs) noexcept;
+
+        
+        
         
     } // namespace Utilities
     

Sources/Utilities/Numeric/Numeric.hxx

@@ -48,6 +48,13 @@ namespace HappyHeart
         }
         
         
+        template<class T>
+        constexpr std::enable_if_t<std::is_floating_point<T>::value, bool> AreEqual(T lhs, T rhs) noexcept
+        {
+            return Yuni::Math::Equals(lhs, rhs);
+        }
+        
+        
         
     } // namespace Utilities