diff --git a/dumux/geomechanics/poroelastic/couplingmanager.hh b/dumux/geomechanics/poroelastic/couplingmanager.hh
index 4d1da76a8d8ada5536a5e0515f03b413b114c394..5939471352cfbb5cd729485d3811e091189babf3 100644
--- a/dumux/geomechanics/poroelastic/couplingmanager.hh
+++ b/dumux/geomechanics/poroelastic/couplingmanager.hh
@@ -16,6 +16,7 @@
 #include <algorithm>
 #include <type_traits>
 
+#include <dune/common/std/type_traits.hh>
 #include <dumux/discretization/method.hh>
 #include <dumux/common/properties.hh>
 #include <dumux/multidomain/couplingmanager.hh>
@@ -94,6 +95,15 @@ class PoroMechanicsCouplingManager : public virtual CouplingManager< MDTraits >
         std::unique_ptr< ElementVolumeVariables<PMFlowId> > pmFlowElemVolVars;
     };
 
+    template<typename T>
+    using Storage = decltype(std::declval<T>().localResidual().evalStorage(
+        std::declval<T>().fvGeometry(),
+        std::declval<T>().curElemVolVars()
+    ));
+
+    template<typename LA>
+    static constexpr bool hasExperimentalEvalStorage = Dune::Std::is_detected_v<Storage, LA>;
+
 public:
 
     // export the domain ids
@@ -305,10 +315,16 @@ public:
 
         // If the residual instationary, evaluate storage
         if (!pmFlowLocalAssembler.localResidual().isStationary())
-            res += pmFlowLocalAssembler.localResidual().evalStorage(pmFlowLocalAssembler.element(),
-                                                                    pmFlowLocalAssembler.fvGeometry(),
-                                                                    pmFlowLocalAssembler.prevElemVolVars(),
-                                                                    pmFlowLocalAssembler.curElemVolVars());
+        {
+            if constexpr (hasExperimentalEvalStorage<PMFlowLocalAssembler>)
+                res += pmFlowLocalAssembler.localResidual().evalStorage(pmFlowLocalAssembler.fvGeometry(),
+                                                                        pmFlowLocalAssembler.curElemVolVars());
+            else
+                res += pmFlowLocalAssembler.localResidual().evalStorage(pmFlowLocalAssembler.element(),
+                                                                        pmFlowLocalAssembler.fvGeometry(),
+                                                                        pmFlowLocalAssembler.prevElemVolVars(),
+                                                                        pmFlowLocalAssembler.curElemVolVars());
+        }
 
         return res;
     }
diff --git a/test/multidomain/poromechanics/el2p/main.cc b/test/multidomain/poromechanics/el2p/main.cc
index 45b46852335925f59f38bebc01e9ed4026d1da3e..49451ff301275492ce964ce7da1f957e2e2a6239 100644
--- a/test/multidomain/poromechanics/el2p/main.cc
+++ b/test/multidomain/poromechanics/el2p/main.cc
@@ -25,13 +25,18 @@
 #include <dumux/linear/linearsolvertraits.hh>
 #include <dumux/linear/linearalgebratraits.hh>
 
+
 #include <dumux/multidomain/newtonsolver.hh>
-#include <dumux/multidomain/fvassembler.hh>
 #include <dumux/multidomain/traits.hh>
 
 #include <dumux/io/vtkoutputmodule.hh>
 #include <dumux/io/grid/gridmanager_yasp.hh>
 
+#include <dumux/experimental/assembly/multistagemultidomainfvassembler.hh>
+#include <dumux/experimental/timestepping/timelevel.hh>
+#include <dumux/experimental/timestepping/multistagemethods.hh>
+#include <dumux/experimental/timestepping/multistagetimestepper.hh>
+
 #include "properties.hh"
 
 int main(int argc, char** argv)
@@ -138,12 +143,27 @@ int main(int argc, char** argv)
     auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDT);
 
