#820 Second solve of dH implemented but still to debug (doesn't yield the right result).

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

@@ -15,7 +15,7 @@
 # include <vector>
 
 # include "Utilities/InputParameterList/OpsFunction.hpp"
-
+# include "ThirdParty/Wrappers/Petsc/Solver/Snes.hpp"
 # include "Geometry/Domain/Domain.hpp"
 
 # include "Operators/GlobalVariationalOperatorInstances/BilinearForm/Mass.hpp"
@@ -136,6 +136,7 @@ namespace HappyHeart
                                                     const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
                                                     const Parameter<ParameterNS::Type::vector, ParameterNS::TimeDependencyNS::None>& inlet_pressure,
                                                     const GlobalVector& solid_displacement,
+                                                    const GlobalMatrix& tangent_solid_varf,
                                                     DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
                                                     VariableHolder& variable_holder);
                     
@@ -592,6 +593,9 @@ namespace HappyHeart
                     //! Non constant accessor to the deltaDarcyVector (	\todo #820 Rename and comment!).
                     GlobalVector& GetNonCstMonolithicDeltaDarcyVector() noexcept;
 
+                    //! Accessor to the tangent of solid varf.
+                    const GlobalMatrix& GetTangentSolidVarf() const noexcept;
+
                     
                     ///@}
                     
@@ -647,12 +651,18 @@ namespace HappyHeart
                     GlobalVector& GetNonCstSolidDifferentialVelocityOnRobinInterface() noexcept;
 
                     
-                    //! Constant accessor to the deltaResidual (	odo #820 Rename and comment!)
+                    //! Constant accessor to the deltaResidual (	\todo #820 Rename and comment!)
                     const GlobalVector& GetDeltaResidual() const noexcept;
                     
-                    //! Non constant accessor to the deltaResidual (	odo #820 Rename and comment!)
+                    //! Non constant accessor to the deltaResidual (	\todo #820 Rename and comment!)
                     GlobalVector& GetNonCstDeltaResidual() noexcept;
 
+                    //! Constant accessor to the deltaSolution (	\todo #820 Rename and comment!)
+                    const GlobalVector& GetDeltaSolution() const noexcept;
+                    
+                    //! Non constant accessor to the deltaSolution (	\todo #820 Rename and comment!)
+                    GlobalVector& GetNonCstDeltaSolution() noexcept;
+                    
                     
                 private:
                     
@@ -976,6 +986,8 @@ namespace HappyHeart
                     //! deltaResidual (\todo #820 Rename and comment!)
                     GlobalVector::unique_ptr delta_residual_ = nullptr;
                     
+                    //! deltaSolution (\todo #820 Rename and comment!)   dispDeltaNew in Freefem
+                    GlobalVector::unique_ptr delta_solution_ = nullptr;
                     
                     ///@}
                     
@@ -1016,6 +1028,12 @@ namespace HappyHeart
                
                     //! Value of the fluid mass at the end of previous time step.
                     const GlobalVector& previous_time_step_fluid_mass_on_solid_mesh_;
+                    
+                    //! Tangent of the solid varf. \todo #820 This is rather inelegant; model should be reordered heavily...
+                    const GlobalMatrix& tangent_solid_varf_;
+                    
+                    //! Solver used to solve solid varf.
+                    Wrappers::Petsc::Snes::unique_ptr delta_solid_varf_solver_ = nullptr;
                
                 };
                 

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

@@ -408,6 +408,36 @@ namespace HappyHeart
                                                __FILE__, __LINE__,
                                                PETSC_VIEWER_ASCII_MATLAB);
                         
