diff --git a/dumux/boxmodels/common/boxfvelementgeometry.hh b/dumux/boxmodels/common/boxfvelementgeometry.hh
index 4c4f33f258b1bd3178c525f7a81a86c104b2952f..8ac1925afe60518843ffd45d62626d6062ef301c 100644
--- a/dumux/boxmodels/common/boxfvelementgeometry.hh
+++ b/dumux/boxmodels/common/boxfvelementgeometry.hh
@@ -692,7 +692,7 @@ public:
 
             // calculate the shape function gradients
             //typedef Dune::FieldVector< Dune::FieldVector< CoordScalar, dim >, 1 > ShapeJacobian;
-            typedef Dune::FieldVector< CoordScalar, 1 >                           ShapeValue;
+            typedef Dune::FieldVector< CoordScalar, 1 > ShapeValue;
             std::vector<ShapeJacobian> localJac;
             std::vector<ShapeValue>    shapeVal;
             localFiniteElement.localBasis().evaluateJacobian(subContVolFace[k].ipLocal, localJac);
@@ -816,17 +816,13 @@ public:
                         }
                     }
 
-    //                typedef Dune::FieldVector< CoordScalar, 1 > ShapeValue;
                     std::vector<ShapeJacobian> localJac;
-    //                std::vector<ShapeValue>    shapeVal;
                     localFiniteElement.localBasis().evaluateJacobian(subContVol[vert].localCenter, localJac);
-    //                localFiniteElement.localBasis().evaluateFunction(subContVol[vert].ipLocal, shapeVal);
 
                     Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT =
                             geometry.jacobianInverseTransposed(subContVol[vert].localCenter);
                     for (int vert = 0; vert < numVertices; vert++)
                         jacInvT.mv(localJac[vert][0], subContVol[vert].gradCenter[vert]);
-    //                    subContVol[vert].shapeValue[vert] = Scalar(shapeVal[vert]);
                 }
         }
     }
diff --git a/dumux/common/boundarytypes.hh b/dumux/common/boundarytypes.hh
index 520346b69c8eb2fe1c5cf184b7726170d57b0942..af4994c0f61b61645f774f7726be26d84a9d0e76 100644
--- a/dumux/common/boundarytypes.hh
+++ b/dumux/common/boundarytypes.hh
@@ -51,6 +51,7 @@ public:
 
             boundaryInfo_[i].isDirichlet = 0;
             boundaryInfo_[i].isNeumann = 0;
+            boundaryInfo_[i].isOutflow = 0;
 
             eq2pvIdx_[i] = i;
             pv2eqIdx_[i] = i;
@@ -85,6 +86,7 @@ public:
             boundaryInfo_[eqIdx].visited = 1;
             boundaryInfo_[eqIdx].isDirichlet = 0;
             boundaryInfo_[eqIdx].isNeumann = 1;
+            boundaryInfo_[eqIdx].isOutflow = 0;
 
             Valgrind::SetDefined(boundaryInfo_[eqIdx]);
         }
@@ -99,6 +101,7 @@ public:
             boundaryInfo_[eqIdx].visited = 1;
             boundaryInfo_[eqIdx].isDirichlet = 1;
             boundaryInfo_[eqIdx].isNeumann = 0;
+            boundaryInfo_[eqIdx].isOutflow = 0;
 
             eq2pvIdx_[eqIdx] = eqIdx;
             pv2eqIdx_[eqIdx] = eqIdx;
@@ -107,6 +110,21 @@ public:
         }
     }
 
+    /*!
+     * \brief Set all boundary conditions to neuman.
+     */
+    void setAllOutflow()
+    {
+        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+            boundaryInfo_[eqIdx].visited = 1;
+            boundaryInfo_[eqIdx].isDirichlet = 0;
+            boundaryInfo_[eqIdx].isNeumann = 0;
+            boundaryInfo_[eqIdx].isOutflow = 1;
+
+            Valgrind::SetDefined(boundaryInfo_[eqIdx]);
+        }
+    }
+
     /*!
      * \brief Set a neumann boundary condition for a single a single
      *        equation.
@@ -114,6 +132,9 @@ public:
     void setNeumann(int eqIdx)
     {
         boundaryInfo_[eqIdx].visited = 1;
+        boundaryInfo_[eqIdx].isDirichlet = 0;
+        boundaryInfo_[eqIdx].isNeumann = 1;
+        boundaryInfo_[eqIdx].isOutflow = 0;
 
         Valgrind::SetDefined(boundaryInfo_[eqIdx]);
     }
@@ -133,6 +154,7 @@ public:
         boundaryInfo_[eqIdx].visited = 1;
         boundaryInfo_[eqIdx].isDirichlet = 1;
         boundaryInfo_[eqIdx].isNeumann = 0;
+        boundaryInfo_[eqIdx].isOutflow = 0;
 
         // update the equation <-> primary variable mapping
         eq2pvIdx_[eqIdx] = pvIdx;
@@ -141,6 +163,20 @@ public:
         Valgrind::SetDefined(boundaryInfo_[eqIdx]);
     }
 
+    /*!
+     * \brief Set a neumann boundary condition for a single a single
+     *        equation.
+     */
+    void setOutflow(int eqIdx)
+    {
+        boundaryInfo_[eqIdx].visited = 1;
+        boundaryInfo_[eqIdx].isDirichlet = 0;
+        boundaryInfo_[eqIdx].isNeumann = 0;
+        boundaryInfo_[eqIdx].isOutflow = 1;
+
+        Valgrind::SetDefined(boundaryInfo_[eqIdx]);
+    }
+
     /*!
      * \brief Set a dirichlet boundary condition for a single primary
      *        variable
@@ -191,6 +227,25 @@ public:
         return false;
     };
 
+    /*!
+     * \brief Returns true if an equation is used to specify an
+     *        outflow condition.
+     */
+    bool isOutflow(unsigned eqIdx) const
+    { return boundaryInfo_[eqIdx].isOutflow; };
+
+    /*!
+     * \brief Returns true if some equation is used to specify an
+     *        outflow condition.
+     */
+    bool hasOutflow() const
+    {
+        for (int i = 0; i < numEq; ++i)
+            if (boundaryInfo_[i].isOutflow)
+                return true;
+        return false;
+    };
+
     /*!
      * \brief Returns the index of the equation which should be used
      *        for the dirichlet condition of the pvIdx's primary
@@ -213,6 +268,9 @@ private:
         unsigned char visited : 1;
         unsigned char isDirichlet : 1;
         unsigned char isNeumann : 1;
+        unsigned char isOutflow : 1;
+//        unsigned char isCouplingInflow : 1;
+//        unsigned char isCouplingOutflow : 1;
     } boundaryInfo_[numEq];
 
     unsigned char eq2pvIdx_[numEq];