+    const auto timeSteppingMethod = []() -> std::shared_ptr<Experimental::MultiStageMethod<double>>
+    {
+        const auto timeStepScheme = getParam<std::string>("TimeLoop.Scheme", "DIRK3");
+        if (timeStepScheme == "DIRK3")
+            return std::make_shared<Experimental::MultiStage::DIRKThirdOrderAlexander<double>>();
+        else if (timeStepScheme == "ImplicitEuler")
+            return std::make_shared<Experimental::MultiStage::ImplicitEuler<double>>();
+        else if (timeStepScheme == "CrankNicolson")
+            return std::make_shared<Experimental::MultiStage::Theta<double>>(0.5);
+        else
+            DUNE_THROW(Dumux::ParameterException,
+                "Unknown time stepping scheme. Possible values are "
+                << "DIRK3 or ImplicitEuler or CrankNicolson");
+    }();
+
     // the assembler
-    using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric, /*implicit?*/true>;
-    auto assembler = std::make_shared<Assembler>( std::make_tuple(twoPProblem, poroMechProblem),
-                                                  std::make_tuple(twoPFvGridGeometry, poroMechFvGridGeometry),
-                                                  std::make_tuple(twoPGridVariables, poroMechGridVariables),
-                                                  couplingManager, timeLoop, xOld);
+    using Assembler = Experimental::MultiStageMultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(twoPProblem, poroMechProblem),
+                                                 std::make_tuple(twoPFvGridGeometry, poroMechFvGridGeometry),
+                                                 std::make_tuple(twoPGridVariables, poroMechGridVariables),
+                                                 couplingManager, timeSteppingMethod, xOld);
 
     // the linear solver
     using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
@@ -151,13 +171,16 @@ int main(int argc, char** argv)
 
     // the non-linear solver
     using NewtonSolver = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
-    NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+    auto nonLinearSolver = std::make_shared<NewtonSolver>(assembler, linearSolver, couplingManager);
+
+    using TimeStepper = Experimental::MultiStageTimeStepper<NewtonSolver>;
+    TimeStepper timeStepper(nonLinearSolver, timeSteppingMethod);
 
     // time loop
     timeLoop->start(); do
     {
         // solve the non-linear system with time step control
-        nonLinearSolver.solve(x, *timeLoop);
+        timeStepper.step(x, timeLoop->time(), timeLoop->timeStepSize());
 
         // make the new solution the old solution
         xOld = x;
@@ -179,8 +202,9 @@ int main(int argc, char** argv)
         const auto& twoPLocalResidual = assembler->localResidual(twoPId);
         for (const auto& element : elements(leafGridView, Dune::Partitions::interior))
         {
-            auto storageVec = twoPLocalResidual.evalStorage(*twoPProblem, element, *twoPFvGridGeometry, *twoPGridVariables, x[twoPId]);
-            storage += storageVec[0];
+            const auto fvGeometry = localView(*twoPFvGridGeometry).bindElement(element);
+            const auto elemVolVars = localView(twoPGridVariables->curGridVolVars()).bindElement(element, fvGeometry, x[twoPId]);
+            storage += twoPLocalResidual.evalStorage(fvGeometry, elemVolVars)[0];
         }
         std::cout << "time, mass CO2 (kg), mass brine (kg):" << std::endl;
         std::cout << timeLoop->time() << " , " << storage[1] << " , " << storage[0] << std::endl;
@@ -190,7 +214,7 @@ int main(int argc, char** argv)
 
 
     // output some Newton statistics
-    nonLinearSolver.report();
+    nonLinearSolver->report();
 
     timeLoop->finalize(leafGridView.comm());
 
diff --git a/test/multidomain/poromechanics/el2p/params.input b/test/multidomain/poromechanics/el2p/params.input
index bd98c090e5041f1c36e2eac95d3af73d3e035c72..28cce7fc70f80fc17836316eb88fb36df9f833b1 100644
--- a/test/multidomain/poromechanics/el2p/params.input
+++ b/test/multidomain/poromechanics/el2p/params.input
@@ -1,6 +1,7 @@
 [TimeLoop]
 Dt = 10 # [s]
 TEnd = 100 # [s]
+Scheme = ImplicitEuler
 
 [Grid]
 LowerLeft = 0 0 0
diff --git a/test/multidomain/poromechanics/el2p/problem_2p.hh b/test/multidomain/poromechanics/el2p/problem_2p.hh
index efeffbfb1fc7b84f0a6a0553beeb4d63007c6f4a..981645d1096d61547a993fe2925caf6188771243 100644
--- a/test/multidomain/poromechanics/el2p/problem_2p.hh
+++ b/test/multidomain/poromechanics/el2p/problem_2p.hh
@@ -117,9 +117,13 @@ public:
         return values;
     }
 
