diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp
index 194e7f092e4fa2116f9e04ab1b34e0a6e9d95f96..b09492b913f117d37537f4b806b0792fdc130020 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp
@@ -323,12 +323,7 @@ namespace HappyHeart
                     //! Non constant accessor to the \#820 [DEV] Vector defined on fluid mass numbering subset.
                     GlobalVector& GetNonCstWorkFluidMass() noexcept;
                     
-                    //! Constant accessor to the \#820 [DEV] Vector defined on fluid velocity numbering subset.
-                    const GlobalVector& GetWorkFluidVelocity() const noexcept;
-                    
-                    //! Non constant accessor to the \#820 [DEV] Vector defined on fluid velocity numbering subset.
-                    GlobalVector& GetNonCstWorkFluidVelocity() noexcept;
-
+                
                     //! Constant accessor to the mass div Velocity operator.
                     const GlobalVariationalOperatorNS::ScalarDivVectorial& GetMassDivVelocity() const noexcept;
                     
@@ -732,9 +727,6 @@ namespace HappyHeart
                     //! \#820 [DEV] Vector defined on fluid mass numbering subset.
                     GlobalVector::unique_ptr work_fluid_mass_ = nullptr;
                     
-                    //! \#820 [DEV] Vector defined on fluid velocity numbering subset.
-                    GlobalVector::unique_ptr work_fluid_velocity_ = nullptr;
-                    
                     //! \#820 [DEV] Vector defined on the numbering subset for solid on fluid.
                     GlobalVector::unique_ptr work_solid_velocity_ = nullptr;
 
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
index 0a6db232182fad14dd35c5e77b70ec99ca83e263..1f08bb47dc89c34d6145bf8251377a89fec13016 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
@@ -66,6 +66,10 @@ namespace HappyHeart
                                                  corrected_values,
                                                  __FILE__, __LINE__);
                         
+                        auto& vel_contrib_on_vel_ns = variable_holder.GetVelSHalfVhfZeroedOutOfRobinInterface();
+                        
+                        
+                        
                         variable_holder.template DevPrint<DevPhase::newton_fp>(variable_holder.GetVelSHalfVhf(),
                                                                                "dirichlet_just_before_newton_fp_loop");
                         
@@ -589,9 +593,10 @@ namespace HappyHeart
                     GlobalVectorWithCoefficient contribution_with_coeff(contribution, 1.);
                     
                     decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