+                        Wrappers::Petsc::AXPY(-1.,
+                                              delta_solid_vector,
+                                              delta_residual,
+                                              __FILE__, __LINE__);
+                        
+                        delta_residual.Scale(-1., __FILE__, __LINE__);
+                    }
+                    
+                    decltype(auto) tangent_varf = GetTangentSolidVarf();
+                    
+                    delta_solid_varf_solver_->SolveLinear(tangent_varf,
+                                                          tangent_varf,
+                                                          GetDeltaResidual(),
+                                                          GetNonCstDeltaSolution(),
+                                                          __FILE__, __LINE__);
+                    
+                    
+                    
+                    {
+                        decltype(auto) delta_sol = GetDeltaSolution();
+                        const auto filename = parent::GetOutputDirectory(numbering_subset)
+                        + "/"
+                        + DifferentialPreffix(differential::yes)
+                        + "solid_varf_sol_" + GetIterationTag() + ".m";
+                        
+                        delta_sol.View(parent::MpiHappyHeart(),
+                                               filename,
+                                               __FILE__, __LINE__,
+                                               PETSC_VIEWER_ASCII_MATLAB);
+                        
                     }
 
                     

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

@@ -786,6 +786,21 @@ namespace HappyHeart
                 {
                     return const_cast<GlobalVector&>(GetDeltaResidual());
                 }
+                
+                
+                template<class HyperelasticLawT>
+                inline const GlobalVector& VariationalFormulation<HyperelasticLawT>::GetDeltaSolution() const noexcept
+                {
+                    assert(!(!delta_solution_));
+                    return *delta_solution_;
+                }
+                
+                
+                template<class HyperelasticLawT>
+                inline GlobalVector& VariationalFormulation<HyperelasticLawT>::GetNonCstDeltaSolution() noexcept
+                {
+                    return const_cast<GlobalVector&>(GetDeltaSolution());
+                }
 
 
                 template<class HyperelasticLawT>
@@ -796,6 +811,14 @@ namespace HappyHeart
                     return *delta_fluid_mass_param_;
                 }
 
+                
+                template<class HyperelasticLawT>
+                inline const GlobalMatrix& VariationalFormulation<HyperelasticLawT>
+                ::GetTangentSolidVarf() const noexcept
+                {
+                    return tangent_solid_varf_;
+                }
+                
              
             } // namespace InnerLoopNS
                 

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

@@ -38,6 +38,7 @@ namespace HappyHeart
                                          const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
                                          const Parameter<ParameterNS::Type::vector, ParameterNS::TimeDependencyNS::None>& inlet_pressure,
                                          const GlobalVector& solid_displacement,
+                                         const GlobalMatrix& tangent_solid_varf,
                                          DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
                                          VariableHolder& variable_holder)
                 : parent(mpi,
@@ -49,7 +50,8 @@ namespace HappyHeart
                 hyperelastic_law_(hyperelastic_law),
                 solid_displacement_(solid_displacement),
                 inlet_pressure_(inlet_pressure),
-                previous_time_step_fluid_mass_on_solid_mesh_(variable_holder.GetFluidMassVector())
+                previous_time_step_fluid_mass_on_solid_mesh_(variable_holder.GetFluidMassVector()),
+                tangent_solid_varf_(tangent_solid_varf)
                 { }
                 
                 
@@ -187,6 +189,9 @@ namespace HappyHeart
                         delta_residual_ =
                             std::make_unique<GlobalVector>(solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_displacement)));
                         AllocateGlobalVector(solid_god_of_dof, GetNonCstDeltaResidual());
+                        
+                        delta_solution_ = std::make_unique<GlobalVector>(GetDeltaResidual());
+                        
                     }
                     
                     {
@@ -248,6 +253,9 @@ namespace HappyHeart
                     
                     decltype(auto) fluid_pressure =
                         unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_pressure));
+                    
+                    delta_solid_varf_solver_ = InitSolver<EnumUnderlyingType(SolverIndex::main_solver)>(parent::MpiHappyHeart(),
+                                                                                                        input_parameter_data);
 
                     
                     {

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

@@ -94,6 +94,7 @@ namespace HappyHeart
                                                 const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
                                                 const Parameter<ParameterNS::Type::vector, ParameterNS::TimeDependencyNS::None>& inlet_pressure,
                                                 const GlobalVector& solid_displacement,
+                                                const GlobalMatrix& tangent_solid_varf,
                                                 DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
                                                 VariableHolder& variable_holder);
         

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