+    void setTime(Scalar t) const
+    { time_ = t; }
+
 private:
     static constexpr Scalar eps_ = 1.0e-6;
     std::string problemName_;
+    mutable Scalar time_ = 0;
 };
 
 } // end namespace Dumux
diff --git a/test/multidomain/poromechanics/el2p/problem_poroelastic.hh b/test/multidomain/poromechanics/el2p/problem_poroelastic.hh
index e95ced64e36e9570246c5523c3599f55a77d2610..0b8ac1189b5f8ee03cdfeab75f697a236300ea5e 100644
--- a/test/multidomain/poromechanics/el2p/problem_poroelastic.hh
+++ b/test/multidomain/poromechanics/el2p/problem_poroelastic.hh
@@ -102,9 +102,13 @@ public:
                             const SubControlVolume& scv) const
     { return PrimaryVariables(0.0); }
 
+    void setTime(Scalar t) const
+    { time_ = t; }
+
 private:
     static constexpr Scalar eps_ = 3e-6;
     std::string problemName_;
+    mutable Scalar time_ = 0;
 };
 
 } // end namespace Dumux
diff --git a/test/references/test_md_poromechanics_el2p_gravity_2p-reference.vtu b/test/references/test_md_poromechanics_el2p_gravity_2p-reference.vtu
index d92817936501ca1175b9406faf27093869b4a6b7..ae43c34fa0f63372cef53b2e77edd6e1d10a2a3a 100644
--- a/test/references/test_md_poromechanics_el2p_gravity_2p-reference.vtu
+++ b/test/references/test_md_poromechanics_el2p_gravity_2p-reference.vtu
@@ -12,20 +12,20 @@
           1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="p_liq" NumberOfComponents="1" format="ascii">
-          1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07
-          1.50003e+07 1.50003e+07 1.50003e+07 1.50003e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.50016e+07 1.50016e+07 1.5e+07
-          1.5e+07 1.50016e+07 1.50016e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07
-          1.5e+07 1.50016e+07 1.50016e+07 1.5e+07 1.5e+07 1.50016e+07 1.50016e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07
-          1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07 1.5e+07
-          1.5e+07 1.5e+07 1.5e+07 1.5e+07
+          1.33158e+07 1.22817e+07 1.22817e+07 1.33158e+07 1.22817e+07 1.08534e+07 1.08534e+07 1.22817e+07 1.22817e+07 1.08534e+07 1.08534e+07 1.22817e+07
+          1.33158e+07 1.22817e+07 1.22817e+07 1.33158e+07 1.42344e+07 1.39939e+07 1.39939e+07 1.42344e+07 1.39939e+07 1.40106e+07 1.40106e+07 1.39939e+07
+          1.39939e+07 1.40106e+07 1.40106e+07 1.39939e+07 1.42344e+07 1.39939e+07 1.39939e+07 1.42344e+07 1.57732e+07 1.6016e+07 1.6016e+07 1.57732e+07
+          1.6016e+07 1.59998e+07 1.59998e+07 1.6016e+07 1.6016e+07 1.59998e+07 1.59998e+07 1.6016e+07 1.57732e+07 1.6016e+07 1.6016e+07 1.57732e+07
+          1.66941e+07 1.77385e+07 1.77385e+07 1.66941e+07 1.77385e+07 1.91842e+07 1.91842e+07 1.77385e+07 1.77385e+07 1.91842e+07 1.91842e+07 1.77385e+07
+          1.66941e+07 1.77385e+07 1.77385e+07 1.66941e+07
         </DataArray>
         <DataArray type="Float32" Name="rho_liq" NumberOfComponents="1" format="ascii">