+                    decltype(auto) variable_holder = this->GetVariableHolder();
                     
                     {
-                        auto& vel_only = GetNonCstWorkFluidVelocity();
+                        auto& vel_only = variable_holder.GetNonCstWorkFluidVelocity();
                         
                         Wrappers::Petsc::MatMult(interpolator_holder.GetMonolithic2Velocity().GetInterpolationMatrix(),
                                                  GetDarcyOperator().GetVelocitySolutionParam().GetGlobalVector(),
@@ -647,18 +652,18 @@ namespace HappyHeart
                 void VariationalFormulation<HyperelasticLawT>::SetSolidVelocityForDifferential(Vec petsc_solid_velocity)
                 {
                     auto& delta_displacement = GetNonCstDeltaDisplacement();
-                    auto& work_vel_fluid = GetNonCstWorkFluidVelocity();
+                    
+                    decltype(auto) variable_holder = this->GetNonCstVariableHolder();
+                    
+                    auto& work_vel_fluid = variable_holder.GetNonCstWorkFluidVelocity();
 
                     delta_displacement.SetFromPetscVec(petsc_solid_velocity, __FILE__, __LINE__);
                     
-                    decltype(auto) variable_holder = this->GetVariableHolder();
-                    
                     {
                         variable_holder.template DevPrint<DevPhase::dH>(delta_displacement,
                                                                         "initial_disp_delta");
 
                     }
-
                     
                     Wrappers::Petsc::MatMult(variable_holder_parent::GetVariableHolder().GetSolidOnFluidMesh().GetMatrixSolidToFluid(),
                                              delta_displacement,
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationAccessors.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationAccessors.hxx
index 5e08715953810a0fbed15f48a8a0f986cbc34072..f241ee2b2e86cab227fbe6745ba105236e2602d9 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationAccessors.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationAccessors.hxx
@@ -120,22 +120,7 @@ namespace HappyHeart
                 {
                     return const_cast<GlobalVector&>(GetWorkFluidMass());
                 }
-                
-                
-                template<class HyperelasticLawT>
-                inline const GlobalVector& VariationalFormulation<HyperelasticLawT>::GetWorkFluidVelocity() const noexcept
-                {
-                    assert(!(!work_fluid_velocity_));
-                    return *work_fluid_velocity_;
-                }
-                
-                
-                template<class HyperelasticLawT>
-                inline GlobalVector& VariationalFormulation<HyperelasticLawT>::GetNonCstWorkFluidVelocity() noexcept
-                {
-                    return const_cast<GlobalVector&>(GetWorkFluidVelocity());
-                }
-                
+          
                 
                 template<class HyperelasticLawT>
                 inline const GlobalVariationalOperatorNS::ScalarDivVectorial& VariationalFormulation<HyperelasticLawT>
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationInit.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationInit.hxx
index d81d28a9dfd624db0d818592053127eb9aef0903..505e8df99dacdbf982d38d028fec6325df20669a 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationInit.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationInit.hxx
@@ -145,10 +145,6 @@ namespace HappyHeart
                     
                     AllocateGlobalVector(god_of_dof, GetNonCstWorkFluidMass());
                     
-                    work_fluid_velocity_ =
-                        std::make_unique<GlobalVector>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
-                    
-                    AllocateGlobalVector(god_of_dof, GetNonCstWorkFluidVelocity());
                     
                     delta_fluid_mass_ = std::make_unique<GlobalVector>(GetWorkFluidMass());
                     
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
index 68e9918290377b6882260307251ff63dd4b4f5a3..0e7d91685ac5089b380390701202dcc609038e2c 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
@@ -78,6 +78,8 @@ namespace HappyHeart
                 double absolute_error = 1.e20;
                 double relative_error = 1.e20;
                 
+                decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
+                
                 {
                     auto& velSHalfVhf = variable_holder.GetNonCstVelSHalfVhf();
                     
@@ -89,15 +91,17 @@ namespace HappyHeart
 
                     variable_holder.GetNonCstSolidOnFluidMesh().ComputeSolidVelocity(velSHalfVhf);
                     
-                    // \todo #820 Improve this!
-                    // We get value on the whole mesh; we want to zero all values not on the Robin interface.
-                    // Current stupid method is to interpolate back-and-forth; it's obviously not very efficient.
-                    
-                    
+                    auto& velSHalfVhf_nullified_out_of_robin_interface
+                        = variable_holder.GetNonCstVelSHalfVhfZeroedOutOfRobinInterface();
                     
+                    // We get value on the whole mesh; we want to zero all values not on the Robin interface.
+                    Wrappers::Petsc::MatMult(interpolator_holder.GetVelocityNullifierOutOfRobinInterface(),
+                                             velSHalfVhf,
+                                             velSHalfVhf_nullified_out_of_robin_interface,
+                                             __FILE__, __LINE__);
                     
-                    variable_holder.template DevPrint<DevPhase::none>(velSHalfVhf,
-                                                                      "velSHalfVhfr_at_begin");
+                    variable_holder.template DevPrint<DevPhase::none>(velSHalfVhf_nullified_out_of_robin_interface,
+                                                                      "velSHalfVhfr_non_robin_zero_at_begin");
                 }
 
                 
@@ -110,9 +114,6 @@ namespace HappyHeart
                     ++newton_fp_index;
                 }
                 
-                decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
-                
-
                 // Update fluidmass, velocity fluid and pressure.
                 Wrappers::Petsc::MatMult(interpolator_holder.GetMonolithic2Mass().GetInterpolationMatrix(),
                                          inner_varf.GetCorrectedValues(),
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.cpp b/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.cpp
index 0e52f314d0f293006cd01fc0838d21d6ac81a1b2..1c7889fcbe8a69220ad051991ab0983c23cfdf64 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.cpp
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.cpp
@@ -134,10 +134,18 @@ namespace HappyHeart
                 std::make_unique<type>(fluid_felt_space,
                                        fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
                                        robin_felt_space,
-                                       fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
+                                       fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)));
                 
                 velocity_to_robin_interface_->Init();
                 
+                velocity_robin_interface_to_full_mesh_ =
+                std::make_unique<type>(robin_felt_space,
+                                       fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)),
+                                       fluid_felt_space,
+                                       fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
+                
+                velocity_robin_interface_to_full_mesh_->Init();
+                
                 monolithic_2_robin_velocity_ =
                 std::make_unique<GlobalMatrix>(fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
                                                fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure)));
@@ -148,6 +156,22 @@ namespace HappyHeart
                                             __FILE__, __LINE__);
             }
             
+            {
+                decltype(auto) numbering_subset =
+                    fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
+                
+                velocity_nullifier_out_of_robin_interface_ = std::make_unique<GlobalMatrix>(numbering_subset,
+                                                                                            numbering_subset);
+                
+                Wrappers::Petsc::MatMatMult(GetVelocityRobinInterfaceToFullMesh().GetInterpolationMatrix(),
+                                            GetVelocityToRobinInterface().GetInterpolationMatrix(),
+                                            *velocity_nullifier_out_of_robin_interface_,
+                                            __FILE__, __LINE__);
+                
+                
+            }
+            
+            
             
         }
 
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hpp b/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hpp
index 10a2dda0cffb1d648e82f9b8955857337228b44e..631dd3b2bd052dd57411afb1c01a1f6472adcc7e 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hpp
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hpp
@@ -113,6 +113,9 @@ namespace HappyHeart
             //! Constant accessor to the matrix to go from (whole mesh, monolithic) to (robin interface, velocity).
             const GlobalMatrix& GetMonolithic2RobinVelocity() const noexcept;
             
