#820 Start writing dH and remove some debug lines.

Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hpp

@@ -179,11 +179,13 @@ namespace HappyHeart
 
         private:
             
-            //! Enum class to tag each GlobalVector related tgo the solid.
+            //! Enum class to tag each GlobalVector related to the solid.
             enum class SolidIndex : std::size_t
             {
                 displacement_k0 = 0,
                 H0,
+                delta_displacement,
+                dH_input_vector,
                 Nvector
             };
 
@@ -494,6 +496,9 @@ namespace HappyHeart
              */
             std::array<GlobalVector::unique_ptr, EnumUnderlyingType(SolidIndex::Nvector)> solid_vector_list_
             = Utilities::NullptrArray<GlobalVector::unique_ptr, EnumUnderlyingType(SolidIndex::Nvector)>();
+                
+            //! Index of gmres iteration. \todo #820 For dev purpose!
+            unsigned int gmres_iteration_index_ = 0u;
 
 
             

Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx

@@ -87,6 +87,18 @@ namespace HappyHeart
                     
                 }
                 
+                std::cout << "Gmres call" << std::endl;
+                auto& delta_displacement = GetNonCstVector<SolidIndex::delta_displacement>();
+                
+                delta_displacement.ZeroEntries(__FILE__, __LINE__);
+                
+                gmres_iteration_index_ = static_cast<unsigned int>(-1);
+                dH_solver_->SolveLinear(gmres_shell_matrix, gmres_shell_matrix,
+                                        GetVector<SolidIndex::H0>(),
+                                        delta_displacement,
+                                        __FILE__, __LINE__);
+
+                
                 
                 
                 // \todo #820 Adapt following block to poro case (it was directly lifted from FSI!);
@@ -228,7 +240,21 @@ namespace HappyHeart
         void Model<SolidVariationalFormulationPolicyT>
         ::ShellMatrixOperation(Vec petsc_input_vector, Vec& petsc_output_vector)
         {
-            // \todo #820 Implement it!
+            ++gmres_iteration_index_;
+            
+            auto& input_vector = GetNonCstVector<SolidIndex::dH_input_vector>();
+            input_vector.SetFromPetscVec(petsc_input_vector, __FILE__, __LINE__);
+
+            input_vector.Scale(parent::GetTimeManager().GetInverseTimeStep(), __FILE__, __LINE__);
+            
+            
+            // We can't take directly the internal pointer of a wrapped vector which pattern would be similar to input_displ_vector
+            // : residual computation is wrong if we do.
+//            auto error_code = VecCopy(*HappyHeart vector*, petsc_output_vector);
+//            
+//            if (error_code)
+//                throw Wrappers::Petsc::ExceptionNS::Exception(error_code, "VecCopy",
+//                                                              __FILE__, __LINE__);
         }
 
         

Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx

@@ -97,7 +97,8 @@ namespace HappyHeart
                 
                 GetNonCstVectorPtr<SolidIndex::displacement_k0>() = std::make_unique<GlobalVector>(solution);
                 GetNonCstVectorPtr<SolidIndex::H0>() = std::make_unique<GlobalVector>(solution);
-                
+                GetNonCstVectorPtr<SolidIndex::delta_displacement>() = std::make_unique<GlobalVector>(solution);
+                GetNonCstVectorPtr<SolidIndex::dH_input_vector>() = std::make_unique<GlobalVector>(solution);
             }
             
             

Sources/ModelInstances/UnderDevelopment/Poromechanics/SolidVariationalFormulationPolicy/Hyperelasticity/VariationalFormulation.hxx

@@ -228,8 +228,6 @@ namespace HappyHeart
             inline void VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>
             ::ComputeResidual()
             {
-                std::cout << "========= COMPUTE RESIDUAL" << std::endl;
-                
                 const auto& numbering_subset = GetNumberingSubset();
                 
                 auto& vectors_and_matrices = this->GetNonCstVectorsAndMatrices();
@@ -277,8 +275,6 @@ namespace HappyHeart
                 + std::to_string(parent::GetTimeManager().NtimeModified()) + "_it_1_"
                 + "it_newton_" + std::to_string(parent::GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".m";
                 
-                std::cout << "WRITE " << filename << std::endl;
-                
                 rhs.View(parent::MpiHappyHeart(),
                          filename,
                          __FILE__, __LINE__,
@@ -294,8 +290,6 @@ namespace HappyHeart
             inline void VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>
             ::ComputeTangent()
             {
-                std::cout << "========= COMPUTE TANGENT" << std::endl;
-                
                 const auto& numbering_subset = GetNumberingSubset();
 
                 const auto& vectors_and_matrices = this->GetVectorsAndMatrices();
@@ -320,8 +314,6 @@ namespace HappyHeart
                 + std::to_string(parent::GetTimeManager().NtimeModified()) + "_it_1_"
                 + "it_newton_" + std::to_string(parent::GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".m";
                 
-                std::cout << "WRITE " << filename << std::endl;
-                
                 system_matrix.View(parent::MpiHappyHeart(),
                                    filename,
                                    __FILE__, __LINE__,
@@ -338,8 +330,6 @@ namespace HappyHeart
             void VariationalFormulationHyperElasticity<HyperelasticLawT, TimeSchemeT>
             ::UpdateVectorsAndMatrices(const Vec& a_evaluation_state)
             {
-                std::cout << "============ UPDATE VECTOR AND MATRIX " << std::endl;
-                
                 auto& vm = this->GetNonCstVectorsAndMatrices();
                 
                 vm.UpdateEvaluationState(a_evaluation_state);
@@ -362,9 +352,6 @@ namespace HappyHeart
                                                                   false,
                                                                   vm);
                 }
-                
-                vm.GetVectorNewStiffness().View(parent::MpiHappyHeart(), __FILE__, __LINE__);
-                
             }