-          1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64
-          1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64
-          1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64
-          1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64
-          1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64 1077.64
-          1077.64 1077.64 1077.64 1077.64
+          1076.87 1076.39 1076.39 1076.87 1076.39 1075.72 1075.72 1076.39 1076.39 1075.72 1075.72 1076.39
+          1076.87 1076.39 1076.39 1076.87 1077.29 1077.18 1077.18 1077.29 1077.18 1077.19 1077.19 1077.18
+          1077.18 1077.19 1077.19 1077.18 1077.29 1077.18 1077.18 1077.29 1078 1078.11 1078.11 1078
+          1078.11 1078.1 1078.1 1078.11 1078.11 1078.1 1078.1 1078.11 1078 1078.11 1078.11 1078
+          1078.42 1078.9 1078.9 1078.42 1078.9 1079.56 1079.56 1078.9 1078.9 1079.56 1079.56 1078.9
+          1078.42 1078.9 1078.9 1078.42
         </DataArray>
         <DataArray type="Float32" Name="mob_liq" NumberOfComponents="1" format="ascii">
           692.437 692.437 692.437 692.437 692.437 692.437 692.437 692.437 692.437 692.437 692.437 692.437
@@ -37,27 +37,27 @@
         </DataArray>
         <DataArray type="Float32" Name="S_gas" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 6.7843e-07 6.7843e-07 0
-          0 6.7843e-07 6.7843e-07 0 0 0 0 0 0 0 0 0
-          0 6.79485e-07 6.79485e-07 0 0 6.79485e-07 6.79485e-07 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 6.83979e-07 6.83979e-07 0
+          0 6.83979e-07 6.83979e-07 0 0 0 0 0 0 0 0 0
+          0 6.7429e-07 6.7429e-07 0 0 6.7429e-07 6.7429e-07 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0
         </DataArray>
         <DataArray type="Float32" Name="p_gas" NumberOfComponents="1" format="ascii">
-          1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07
-          1.50179e+07 1.50179e+07 1.50179e+07 1.50179e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50192e+07 1.50192e+07 1.50176e+07
-          1.50176e+07 1.50192e+07 1.50192e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07
-          1.50176e+07 1.50192e+07 1.50192e+07 1.50176e+07 1.50176e+07 1.50192e+07 1.50192e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07
-          1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07
-          1.50176e+07 1.50176e+07 1.50176e+07 1.50176e+07
+          1.33334e+07 1.22993e+07 1.22993e+07 1.33334e+07 1.22993e+07 1.08709e+07 1.08709e+07 1.22993e+07 1.22993e+07 1.08709e+07 1.08709e+07 1.22993e+07
+          1.33334e+07 1.22993e+07 1.22993e+07 1.33334e+07 1.4252e+07 1.40115e+07 1.40115e+07 1.4252e+07 1.40115e+07 1.40282e+07 1.40282e+07 1.40115e+07
+          1.40115e+07 1.40282e+07 1.40282e+07 1.40115e+07 1.4252e+07 1.40115e+07 1.40115e+07 1.4252e+07 1.57908e+07 1.60336e+07 1.60336e+07 1.57908e+07
+          1.60336e+07 1.60174e+07 1.60174e+07 1.60336e+07 1.60336e+07 1.60174e+07 1.60174e+07 1.60336e+07 1.57908e+07 1.60336e+07 1.60336e+07 1.57908e+07
+          1.67117e+07 1.7756e+07 1.7756e+07 1.67117e+07 1.7756e+07 1.92017e+07 1.92017e+07 1.7756e+07 1.7756e+07 1.92017e+07 1.92017e+07 1.7756e+07
+          1.67117e+07 1.7756e+07 1.7756e+07 1.67117e+07
         </DataArray>
         <DataArray type="Float32" Name="rho_gas" NumberOfComponents="1" format="ascii">
