diff --git a/dune/dorie/solver/flux_reconstruction/localengeine.hh b/dune/dorie/solver/flux_reconstruction/localengeine.hh
index c6de3fe8829d5b09586b6f1c0b68629938dacf9d..09cb159c8c20ee6bca112ab5684412d89850575f 100644
--- a/dune/dorie/solver/flux_reconstruction/localengeine.hh
+++ b/dune/dorie/solver/flux_reconstruction/localengeine.hh
@@ -181,37 +181,29 @@ public:
   template<typename EG, typename LFSWC>
   void onBindLFSV(const EG & eg, const LFSWC & lfsw_cache)
   {
+    lfsv_volume.bind( eg.entity() );
+
+    std::cout << "lfsv_volume.size(): " << lfsv_volume.size() << std::endl;
+    std::cout << "lfsv_volume.debug():" << std::endl; lfsv_volume.debug();
+    
     global_pl_view.bind(lfsw_cache);
     pl.resize(lfsw_cache.size());
 
-    lfsv_volume.bind( eg.entity() );
 
     rl.assign(lfsv_volume.size(),0.0);
   }
 
   template<typename IG, typename LFSUC, typename LFSWC>
   void onBindLFSUVInside(const IG & ig, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
-  {
-    // This is never called
-  }
-
-  template<typename IG, typename LFSUC, typename LFSWC>
-  void onBindLFSUVOutside(const IG & ig,
-                          const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
-                          const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
   {
     lfsv_skeleton.bind( ig.intersection() );
-    lfsv_skeleton_n.bind( ig.intersection() );
 
     // inside part
-    std::cout << "lfsu_s_cache.size(): " << lfsu_s_cache.size() << std::endl;
-    std::cout << "lfsu_s.debug():" << std::endl; lfsu_s_cache.localFunctionSpace().debug();
-
-    std::cout << "lfsw_s_cache.size(): " << lfsw_s_cache.size() << std::endl;
-    std::cout << "lfsw_s.debug():" << std::endl; lfsw_s_cache.localFunctionSpace().debug();
+    std::cout << "lfsu_cache.size(): " << lfsu_cache.size() << std::endl;
+    std::cout << "lfsu.debug():" << std::endl; lfsu_cache.localFunctionSpace().debug();
 
-    std::cout << "lfsv_volume.size(): " << lfsv_volume.size() << std::endl;
-    std::cout << "lfsv_volume.debug():" << std::endl; lfsv_volume.debug();
+    std::cout << "lfsw_cache.size(): " << lfsw_cache.size() << std::endl;
+    std::cout << "lfsw.debug():" << std::endl; lfsw_cache.localFunctionSpace().debug();
 
     std::cout << "lfsv_skeleton.size(): " << lfsv_skeleton.size() << std::endl;
     std::cout << "lfsv_skeleton.debug():" << std::endl; lfsv_skeleton.debug();
@@ -219,16 +211,18 @@ public:
     std::cout << "ig.entity().subEntities(1): " << ig.inside().subEntities(1) << std::endl;
 
     // Check that local matrix will be conforming
-    assert(lfsu_s_cache.size() == lfsw_s_cache.size());
-    assert(lfsu_s_cache.size() == (lfsv_volume.size() + lfsv_skeleton.size() * ig.inside().subEntities(1)));
-
-    global_sl_view.bind(lfsu_s_cache);
-    xl.assign(lfsu_s_cache.size(),0.0);
+    assert(lfsu_cache.size() == lfsw_cache.size());
+  }
 
-    global_pl_view.bind(lfsw_s_cache);
-    pl.resize(lfsw_s_cache.size());
+  template<typename IG, typename LFSUC, typename LFSWC>
+  void onBindLFSUVOutside(const IG & ig,
+                          const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
+                          const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
+  {
+    lfsv_skeleton_n.bind( ig.intersection() ); // TODO switch inside/outside
 
-    // outside part
+    // Check that local matrix will be conforming
+    assert(lfsu_s_cache.size() == (lfsv_volume.size() + lfsv_skeleton.size() * ig.inside().subEntities(1)));
     assert(lfsv_skeleton.size() == lfsv_skeleton_n.size());
 
     global_sn_view.bind(lfsu_n_cache);
@@ -241,7 +235,8 @@ public:
   template<typename IG, typename LFSWC>
   void onBindLFSVInside(const IG & ig, const LFSWC & lfsw_cache)
   {
-    // This is never called
+    // global_rl_view.bind(lfsw_s_cache);
+    rl.assign(lfsw_cache.size(),0.0);
   }
 
   template<typename IG, typename LFSWC>
@@ -249,8 +244,6 @@ public:
                          const LFSWC & lfsw_s_cache,
                          const LFSWC & lfsw_n_cache)
   {
-    // global_rl_view.bind(lfsw_s_cache);
-    rl.assign(lfsw_s_cache.size(),0.0);
     // global_rn_view.bind(lfsw_n_cache);
     rn.assign(lfsw_n_cache.size(),0.0);
   }
@@ -262,11 +255,19 @@ public:
                      const LFSUC & lfsu_cache, 
                      const LFSVC & lfsv_cache)
   {
+    std::cout << "r_vec" << std::endl;
+    std::cout << r_vec << std::endl;
+    std::cout << "m_matrix" << std::endl;
+    std::cout << m_matrix << std::endl;
+
+    x_vec.resize(lfsu_cache.size(),0.0);
+    m_matrix.solve(x_vec,r_vec);
+
     auto& lfsu = lfsu_cache.localFunctionSpace();
 
-    // for (unsigned int i = 0; i < lfsu.size(); ++i)
-    //   xl(lfsu,i) = x_vec[i];
-    // global_sl_view.add(xl);
+    for (unsigned int i = 0; i < lfsu.size(); ++i)
+      xl(lfsu,i) = x_vec[i];
+    global_sl_view.write(xl);
   }
   //! @}
 
@@ -381,8 +382,15 @@ public:
 
       for (unsigned int i = 0; i<lfsv_volume.size(); i++)
         for (unsigned int j = 0; j<lfsu.size(); j++)
-          m_matrix[i][j] = (phiu[j]*phiv[i])/detB * factor/detB;
+          m_matrix[i][j] += (phiu[j]*phiv[i])/detB * factor/detB;
     }
+
+    std::cout << "assembleUVVolume" << std::endl;
+    std::cout << "r_vec" << std::endl;
+    std::cout << r_vec << std::endl;
+    std::cout << "m_matrix" << std::endl;
+    std::cout << m_matrix << std::endl;
+
   }
 
   template<typename EG, typename LFSWC>
