#531 Fix the non monolithic case; efficiency has been let aside at the moment.

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

@@ -64,7 +64,13 @@ namespace HappyHeart
                                                                LocalVariationalOperatorT& local_variational_operator,
                                                                const GlobalVectorWithCoefficient& global_vector,
                                                                double previous_coefficient);
-
+                
+                
+                // \todo #531 Improve efficiency (copy is stupid here...)
+                template<class LocalVariationalOperatorT>
+                static LocalMatrix ComputeLocalMatrix(const LocalVariationalOperatorT& local_variational_operator,
+                                                      const NumberingSubset& row_numbering_subset,
+                                                      const NumberingSubset& col_numbering_subset);
                 
             };
             

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

@@ -84,7 +84,15 @@ namespace HappyHeart
                 auto& global_matrix = global_matrix_with_coefficient.first;
                 
                 auto& elementary_data = local_variational_operator.GetNonCstElementaryData();
-                auto& local_matrix = elementary_data.GetNonCstMatrixResult();
+//                auto& local_matrix = elementary_data.GetNonCstMatrixResult();
+                
+                
+                const auto& row_numbering_subset = global_matrix.GetRowNumberingSubset();
+                const auto& col_numbering_subset = global_matrix.GetColNumberingSubset();
+                
+                auto local_matrix = ComputeLocalMatrix(local_variational_operator,
+                                                       row_numbering_subset,
+                                                       col_numbering_subset);
                 
                 // \todo #531 We might need here to reduce the local matrix!
                 
@@ -98,9 +106,7 @@ namespace HappyHeart
                 
                 if (!Utilities::AreEqual(factor, 1.))
                     local_matrix *= factor;
-                
-                const auto& row_numbering_subset = global_matrix.GetRowNumberingSubset();
-                const auto& col_numbering_subset = global_matrix.GetColNumberingSubset();
+              
                 
                 const auto& unknown_and_numbering_subset_list =
                     local_variational_operator.GetUnknownAndNumberingSubsetList();
@@ -140,6 +146,69 @@ namespace HappyHeart
             }
             
             
+            template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
+            template<class LocalVariationalOperatorT>
+            LocalMatrix Recursivity<LinearAlgebraTupleT, I, TupleSizeT>
+            ::ComputeLocalMatrix(const LocalVariationalOperatorT& local_variational_operator,
+                                 const NumberingSubset& row_numbering_subset,
+                                 const NumberingSubset& col_numbering_subset)
+            {
+                const auto& elementary_data = local_variational_operator.GetElementaryData();
+                
+                const auto& unknown_and_numbering_subset_list =
+                    local_variational_operator.GetUnknownAndNumberingSubsetList();
+                
+
+                Unknown::vector_const_shared_ptr row_unknown_list;
+                Unknown::vector_const_shared_ptr col_unknown_list;
+                
+                for (const auto& unknown_and_numbering_subset_ptr : unknown_and_numbering_subset_list)
+                {
+                    assert(!(!unknown_and_numbering_subset_ptr));
+                    const auto& unknown_and_numbering_subset = *unknown_and_numbering_subset_ptr;
+                    
+                    if (unknown_and_numbering_subset.GetNumberingSubset() == row_numbering_subset)
+                        row_unknown_list.push_back(unknown_and_numbering_subset.GetUnknownPtr());
+                    
+                    if (unknown_and_numbering_subset.GetNumberingSubset() == col_numbering_subset)
+                        col_unknown_list.push_back(unknown_and_numbering_subset.GetUnknownPtr());
+                }
+                
+                // \todo #513 Here I assumed clearly exactly one unknown is expected for row and columm!
+                assert(row_unknown_list.size() == 1 && "See #513");
+                assert(col_unknown_list.size() == 1 && "See #513");
+                const auto& unknown1 = *(row_unknown_list.front());
+                const auto& unknown2 = *(col_unknown_list.front());
+                
+//                auto it = reduced_local_matrix_.find(elementary_data.GetRefGeomEltPtr());
+//                assert(it != reduced_local_matrix_.cend());
+//                `auto& ret = it->second;
+                
+                // #197 const_cast to remove! Very ugly, but I guarantee the underlying value is not modified...
+                auto& local_matrix = const_cast<LocalMatrix&>(elementary_data.GetMatrixResult());
+            
+                auto block = elementary_data.ExtractBlockFromLocalMatrix(elementary_data.GetRefFElt(unknown1),
+                                                                         elementary_data.GetRefFElt(unknown2),
+                                                                         local_matrix);
+                
+                const int Nrow_reduced = block.GetM();
+                const int Ncol_reduced = block.GetN();
+                
+                assert(Nrow_reduced <= local_matrix.GetM());
+                assert(Ncol_reduced <= local_matrix.GetN());
+                
+                LocalMatrix ret(Nrow_reduced, Ncol_reduced);
+                
+                for (int i = 0; i < Nrow_reduced; ++i)
+                    for (int j = 0; j < Ncol_reduced; ++j)
+                        ret(i, j) = block(i, j);
+                
+                return ret;
+            }
+            
+
+            
+            
             template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
             template<class LocalVariationalOperatorT>
             void Recursivity<LinearAlgebraTupleT, I, TupleSizeT>
@@ -237,6 +306,7 @@ namespace HappyHeart
             
             
             
+            
         } //namespace Impl
         
         

Sources/FiniteElement/RefFiniteElement/Private/ElementaryDataImpl.hpp

@@ -148,8 +148,6 @@ namespace HappyHeart
             
             static BlockLocalVector ExtractBlockFromLocalVector(const RefFEltInLocalNumbering& ref_felt,
                                                                 LocalVector& vector);
-
-
             
             const QuadratureRule& GetQuadratureRule() const;
             

Sources/FiniteElement/Unknown/UnknownAndNumberingSubset.hpp

@@ -60,6 +60,9 @@ namespace HappyHeart
         //! Get the underlying unknown.
         const Unknown& GetUnknown() const;
         
+        //! Get the underlying unknown.
+        const Unknown::const_shared_ptr& GetUnknownPtr() const;
+        
         //! Get the underlying numbering subset.
         const NumberingSubset& GetNumberingSubset() const;
         

Sources/FiniteElement/Unknown/UnknownAndNumberingSubset.hxx

@@ -20,6 +20,13 @@ namespace HappyHeart
         return *unknown_;
     }
     
+    
+    inline const Unknown::const_shared_ptr& UnknownAndNumberingSubset::GetUnknownPtr() const
+    {
+        assert(!(!unknown_));
+        return unknown_;
+    }
+    
 
     inline const NumberingSubset& UnknownAndNumberingSubset::GetNumberingSubset() const
     {