-          920.521 920.521 920.521 920.521 920.521 920.521 920.521 920.521 920.521 920.521 920.521 920.521
-          920.521 920.521 920.521 920.521 920.518 920.518 920.518 920.518 920.518 920.529 920.529 920.518
-          920.518 920.529 920.529 920.518 920.518 920.518 920.518 920.518 920.518 920.518 920.518 920.518
-          920.518 920.529 920.529 920.518 920.518 920.529 920.529 920.518 920.518 920.518 920.518 920.518
-          920.518 920.518 920.518 920.518 920.518 920.518 920.518 920.518 920.518 920.518 920.518 920.518
-          920.518 920.518 920.518 920.518
+          908.08 899.669 899.669 908.08 899.669 886.802 886.802 899.669 899.669 886.802 886.802 899.669
+          908.08 899.669 899.669 908.08 915.039 913.259 913.259 915.039 913.259 913.384 913.384 913.259
+          913.259 913.384 913.384 913.259 915.039 913.259 913.259 915.039 925.794 927.403 927.403 925.794
+          927.403 927.296 927.296 927.403 927.403 927.296 927.296 927.403 925.794 927.403 927.403 925.794
+          931.778 938.211 938.211 931.778 938.211 946.58 946.58 938.211 938.211 946.58 946.58 938.211
+          931.778 938.211 938.211 931.778
         </DataArray>
         <DataArray type="Float32" Name="mob_gas" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -76,12 +76,12 @@
           17569.9 17569.9 17569.9 17569.9
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
-          0.20019 0.200355 0.200355 0.20019 0.200355 0.200655 0.200655 0.200355 0.200355 0.200655 0.200655 0.200355
-          0.20019 0.200355 0.200355 0.20019 0.200073 0.200114 0.200114 0.200073 0.200114 0.200156 0.200156 0.200114
-          0.200114 0.200156 0.200156 0.200114 0.200073 0.200114 0.200114 0.200073 0.199927 0.199886 0.199886 0.199927
-          0.199886 0.199845 0.199845 0.199886 0.199886 0.199845 0.199845 0.199886 0.199927 0.199886 0.199886 0.199927
-          0.19981 0.199645 0.199645 0.19981 0.199645 0.199343 0.199343 0.199645 0.199645 0.199343 0.199343 0.199645
-          0.19981 0.199645 0.199645 0.19981
+          0.200144 0.200233 0.200233 0.200144 0.200233 0.200357 0.200357 0.200233 0.200233 0.200357 0.200357 0.200233
+          0.200144 0.200233 0.200233 0.200144 0.200066 0.200086 0.200086 0.200066 0.200086 0.200085 0.200085 0.200086
+          0.200086 0.200085 0.200085 0.200086 0.200066 0.200086 0.200086 0.200066 0.199934 0.199913 0.199913 0.199934
+          0.199913 0.199915 0.199915 0.199913 0.199913 0.199915 0.199915 0.199913 0.199934 0.199913 0.199913 0.199934
+          0.199855 0.199767 0.199767 0.199855 0.199767 0.199644 0.199644 0.199767 0.199767 0.199644 0.199644 0.199767
+          0.199855 0.199767 0.199767 0.199855
         </DataArray>
         <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
