diff --git a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc
index b581a125708f44aea0c2a8f03897b9569532f07f..5c6054e103ae21ddf40e40a32c39ad4527ece53e 100644
--- a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc
+++ b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc
@@ -1,4 +1,3 @@
-
 #include "config.h"
 #include "dune/common/parallel/mpihelper.hh"
 #include "dune/pdelab/stationary/linearproblem.hh"
@@ -94,6 +93,25 @@
 #include "sumfact_poisson_dg_3d_facedir_facemod_variation_localoperator_msh35.hh"
 #endif
 
+template <typename T>
+bool isclose(T a, T b, double rel_tol=1e-9, double abs_tol=0.0){
+  // Basically a python isclose (PEP 485) without error checking...
+  return std::abs(a-b) <= std::max(rel_tol*std::max(std::abs(a), std::abs(b)), abs_tol);
+}
+
+template<typename T>
+bool fuzzy_is_permutation(T a, T b, double rel_tol=1e-9, double abs_tol=0.0){
+  assert (a.size()==b.size());
+  std::sort(a.begin(), a.end());
+  std::sort(b.begin(), b.end());
+
+  for (std::size_t i=0; i<a.size(); ++i){
+    if (!isclose(a[i], b[i], rel_tol, abs_tol))
+      return false;
+  }
+  return true;
+}
+
 int main(int argc, char** argv){
   try
   {
@@ -180,7 +198,24 @@ int main(int argc, char** argv){
       Dune::PDELab::interpolate(interpolate,dg1_gfs_,x_r);
       r=0.0;
       go_r.residual(x_r, r);
+      // Dune::printvector(std::cout, native(x_r), "input vector", "row");
+      std::cout.precision(17);
+
       Dune::printvector(std::cout, native(r), "residual vector", "row");
+
+      std::vector<RangeType> solution {0, -28.738355599166038, 0, -39.789562942352646, 0, -28.634891475085336, 0, -39.403880102483868, 28.738355599166038, 0, 39.789562942352646, 0, 28.634891475085336, 0, 39.403880102483868, 0};
+
+      std::vector<RangeType> residual(16);
+      for (std::size_t i=0; i<16; ++i){
+        residual[i] = native(r)[i];
+      }
+
+      std::cout << "palpo is_permuation: "
+                << fuzzy_is_permutation(solution, residual) << std::endl;
+
+      if (!fuzzy_is_permutation(solution, residual))
+        return 1;
+
     }
     if (initree.get<bool>("printmatrix", false)) {
       using Dune::PDELab::Backend::native;