#820 Add missing part of Darcy operator. Add also a missing ZeroEntries().

Sources/Core/LinearAlgebra/GlobalVector.hpp

@@ -13,6 +13,7 @@
 
 # include <memory>
 # include <vector>
+# include <array>
 
 # include "ThirdParty/Wrappers/Petsc/Vector/Vector.hpp"
 
@@ -55,6 +56,10 @@ namespace HappyHeart
         //! Alias to vector of unique pointers.
         using vector_unique_ptr = std::vector<unique_ptr>;
         
+        //! Alias to array of unique pointers.
+        template<std::size_t N>
+        using array_unique_ptr = std::array<unique_ptr, N>;
+        
         //! Alias to parent.
         using parent = Wrappers::Petsc::Vector;
     

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx

@@ -70,7 +70,7 @@ namespace HappyHeart
                         // is probably not required publicly in the process; we could probably compute directly contrib to monolithic
                         // from the original VelSHalfVhf.
                         auto& vel_contrib_on_vel_ns = variable_holder.GetVelSHalfVhfZeroedOutOfRobinInterface();
-                        auto& vel_contrib_on_monolithic = variable_holder.GetNonCstWorkFluidMonolithic();
+                        auto& vel_contrib_on_monolithic = variable_holder.template GetNonCstWorkFluidMonolithic<0>();
                         Wrappers::Petsc::MatMult(interpolator_holder.GetExpandVelocity().GetInterpolationMatrix(),
                                                  vel_contrib_on_vel_ns,
                                                  vel_contrib_on_monolithic,
@@ -518,7 +518,6 @@ namespace HappyHeart
                         
                         variable_holder.template DevPrint<DevPhase::newton_fp>(contribution,
                                                                           "pressure_contrib_to_rhs");
-                        
                     }
                     
                     {
@@ -555,18 +554,38 @@ namespace HappyHeart
                             
                             Wrappers::Petsc::MatMult(interpolator_holder.GetMass2Monolithic().GetInterpolationMatrix(),
                                                      contrib_to_rhs,
-                                                     contribution, // replace buf_for_test by rhs!
+                                                     contribution,
                                                      __FILE__, __LINE__);
+                            
+                            variable_holder.template DevPrint<DevPhase::newton_fp>(contribution,
+                                                                                   "fluid_mass_contrib_to_rhs_matrix_part");
                         }
                         
+                        
                         {
-                            GlobalVectorWithCoefficient contribution_with_coeff(contribution, 1.); // replace buf_for_test by rhs!
+                            auto& vel_term = variable_holder.template GetNonCstWorkFluidMonolithic<0>();
+                            vel_term.Copy(GetCorrectedValues(), __FILE__, __LINE__);
+                            
+                            auto& vel_SHalfVhf_contrib = variable_holder.template GetNonCstWorkFluidMonolithic<1>();
+                            
+                            Wrappers::Petsc::MatMult(interpolator_holder.GetExpandVelocity().GetInterpolationMatrix(),
+                                                     variable_holder.GetVelSHalfVhf(),
+                                                     vel_SHalfVhf_contrib,
+                                                     __FILE__, __LINE__);
+                                                     
+                            
+                            Wrappers::Petsc::AXPY(-1.,
+                                                  vel_SHalfVhf_contrib,
+                                                  vel_term,
+                                                  __FILE__, __LINE__);
+                            
+                            GlobalVectorWithCoefficient contribution_with_coeff(contribution, 1.);
                             GetHybridVectorOperator().Assemble(std::make_tuple(std::ref(contribution_with_coeff)),
-                                                               GetCorrectedValues());
+                                                               vel_term);
                         }
                         
                         variable_holder.template DevPrint<DevPhase::newton_fp>(contribution,
-                                                                          "fluid_mass_contrib_to_rhs");
+                                                                               "fluid_mass_contrib_to_rhs");
 
                         
                         Wrappers::Petsc::AXPY(1.,

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/LocalVariationalOperatorInstances/HybridVector.cpp

@@ -104,6 +104,9 @@ namespace HappyHeart
                             {
                                 // Equivalent in Freefem:
                                 // - int2d(ThF)( weighting(x) * rhoF * divPhi(velFr,velFz,phi) * qtest )
+                                // + int2d(ThF)( weighting(x) * rhoF * divPhi(velSHalfVhfr,velSHalfVhfz,phi) * qtest )
+                                // NOTE: it is velFr - velSHalfVhf that is givemn to Assemble (\todo #820 document that!!!)
+
 
                                 assert(vel_dof_index < local_velocity.size());
                                 
@@ -119,7 +122,10 @@ namespace HappyHeart
                     }
                     // \todo #820 Still missing.
                     
-                    // int2d(ThF)( weighting(x) * rhoF * divPhi(velSHalfVhfr,velSHalfVhfz,phi) * qtest )
+                    
+                    // Note: term
+                    
+                    
                     // int2d(ThF)( weighting(x) * beta * rhoF * pOld * qtest) ignored as beta = 0
 
                     

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx

@@ -149,12 +149,12 @@ namespace HappyHeart
                 
                 {
                     decltype(auto) darcy_vector = GetNonCstDarcyVector();
+                    darcy_vector.ZeroEntries(__FILE__, __LINE__);
                     
                     inner_varf.ComputeDarcy(-1., IsFullDarcy::no, darcy_vector);
                     
-                    
                     variable_holder.template DevPrint<DevPhase::none>(darcy_vector,
-                                                             "darcy_vector_after_FP_mixt");
+                                                                      "darcy_vector_after_newton_FP_mixt");
                 }
             }
             

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp

@@ -121,8 +121,18 @@ namespace HappyHeart
             }
             
             {
-                work_fluid_monolithic_ = std::make_unique<GlobalVector>(fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure)));
-                AllocateGlobalVector(fluid_god_of_dof, GetNonCstWorkFluidMonolithic());
+                constexpr auto size = Utilities::ArraySize<decltype(work_fluid_monolithic_)>::GetValue();
+                
+                static_assert(size > 0ul, "");
+                
+                work_fluid_monolithic_[0] = std::make_unique<GlobalVector>(fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure)));
+                
+                auto& first_item = GetNonCstWorkFluidMonolithic<0>();
+                AllocateGlobalVector(fluid_god_of_dof, first_item);
+                
+                for (auto i = 1ul; i < size; ++i)
+                    work_fluid_monolithic_[i] = std::make_unique<GlobalVector>(first_item);
+                
             }
             
         }

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp

@@ -14,6 +14,8 @@
 # include <memory>
 # include <vector>
 
+# include "Utilities/Containers/Array.hpp"
+
 # include "Core/LinearAlgebra/GlobalMatrix.hpp"
 # include "Core/LinearAlgebra/GlobalVector.hpp"
 
@@ -163,6 +165,7 @@ namespace HappyHeart
             GlobalVector& GetNonCstWorkFluidVelocity() const noexcept;
             
             //! Non constant accessor to the \#820 [DEV] Vector defined on monolithic numbering subset.
+            template<std::size_t I>
             GlobalVector& GetNonCstWorkFluidMonolithic() const noexcept;
             
             
@@ -234,11 +237,12 @@ namespace HappyHeart
         private:
             
             
-            //! \#820 [DEV] Vector defined on fluid velocity numbering subset.
-            mutable GlobalVector::unique_ptr work_fluid_velocity_ = nullptr;
+            //! Vector defined on fluid monolithic numbering subset.
+            mutable GlobalVector::array_unique_ptr<2> work_fluid_monolithic_
+                = Utilities::NullptrArray<GlobalVector::unique_ptr, 2>();
             
-            //! \#820 [DEV] Vector defined on fluid monolithic numbering subset.
-            mutable GlobalVector::unique_ptr work_fluid_monolithic_ = nullptr;
+            //! Vector defined on fluid velocity numbering subset.
+            mutable GlobalVector::unique_ptr work_fluid_velocity_ = nullptr;
          
         };
         

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx

@@ -181,10 +181,14 @@ namespace HappyHeart
         }
         
         
+        template<std::size_t I>
         inline GlobalVector& VariableHolder::GetNonCstWorkFluidMonolithic() const noexcept
         {
-            assert(!(!work_fluid_monolithic_));
-            return *work_fluid_monolithic_;
+            static_assert(I < Utilities::ArraySize<decltype(work_fluid_monolithic_)>::GetValue(),
+                          "Index should be in range!");
+
+            assert(!(!work_fluid_monolithic_[I]));
+            return *work_fluid_monolithic_[I];
         }