diff --git a/test/references/test_md_poromechanics_el2p_gravity_poroelastic-reference.vtu b/test/references/test_md_poromechanics_el2p_gravity_poroelastic-reference.vtu
index fb72dcd8d38b3835d99b47602179b9425ea59165..c7af1bfd842f020237f8e0dbc7adeb62cf65b32c 100644
--- a/test/references/test_md_poromechanics_el2p_gravity_poroelastic-reference.vtu
+++ b/test/references/test_md_poromechanics_el2p_gravity_poroelastic-reference.vtu
@@ -4,43 +4,43 @@
     <Piece NumberOfCells="64" NumberOfPoints="125">
       <PointData Scalars="porosity" Vectors="u">
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
-          0.200047 0.200175 0.200175 0.20064 0.200035 0.200083 0.200083 0.200091 0.200187 0.200688 0.200093 0.200119
-          0.200134 0.200492 0.200079 0.200146 0.200047 0.200175 0.200035 0.200083 0.200187 0.200688 0.200093 0.200119
-          0.200745 0.200152 0.200526 0.200172 0.200187 0.200093 0.200134 0.200492 0.200079 0.200146 0.200526 0.200172
-          0.200378 0.200165 0.200134 0.200079 0.200047 0.200175 0.200035 0.200083 0.200187 0.200093 0.200134 0.200079
-          0.200047 0.200035 0.199982 0.199943 0.199943 0.199842 0.199935 0.199811 0.199953 0.199863 0.199982 0.199943
-          0.199935 0.199811 0.199776 0.19984 0.199935 0.199953 0.199863 0.19984 0.199884 0.199953 0.199982 0.199943
-          0.199935 0.199953 0.199982 0.199936 0.199799 0.199799 0.199426 0.199784 0.199381 0.199833 0.199498 0.199936
-          0.199799 0.199784 0.199381 0.199326 0.199461 0.199784 0.199833 0.199498 0.199461 0.199573 0.199833 0.199936
-          0.199799 0.199784 0.199833 0.199936 0.199953 0.199825 0.199825 0.199358 0.199812 0.199311 0.199866 0.199507
-          0.199953 0.199825 0.199812 0.199311 0.199254 0.199473 0.199812 0.199866 0.199507 0.199473 0.199622 0.199866
-          0.199953 0.199825 0.199812 0.199866 0.199953
+          0.200036 0.20012 0.20012 0.200385 0.200037 0.200075 0.200075 0.200009 0.200123 0.200397 0.200084 0.200042
+          0.200098 0.20032 0.20008 0.200102 0.200036 0.20012 0.200037 0.200075 0.200123 0.200397 0.200084 0.200042
+          0.200414 0.200081 0.200328 0.200126 0.200123 0.200084 0.200098 0.20032 0.20008 0.200102 0.200328 0.200126
+          0.200264 0.200148 0.200098 0.20008 0.200036 0.20012 0.200037 0.200075 0.200123 0.200084 0.200098 0.20008
+          0.200036 0.200037 0.199983 0.199956 0.199956 0.199909 0.199947 0.199875 0.199961 0.199908 0.199983 0.199956
+          0.199947 0.199875 0.199835 0.199884 0.199947 0.199961 0.199908 0.199884 0.199913 0.199961 0.199983 0.199956
+          0.199947 0.199961 0.199983 0.199943 0.199849 0.199849 0.199699 0.199845 0.199687 0.199861 0.19967 0.199943
+          0.199849 0.199845 0.199687 0.199671 0.199661 0.199845 0.199861 0.19967 0.199661 0.199675 0.199861 0.199943
+          0.199849 0.199845 0.199861 0.199943 0.199964 0.19988 0.19988 0.199616 0.199877 0.199604 0.199902 0.19968
+          0.199964 0.19988 0.199877 0.199604 0.199587 0.199672 0.199877 0.199902 0.19968 0.199672 0.199736 0.199902
+          0.199964 0.19988 0.199877 0.199902 0.199964
         </DataArray>
         <DataArray type="Float32" Name="u" NumberOfComponents="3" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0.0203967 0.0203967 0.196367