+            //! Constant accessor to the matrix used to zero all velocity terms that aren't on the Robin interface.
+            const GlobalMatrix& GetVelocityNullifierOutOfRobinInterface() const noexcept;
+            
         private:
             
             /*!
@@ -121,6 +124,7 @@ namespace HappyHeart
              */
             GlobalMatrix& GetNonCstMonolithic2RobinVelocity() noexcept;
             
+            
         private:
             
             //! Interpolator fluid mass on solid mesh -> fluid mass on fluid mesh.
@@ -173,7 +177,8 @@ namespace HappyHeart
             //! Matrix to go from (whole mesh, monolithic) to (robin interface, velocity).
             GlobalMatrix::unique_ptr monolithic_2_robin_velocity_ = nullptr;
             
-
+            //! Matrix used to zero all velocity terms that aren't on the Robin interface.
+            GlobalMatrix::unique_ptr velocity_nullifier_out_of_robin_interface_ = nullptr;
             ///@}
 
             
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hxx
index 6808c893dcb5b267f862bb391630faaf83bc442d..8c86d62e3ccf1bb5799766c669beb19e418898ae 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hxx
@@ -131,6 +131,12 @@ namespace HappyHeart
         }
         
 
+        inline const GlobalMatrix& InterpolatorHolder
+        ::GetVelocityNullifierOutOfRobinInterface() const noexcept
+        {
+            assert(!(!velocity_nullifier_out_of_robin_interface_));
+            return *velocity_nullifier_out_of_robin_interface_;
+        }
         
         
     } // namespace PoromechanicsNS
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp b/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp
index fd657a5aa761c79b5868b3ebd2cbf318fb230660..34f62f6b8c1a74c50d9f696811c4066865ef265f 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp
@@ -115,6 +115,9 @@ namespace HappyHeart
                     fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
                 velSHalfVhf_ = std::make_unique<GlobalVector>(numbering_subset);
                 AllocateGlobalVector(fluid_god_of_dof, GetNonCstVelSHalfVhf());
+                
+                work_fluid_velocity_ = std::make_unique<GlobalVector>(GetVelSHalfVhf());
+                velSHalfVhf_zeroed_out_of_robin_interface_ = std::make_unique<GlobalVector>(GetVelSHalfVhf());
             }
             
         }
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp b/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp
index 40b5d02264081ba1fc382939713cafc0dab1b9d0..0cd922c6877a6e487456538d9bd6939f0e567ba9 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp
@@ -153,6 +153,15 @@ namespace HappyHeart
             //! Non constant accessor to the 	odo #820 Rename this (Currently same as Bruno's Freefem)
             GlobalVector& GetNonCstVelSHalfVhf() noexcept;
             
+            //! Constant accessor to the 	odo #820 Rename this (Currently same as Bruno's Freefem)
+            const GlobalVector& GetVelSHalfVhfZeroedOutOfRobinInterface() const noexcept;
+            
+            //! Non constant accessor to the 	odo #820 Rename this (Currently same as Bruno's Freefem)
+            GlobalVector& GetNonCstVelSHalfVhfZeroedOutOfRobinInterface() noexcept;
+            
+            //! Non constant accessor to the \#820 [DEV] Vector defined on fluid velocity numbering subset.
+            GlobalVector& GetNonCstWorkFluidVelocity() const noexcept;
+            
             
             //! Debug indexes...
             unsigned int implicit_loop_index = 0u;
@@ -170,6 +179,8 @@ namespace HappyHeart
             // \todo #820 At the moment Bruno's name is kept!
             GlobalVector::unique_ptr velSHalfVhf_ = nullptr;
             
+            GlobalVector::unique_ptr velSHalfVhf_zeroed_out_of_robin_interface_ = nullptr;
+            
             
             ///@}
             
@@ -216,7 +227,12 @@ namespace HappyHeart
             
             //!
             const TimeManager& time_manager_;
-        
+            
+        private:
+            
+            
+            //! \#820 [DEV] Vector defined on fluid velocity numbering subset.
+            mutable GlobalVector::unique_ptr work_fluid_velocity_ = nullptr;
          
         };
         
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx
index eefe0d5744045c484e5cd37f4a7fbd7d3e0e2327..ce00788c53fbf356e5651708a7ef2c1371021bce 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx
@@ -159,7 +159,27 @@ namespace HappyHeart
         {
             return const_cast<GlobalVector&>(GetVelSHalfVhf());
         }
-
+        
+        
+        inline const GlobalVector& VariableHolder::GetVelSHalfVhfZeroedOutOfRobinInterface() const noexcept
+        {
+            assert(!(!velSHalfVhf_zeroed_out_of_robin_interface_));
+            return *velSHalfVhf_zeroed_out_of_robin_interface_;
+        }
+        
+        
+        inline GlobalVector& VariableHolder::GetNonCstVelSHalfVhfZeroedOutOfRobinInterface() noexcept
+        {
+            return const_cast<GlobalVector&>(GetVelSHalfVhfZeroedOutOfRobinInterface());
+        }
+        
+        
+        inline GlobalVector& VariableHolder::GetNonCstWorkFluidVelocity() const noexcept
+        {
+            assert(!(!work_fluid_velocity_));
+            return *work_fluid_velocity_;
+        }
+        
         
     } // namespace PoromechanicsNS