@@ -34,6 +34,7 @@ namespace HappyHeart
                                      const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
                                      const Parameter<ParameterNS::Type::vector, ParameterNS::TimeDependencyNS::None>& inlet_pressure,
                                      const GlobalVector& solid_displacement,
+                                     const GlobalMatrix& tangent_solid_varf,
                                      DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
                                      VariableHolder& variable_holder)
             : parent(mpi,
@@ -52,6 +53,7 @@ namespace HappyHeart
                                                            porosity,
                                                            inlet_pressure,
                                                            solid_displacement,
+                                                           tangent_solid_varf,
                                                            std::move(copy_bc_list),
                                                            variable_holder);
             }

Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx

@@ -265,12 +265,6 @@ namespace HappyHeart
             
             
             {
-//                auto&& bc_list =
-//                {
-//                    bc_manager.GetDirichletBoundaryConditionPtr(EnumUnderlyingType(BoundaryConditionIndex::fluid_non_homegeneous))
-//                };
-                
-                DirichletBoundaryCondition::vector_shared_ptr bc_list;
                 
                 hyperelastic_law_ =
                     std::make_unique<hyperelastic_law_type>(GetVariableHolder().GetSolid(),
@@ -279,21 +273,10 @@ namespace HappyHeart
                                                             GetInitialPorosityValue(),
                                                             GetFluidMassParameter(),
                                                             GetBulkSolid());
-
-                implicit_step_fluid_varf_ =
-                    std::make_unique<implicit_step_fluid_varf_type>(mpi,
-                                                                    time_manager,
-                                                                    fluid_god_of_dof,
-                                                                    GetHyperelasticLaw(),
-                                                                    porosity_parameter.GetOnFluidMesh(),
-                                                                    GetInletPressure(),
-                                                                    GetSolidOnFluidMesh().GetCurrentDisplacement(),
-                                                                    std::move(bc_list),
-                                                                    GetNonCstVariableHolder());
-                
-                implicit_step_fluid_varf_->Init(input_parameter_data);
             }
             
+
+            
             
             {
                 const auto& solid_god_of_dof = this->GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
@@ -323,7 +306,36 @@ namespace HappyHeart
                                                                            GetSolidVelocityVector(),
                                                                            GetInitialPorosityValue(),
                                                                            GetNonCstVariableHolder());
-            }           
+            }
+            
+            
+            
+            {
+                //                auto&& bc_list =
+                //                {
+                //                    bc_manager.GetDirichletBoundaryConditionPtr(EnumUnderlyingType(BoundaryConditionIndex::fluid_non_homegeneous))
+                //                };
+                
+                DirichletBoundaryCondition::vector_shared_ptr bc_list;
+              
+                const auto& solid_god_of_dof = this->GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
+
+                decltype(auto) solid_ns = solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_displacement));
+                
+                implicit_step_fluid_varf_ =
+                std::make_unique<implicit_step_fluid_varf_type>(mpi,
+                                                                time_manager,
+                                                                fluid_god_of_dof,
+                                                                GetHyperelasticLaw(),
+                                                                porosity_parameter.GetOnFluidMesh(),
+                                                                GetInletPressure(),
+                                                                GetSolidOnFluidMesh().GetCurrentDisplacement(),
+                                                                SolidVariationalFormulationPolicyT::GetVariationalFormulation().GetSystemMatrix(solid_ns, solid_ns),
+                                                                std::move(bc_list),
+                                                                GetNonCstVariableHolder());
+                
+                implicit_step_fluid_varf_->Init(input_parameter_data);
+            }
         }
         
         

Sources/ThirdParty/Wrappers/Petsc/Solver/Snes.hpp

@@ -311,6 +311,8 @@ namespace HappyHeart
     
 } // namespace HappyHeart
 
+
 # include "ThirdParty/Wrappers/Petsc/Solver/Snes.hxx"
 
+
 #endif // HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_SOLVER_x_SNES_HPP_