@@ -393,6 +401,9 @@ public:
   void assembleUVSkeleton(const IG & ig, const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
                           const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
   {
+
+    std::cout << "assembleUVSkeleton" << std::endl;
+
     rl_view.setWeight(local_assembler.weight());
     rn_view.setWeight(local_assembler.weight());
     Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaSkeleton>::
@@ -404,8 +415,10 @@ public:
     unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
     
     // Load residual view into the local residual vector
-    for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i)
-      r_vec[i+offset] += rl(lfsv_skeleton,i) + rn(lfsv_skeleton_n,i);
+    for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
+      std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << " + " << rn(lfsv_skeleton_n,i) << std::endl;
+      r_vec[i+offset] += rl(lfsv_skeleton,i);// + rn(lfsv_skeleton_n,i);
+    }
 
     auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
 
@@ -451,10 +464,16 @@ public:
 
       for (unsigned int i = 0; i<lfsv_skeleton.size(); i++)
         for (unsigned int j = 0; j<lfsu_s.size(); j++){
-          auto xx = y*phiu_s[j]/detB * phiv_s[i+offset] * factor;
-          m_matrix[i+offset][j] += xx;
+          m_matrix[i+offset][j] += y*phiu_s[j]/detB * phiv_s[i] * factor;
         }
     }
+
+
+    std::cout << "r_vec" << std::endl;
+    std::cout << r_vec << std::endl;
+    std::cout << "m_matrix" << std::endl;
+    std::cout << m_matrix << std::endl;
+
   }
 
   template<typename IG, typename LFSWC>
@@ -468,10 +487,15 @@ public:
     Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
       alpha_boundary(lop,ig,lfsw_s_cache.localFunctionSpace(),pl,lfsv_skeleton,rl_view);
 
+    std::cout << "assembleUVBoundary" << std::endl;
+
+
     unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
     // Load residual view into the local residual vector
-    for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i)
+    for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
+      std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << std::endl;
       r_vec[i+offset] += rl(lfsv_skeleton,i);
+    }
 
     auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
 
@@ -517,8 +541,13 @@ public:
 
       for (unsigned int i = 0; i<lfsv_skeleton.size(); i++)
         for (unsigned int j = 0; j<lfsu_s.size(); j++)
-          m_matrix[i+offset][j] += (y*phiu_s[j])/detB * phiv_s[i+offset] * factor;
+          m_matrix[i+offset][j] += (y*phiu_s[j])/detB * phiv_s[i] * factor;
     }
+
+    std::cout << "r_vec" << std::endl;
+    std::cout << r_vec << std::endl;
+    std::cout << "m_matrix" << std::endl;
+    std::cout << m_matrix << std::endl;
   }
 
   template<typename IG, typename LFSWC>
@@ -549,9 +578,6 @@ public:
     rl_view.setWeight(local_assembler.weight());
     Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolumePostSkeleton>::
       alpha_volume_post_skeleton(lop,eg,lfsw_cache.localFunctionSpace(),pl,lfsv_volume,rl_view);