-          0 0 0 0 0 0 0 0 0 -5.14796e-12 0.0232633 0.223806
-          0 0 0 0 0 0 0 0 0 -0.0203967 0.0203967 0.196367
+          0 0 0 0 0 0 0 0 0 0.0255126 0.0255126 0.129579
+          0 0 0 0 0 0 0 0 0 1.89588e-13 0.0262156 0.136124
+          0 0 0 0 0 0 0 0 0 -0.0255126 0.0255126 0.129579
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0.0232633 -1.12797e-13 0.223806
-          0 0 0 -6.92946e-12 5.21607e-13 0.263098 0 0 0 -0.0232633 5.73196e-13 0.223806
+          0 0 0 0 0 0 0 0 0 0.0262156 -4.03665e-13 0.136124
+          0 0 0 3.62161e-14 -3.33315e-13 0.147887 0 0 0 -0.0262156 -9.36753e-14 0.136124
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0.0203967 -0.0203967 0.196367 0 0 0 -4.19502e-12 -0.0232633 0.223806
-          0 0 0 -0.0203967 -0.0203967 0.196367 0 0 0 0 0 0
+          0 0 0 0.0255126 -0.0255126 0.129579 0 0 0 1.95653e-13 -0.0262156 0.136124
+          0 0 0 -0.0255126 -0.0255126 0.129579 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 -9.5166e-06 -9.5166e-06 0.247008 0 0 0 -8.22174e-12 -1.94599e-05 0.292708
-          0 0 0 9.51661e-06 -9.5166e-06 0.247008 0 0 0 0 0 0
-          0 0 0 -1.94599e-05 1.06135e-12 0.292708 -8.84654e-12 4.37292e-12 0.356343 1.94599e-05 3.83762e-12 0.292708
-          0 0 0 0 0 0 -9.5166e-06 9.5166e-06 0.247008 -5.53219e-12 1.94599e-05 0.292708
-          9.5166e-06 9.5166e-06 0.247008 0 0 0 0 0 0 0 0 0
+          0 0 0 -3.51365e-05 -3.51365e-05 0.160592 0 0 0 3.63365e-14 -3.26044e-05 0.186667
+          0 0 0 3.51365e-05 -3.51365e-05 0.160592 0 0 0 0 0 0
+          0 0 0 -3.26044e-05 1.88149e-13 0.186667 5.19604e-13 1.10238e-12 0.225174 3.26044e-05 5.37361e-13 0.186667
+          0 0 0 0 0 0 -3.51365e-05 3.51365e-05 0.160592 4.46337e-13 3.26044e-05 0.186667
+          3.51365e-05 3.51365e-05 0.160592 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 -0.0204115 -0.0204115 0.196421 0 0 0
-          -5.43843e-12 -0.0232888 0.223877 0 0 0 0.0204115 -0.0204115 0.196421 0 0 0
-          0 0 0 0 0 0 -0.0232888 1.35531e-12 0.223877 -5.67541e-12 4.58532e-12 0.263199
-          0.0232888 3.16176e-12 0.223877 0 0 0 0 0 0 -0.0204115 0.0204115 0.196421
-          -4.04598e-12 0.0232888 0.223877 0.0204115 0.0204115 0.196421 0 0 0 0 0 0
+          0 0 0 0 0 0 -0.0256535 -0.0256535 0.129392 0 0 0
+          6.0197e-14 -0.0263612 0.135815 0 0 0 0.0256535 -0.0256535 0.129392 0 0 0
+          0 0 0 0 0 0 -0.0263612 1.89377e-13 0.135815 2.22009e-13 9.20984e-13 0.147416
+          0.0263612 6.53446e-13 0.135815 0 0 0 0 0 0 -0.0256535 0.0256535 0.129392
+          5.18473e-13 0.0263612 0.135815 0.0256535 0.0256535 0.129392 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
diff --git a/test/references/test_md_poromechanics_el2p_poroelastic-reference.vtu b/test/references/test_md_poromechanics_el2p_poroelastic-reference.vtu
index a8b7bb53a5a8967cc29bbb3f4af55d2b90fdc4ea..f2cce2c00048bc1675de1abf3e7e3d8d0e9326b1 100644
--- a/test/references/test_md_poromechanics_el2p_poroelastic-reference.vtu
+++ b/test/references/test_md_poromechanics_el2p_poroelastic-reference.vtu
@@ -18,29 +18,29 @@
         </DataArray>
         <DataArray type="Float32" Name="u" NumberOfComponents="3" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 -7.41448e-06 -7.41448e-06 -7.41448e-06
