#1022 Toward fixing Poromechanics issue: interpolators for T21 to go to monolithic must be changed.

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

@@ -501,7 +501,10 @@ namespace HappyHeart
 
                     GlobalMatrix& GetNonCstIntermediateT12() noexcept;
 
-                    GlobalMatrix& GetNonCstT12Block() noexcept;
+                    GlobalMatrix& GetNonCstT12Block() noexcept;                    
+                    
+                    GlobalMatrix& GetNonCstIntermediateT21Monolithic() noexcept;
+
 
                     /*!
                      * \brief Work matrix that gets the exact same structure as system matrix.

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/VariationalFormulationAccessors.hxx

@@ -165,6 +165,16 @@ namespace HappyHeart
                     return *intermediate_t12_;
 
                 }
+                
+                
+                template<class HyperelasticLawT>
+                inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
+                ::GetNonCstIntermediateT21Monolithic() noexcept
+                {
+                    assert(!(!intermediate_t21_monolithic_));
+                    return *intermediate_t21_monolithic_;
+                    
+                }
 
 
                 template<class HyperelasticLawT>
@@ -175,7 +185,6 @@ namespace HappyHeart
                 }
 
 
-
                 template<class HyperelasticLawT>
                 inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>::GetNonCstT22Block() noexcept
                 {

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/VariationalFormulationComputeMatrixBlock.hxx

@@ -272,6 +272,9 @@ namespace HappyHeart
                     
                     decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
                     
+                    std::cout << "COMPUTE T21; DO REUSE = "
+                    << (do_reuse_matrix == Wrappers::Petsc::DoReuseMatrix::yes ? "yes" : "no") << std::endl;
+                    
                     switch(do_recompute_block)
                     {
                         case DoRecomputeBlock::no:
@@ -302,7 +305,7 @@ namespace HappyHeart
                             decltype(auto) monolithic_to_mass =
                             interpolator_holder.GetMatrixMonolithicToFluidmassOnFluid();
                             
-                            auto& work_var = *intermediate_t21_monolithic_;
+                            auto& work_var = GetNonCstIntermediateT21Monolithic();
                             
                             Wrappers::Petsc::MatMatMult(t21_on_solid_mesh,
                                                         monolithic_to_mass,

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/VariationalFormulationInit.hxx

@@ -101,6 +101,7 @@ namespace HappyHeart
                     t22_in_monolithic_ = std::make_unique<GlobalMatrix>(parent::GetSystemMatrix(numbering_subset, numbering_subset));
 
                     intermediate_t12_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset); // Numbering subset inexact but we don't care here!
+                    
 
                     t12_block_ =
                         std::make_unique<GlobalMatrix>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass)),
@@ -110,6 +111,8 @@ namespace HappyHeart
                         std::make_unique<GlobalMatrix>(solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_solid)),
                                                        solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_on_solid)));
                     
+                    intermediate_t21_monolithic_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset); // Numbering subset inexact but we don't care here!
+                    
                     t22_block_ =
                         std::make_unique<GlobalMatrix>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
                                                        god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));