#820 Compute deltaDarcyVector.

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

@@ -266,9 +266,12 @@ namespace HappyHeart
                      */
                     const GlobalVariationalOperatorNS::Darcy& GetDarcyOperator() const noexcept;
                     
-                    //! Constant accessor to differential Darcy operator.
+                    //! Constant accessor to differential Darcy operator (first one in dH.edp)
                     const GlobalVariationalOperatorNS::DifferentialDarcy& GetDifferentialDarcyOperator() const noexcept;
                     
+                    //! Constant accessor to differential Darcy operator (second one in dH.edp)
+                    const GlobalVariationalOperatorNS::Darcy& GetDeltaDarcyOperator() const noexcept;
+                    
                     //! Constant accessor to the hybrid vector operator (same name as in Freefem at the moment...
                     const GlobalVariationalOperatorNS::HybridVector& GetHybridVectorOperator() const noexcept;
                     
@@ -480,6 +483,9 @@ namespace HappyHeart
                     //! Constant accessor to the parameter that encapsulates velocity part of the current solution.
                     const ParameterAtDof<ParameterNS::Type::vector>::type& GetVelocitySolutionParam() const noexcept;
                     
+                    //! Constant accessor to the parameter that encapsulates delta velocity part of the current solution (in dH.edp) // \todo #820 Improve comment.
+                    const ParameterAtDof<ParameterNS::Type::vector>::type& GetDeltaVelocitySolutionParam() const noexcept;
+                    
                     //! Parameter that encapsulates velSHalfVhf.
                     const ParameterAtDof<ParameterNS::Type::vector>::type& GetVelSHalfVhfParam() const noexcept;
                     
@@ -561,6 +567,12 @@ namespace HappyHeart
                     //! Non constant accessor to the deltaMixtVariables (	odo #820 Rename and comment!).
                     GlobalVector& GetNonCstDeltaMixtVariables() noexcept;
                     
+                    //! Constant accessor to the deltaDarcyVector (	odo #820 Rename and comment!).
+                    const GlobalVector& GetDeltaDarcyVector() const noexcept;
+                    
+                    //! Non constant accessor to the deltaDarcyVector (	odo #820 Rename and comment!).
+                    GlobalVector& GetNonCstDeltaDarcyVector() noexcept;
+                    
                     ///@}
                     
                     
@@ -680,9 +692,12 @@ namespace HappyHeart
                      */
                     GlobalVariationalOperatorNS::Darcy::const_unique_ptr darcy_operator_ = nullptr;
                     
-                    //! Differential Darcy operator.
+                    //! Differential Darcy operator. // \todo #820 Rename; it is the first one used in dH.
                     GlobalVariationalOperatorNS::DifferentialDarcy::const_unique_ptr differential_darcy_operator_ = nullptr;
                     
+                    //! Differential Darcy operator. // \todo #820 Rename; it is the second one used in dH.
+                    GlobalVariationalOperatorNS::Darcy::const_unique_ptr delta_darcy_operator_ = nullptr;
+                    
                     //! Inlet pressure (acts as a source).
                     GlobalVariationalOperatorNS::TransientSource<>::const_unique_ptr inlet_pressure_source_operator_ = nullptr;
                     
@@ -931,6 +946,9 @@ namespace HappyHeart
                     //! Parameter that encapsulates velocity part of the current solution.
                     ParameterAtDof<ParameterNS::Type::vector>::type::const_unique_ptr velocity_solution_param_ = nullptr;
                     
+                    //! Parameter that encapsulates delta velocity part of the current solution #820 Improve comment
+                    ParameterAtDof<ParameterNS::Type::vector>::type::const_unique_ptr delta_velocity_solution_param_ = nullptr;
+                    
                     //! Parameter that encapsulates pressure part of the current solution.
                     ParameterAtDof<ParameterNS::Type::scalar>::type::const_unique_ptr pressure_solution_param_ = nullptr;
                     
@@ -954,6 +972,9 @@ namespace HappyHeart
                     
                     //! deltaMixtVariables (\todo #820 Rename and comment!)
                     GlobalVector::unique_ptr delta_mixt_variables_ = nullptr;
+                    
+                    //! deltaDarcyVector (\todo #820 Rename and comment!)
+                    GlobalVector::unique_ptr delta_darcy_vector_ = nullptr;
                 };
                 
                 

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

@@ -296,10 +296,33 @@ namespace HappyHeart
                         
                         
                     }
-
-
                     
                     
+                    auto& delta_darcy_vector = GetNonCstDeltaDarcyVector();
+                    delta_darcy_vector.ZeroEntries(__FILE__, __LINE__);
+                    
+                    {
+                        GlobalVectorWithCoefficient with_coeff(delta_darcy_vector, -1.);
+                        GetDeltaDarcyOperator().Assemble(std::make_tuple(std::ref(with_coeff)),
+                                                         ImplicitStepFluidNS::IsFullDarcy::no);
+                    }
+                    
+                    {
+                        const auto filename = parent::GetOutputDirectory(numbering_subset)
+                        + "/"
+                        + DifferentialPreffix(differential::yes)
+                        + "deltaDarcyVector_" + GetIterationTag() + ".m";
+                        
+                        delta_darcy_vector.View(parent::MpiHappyHeart(),
+                                                  filename,
+                                                  __FILE__, __LINE__,
+                                                  PETSC_VIEWER_ASCII_MATLAB);
+                        
+                        
+                        
+                    }
+                    
+
                     // \todo #820 At the moment, I accept my norm is not the same as Bruno's; see what Petsc gives when
                     // SNES is applied (system_rhs's is what we got in hyperelastic model).
                     const double norm = parent::GetSystemRhs(numbering_subset).Norm(NORM_2, __FILE__, __LINE__);

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