-          0 0 0 0 0 0 0 0 0 -2.95919e-15 -1.27515e-05 -1.27515e-05
-          0 0 0 0 0 0 0 0 0 7.41448e-06 -7.41448e-06 -7.41448e-06
+          0 0 0 0 0 0 0 0 0 -3.2895e-06 -3.2895e-06 -3.28931e-06
+          0 0 0 0 0 0 0 0 0 2.02808e-14 -5.721e-06 -5.72087e-06
+          0 0 0 0 0 0 0 0 0 3.2895e-06 -3.2895e-06 -3.28931e-06
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 -1.27515e-05 5.54018e-14 -1.27515e-05
-          0 0 0 -2.30069e-15 6.79714e-14 -2.31085e-05 0 0 0 1.27515e-05 2.71535e-14 -1.27515e-05
+          0 0 0 0 0 0 0 0 0 -5.721e-06 -2.18471e-14 -5.72087e-06
+          0 0 0 1.85528e-14 -3.5086e-14 -1.13883e-05 0 0 0 5.721e-06 -2.23616e-14 -5.72087e-06
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 -7.41448e-06 7.41448e-06 -7.41448e-06 0 0 0 -2.33167e-14 1.27515e-05 -1.27515e-05
-          0 0 0 7.41448e-06 7.41448e-06 -7.41448e-06 0 0 0 0 0 0
+          0 0 0 -3.2895e-06 3.2895e-06 -3.28931e-06 0 0 0 -7.92271e-15 5.721e-06 -5.72087e-06
+          0 0 0 3.2895e-06 3.2895e-06 -3.28931e-06 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 -1.27515e-05 -1.27515e-05 4.04532e-14 0 0 0 1.69178e-15 -2.31085e-05 4.08255e-14
-          0 0 0 1.27515e-05 -1.27515e-05 2.5728e-14 0 0 0 0 0 0
-          0 0 0 -2.31085e-05 6.13031e-14 4.08955e-14 3.0484e-14 6.97951e-14 3.50766e-14 2.31085e-05 4.31379e-14 3.74653e-14
-          0 0 0 0 0 0 -1.27515e-05 1.27515e-05 1.21806e-14 -3.49545e-15 2.31085e-05 2.80989e-14
-          1.27515e-05 1.27515e-05 4.02227e-14 0 0 0 0 0 0 0 0 0
+          0 0 0 -5.72109e-06 -5.72109e-06 9.35078e-11 0 0 0 -5.42103e-15 -1.13885e-05 1.2128e-10
+          0 0 0 5.72109e-06 -5.72109e-06 9.35185e-11 0 0 0 0 0 0
+          0 0 0 -1.13885e-05 -1.05518e-14 1.21289e-10 1.54563e-14 -1.09011e-14 1.66959e-10 1.13885e-05 -6.02914e-15 1.21296e-10
+          0 0 0 0 0 0 -5.72109e-06 5.72109e-06 9.35062e-11 9.24641e-15 1.13885e-05 1.21295e-10
+          5.72109e-06 5.72109e-06 9.35101e-11 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 -7.41448e-06 -7.41448e-06 7.41448e-06 0 0 0
-          1.5632e-14 -1.27515e-05 1.27515e-05 0 0 0 7.41448e-06 -7.41448e-06 7.41448e-06 0 0 0
-          0 0 0 0 0 0 -1.27515e-05 1.21722e-14 1.27515e-05 3.82396e-14 1.4547e-14 2.31085e-05
-          1.27515e-05 2.73105e-14 1.27515e-05 0 0 0 0 0 0 -7.41448e-06 7.41448e-06 7.41448e-06
-          2.06957e-14 1.27515e-05 1.27515e-05 7.41448e-06 7.41448e-06 7.41448e-06 0 0 0 0 0 0
+          0 0 0 0 0 0 -3.28954e-06 -3.28954e-06 3.28957e-06 0 0 0
+          1.65277e-14 -5.72109e-06 5.72113e-06 0 0 0 3.28954e-06 -3.28954e-06 3.28957e-06 0 0 0
+          0 0 0 0 0 0 -5.72109e-06 -1.11086e-16 5.72113e-06 1.17338e-14 -9.91828e-15 1.13885e-05
+          5.72109e-06 4.13952e-15 5.72113e-06 0 0 0 0 0 0 -3.28954e-06 3.28954e-06 3.28957e-06
+          2.71674e-14 5.72109e-06 5.72113e-06 3.28954e-06 3.28954e-06 3.28957e-06 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0