-
-    x_vec.resize(lfsu_cache.size(),0.0);
-    m_matrix.solve(x_vec,r_vec);
   }
 
   template<typename EG, typename LFSWC>
diff --git a/dune/dorie/solver/operator_DG.hh b/dune/dorie/solver/operator_DG.hh
index e52370df78c49695a9133262afe513422e1a03f1..a6d292dfe2e036dd8a4d81c6105712f15a5e032b 100755
--- a/dune/dorie/solver/operator_DG.hh
+++ b/dune/dorie/solver/operator_DG.hh
@@ -284,6 +284,9 @@ public:
 
         std::vector<Vector> gradphiv;
 
+        // discrete gradient sign.
+        double dg_sign = 1.;
+
         if constexpr (dimRangeTest == dim) {
           // (we assume lfsv = RTk)
           lfsv.finiteElement().localBasis().evaluateFunction(p_local,gradphiv);
@@ -295,6 +298,7 @@ public:
             y /= detB;
             BTransposed.mtv(y,gradphiv[i]);
           }
+          dg_sign = -1.;
         } else {
           // (we assume Galerkin method lfsu = lfsv)
           gradphiv = gradphiu;
@@ -316,7 +320,7 @@ public:
         // update residual
         for (unsigned int i = 0; i<lfsv.size(); i++)
         {
-          r.accumulate(lfsv,i, cond * (gradu*gradphiv[i]) * factor);
+          r.accumulate(lfsv,i, dg_sign * cond * (gradu*gradphiv[i]) * factor);
         }
       }
     }
@@ -401,11 +405,15 @@ public:
         std::vector<Scalar> phiv_s(lfsv_s.size());
         std::vector<Scalar> phiv_n(lfsv_n.size());
 
+        // discrete gradient sign.
+        double dg_sign = 1.;
         // evaluate basis functions of test function
         if constexpr (dimDomainLocalTest == dim-1) {
           // (we assume non-Galerkin method lfsu != lfsv)
           lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
           lfsv_n.finiteElement().localBasis().evaluateFunction(it.position(),phiv_n);
+          // discrete gradient sign.
+          dg_sign = -1.;
         } else {
           // (we assume Galerkin method lfsu = lfsv)
           phiv_s = phiu_s;
@@ -543,23 +551,23 @@ public:
         if (upwinding == RichardsDGUpwinding::none)
         {
           for (unsigned int i = 0; i < lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s, i, theta * omega_s * (tgradphi_s[i] * normal) * relCond_s * satCond_s * jump * factor);
+            r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (tgradphi_s[i] * normal) * relCond_s * satCond_s * jump * factor);
           for (unsigned int i = 0; i < lfsv_n.size(); i++)
-            r_n.accumulate(lfsv_n, i, theta * omega_n * (tgradphi_n[i] * normal) * relCond_n * satCond_n * jump * factor);
+            r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (tgradphi_n[i] * normal) * relCond_n * satCond_n * jump * factor);
         }
         else if (upwinding == RichardsDGUpwinding::semiUpwind)
         {
           for (unsigned int i = 0; i<lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s,i, theta * omega_s * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
+            r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
           for (unsigned int i = 0; i<lfsv_n.size(); i++)
-            r_n.accumulate(lfsv_n,i, theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
+            r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
         }
         else if (upwinding == RichardsDGUpwinding::fullUpwind)
         {
           for (unsigned int i = 0; i < lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s, i, theta * (tgradphi_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
+            r_s.accumulate(lfsv_s, i, dg_sign * theta * (tgradphi_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
           for (unsigned int i = 0; i < lfsv_n.size(); i++)
-            r_n.accumulate(lfsv_n, i, theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
+            r_n.accumulate(lfsv_n, i, dg_sign * theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
         }
       }
     }
@@ -619,10 +627,13 @@ public:
 
         std::vector<Scalar> phiv_s(lfsv_s.size());
 
+        // discrete gradient sign.
+        double dg_sign = 1.;
         // evaluate basis functions of test function
         if constexpr (dimDomainLocalTest == dim-1) {
           // (we assume non-Galerkin method lfsu != lfsv)
           lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
+          dg_sign = -1.;
         } else {
           // (we assume Galerkin method lfsu = lfsv)
           phiv_s = phiu_s;
@@ -819,7 +830,7 @@ public:
           // (non-)symmetric IP term
           // symmetry term
           for (unsigned int i = 0; i<lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s,i, theta * (tgradphi_s[i] * normal) * relCond * satCond_s * jump * factor);
+            r_s.accumulate(lfsv_s,i, dg_sign * theta * (tgradphi_s[i] * normal) * relCond * satCond_s * jump * factor);
 
           continue;
         }