@@ -55,6 +55,16 @@ namespace HappyHeart
                 }
                 
                 
+                template<class HyperelasticLawT>
+                inline const GlobalVariationalOperatorNS::Darcy&
+                VariationalFormulation<HyperelasticLawT>
+                ::GetDeltaDarcyOperator() const noexcept
+                {
+                    assert(!(!delta_darcy_operator_));
+                    return *delta_darcy_operator_;
+                }
+                
+                
                 template<class HyperelasticLawT>
                 inline const Parameter<ParameterNS::Type::vector, ParameterNS::TimeDependencyNS::None>&
                 VariationalFormulation<HyperelasticLawT>::GetInletPressure() const noexcept
@@ -485,6 +495,17 @@ namespace HappyHeart
                 }
                 
                 
+                
+                template<class HyperelasticLawT>
+                inline const ParameterAtDof<ParameterNS::Type::vector>::type& VariationalFormulation<HyperelasticLawT>
+                ::GetDeltaVelocitySolutionParam() const noexcept
+                {
+                    assert(!(!delta_velocity_solution_param_));
+                    return *delta_velocity_solution_param_;
+                }
+
+                
+                
                 template<class HyperelasticLawT>
                 inline const ParameterAtDof<ParameterNS::Type::vector>::type& VariationalFormulation<HyperelasticLawT>
                 ::GetVelSHalfVhfParam() const noexcept
@@ -694,6 +715,21 @@ namespace HappyHeart
                 {
                     return const_cast<GlobalVector&>(GetDeltaMixtVariables());
                 }
+                
+                
+                template<class HyperelasticLawT>
+                inline const GlobalVector& VariationalFormulation<HyperelasticLawT>::GetDeltaDarcyVector() const noexcept
+                {
+                    assert(!(!delta_darcy_vector_));
+                    return *delta_darcy_vector_;
+                }
+                
+                
+                template<class HyperelasticLawT>
+                inline GlobalVector& VariationalFormulation<HyperelasticLawT>::GetNonCstDeltaDarcyVector() noexcept
+                {
+                    return const_cast<GlobalVector&>(GetDeltaDarcyVector());
+                }
 
                 
              

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

@@ -64,12 +64,16 @@ namespace HappyHeart
                     parent::AllocateSystemMatrix(numbering_subset, numbering_subset);
                     parent::AllocateSystemVector(numbering_subset);
                     
-                    work_vector_rhs_like_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
-                    corrected_values_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
-                    solid_differential_velocity_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
-                    non_homogeneous_robin_vector_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
-                    dirichlet_load_state_correction_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
-                    delta_mixt_variables_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
+                    {
+                        decltype(auto) system_solution = parent::GetSystemSolution(numbering_subset);
+                        work_vector_rhs_like_ = std::make_unique<GlobalVector>(system_solution);
+                        corrected_values_ = std::make_unique<GlobalVector>(system_solution);
+                        solid_differential_velocity_ = std::make_unique<GlobalVector>(system_solution);
+                        non_homogeneous_robin_vector_ = std::make_unique<GlobalVector>(system_solution);
+                        dirichlet_load_state_correction_ = std::make_unique<GlobalVector>(system_solution);
+                        delta_mixt_variables_ = std::make_unique<GlobalVector>(system_solution);
+                        delta_darcy_vector_ = std::make_unique<GlobalVector>(system_solution);
+                    }
                     
                     solid_differential_velocity_on_robin_interface_ =
                         std::make_unique<GlobalVector>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)));
@@ -243,6 +247,12 @@ namespace HappyHeart
                                                                                         fluid_velocity,
                                                                                         GetCorrectedValues());
                         
+                        delta_velocity_solution_param_ = std::make_unique<vector_at_dof_type>("Fluid delta velocity solution",
+                                                                                              mesh,
+                                                                                              monolithic_felt_space,
+                                                                                              fluid_velocity,
+                                                                                              GetDeltaMixtVariables());
+                        
                         velSHalfVhf_param_ = std::make_unique<vector_at_dof_type>("velSHalfVhf",
                                                                                   mesh,
                                                                                   monolithic_felt_space,
@@ -269,6 +279,19 @@ namespace HappyHeart
                                                      parent::GetTimeManager(),
                                                      DoComputeProcessorWiseLocal2Global::yes); // \todo #820 Check if no is ok.
                     
+                    delta_darcy_operator_ =
+                        std::make_unique<GVO::Darcy>(monolithic_felt_space,
+                                                     fluid_velocity,
+                                                     mesh_dimension,
+                                                     GetFluidDensity(),
+                                                     porosity,
+                                                     GetDeltaVelocitySolutionParam(),
+                                                     GetPressureSolutionParam(),
+                                                     GetInternalFriction(),
+                                                     parent::GetTimeManager(),
+                                                     DoComputeProcessorWiseLocal2Global::yes); // \todo #820 Check if no is ok.
+
+                    
                     
                     differential_darcy_operator_ =
                         std::make_unique<GVO::DifferentialDarcy>(monolithic_felt_space,

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/LocalVariationalOperatorInstances/DifferentialDarcy.hpp

@@ -41,6 +41,7 @@ namespace HappyHeart
              *
              * \todo Improve the comment by writing its mathematical definition!
              * \todo #559 Time dependency has been momentarily dropped.
+             *
              */
             class DifferentialDarcy final
             : public LinearLocalVariationalOperator<LocalVector>,