diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index fb2e8c701669a07c8c6320ae9f07d3d3be511efa..35acb2579cd8accc343403988b1fb2996dc4cb6a 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -79,6 +79,8 @@ private:
     const QRuleEdge& rule_Kplus1_Edge =
       QRulesEdge::rule(edge, problem_.quadratureOrder_K() + 1);
     const QRuleEdge& rule_2ndOrder = QRulesEdge::rule(edge, 2);
+    const QRuleEdge& rule_boundary =
+      QRulesEdge::rule(edge, problem_.quadratureOrderBoundary());
 
     //we use the basis
     // phi_elem,0 (x) = indicator(elem);
@@ -95,7 +97,7 @@ private:
       //basis function phi_elem,i
       const int elemIdxSLE = (dim + 1)*elemIdx;
 
-      //fill load vector b using a quadrature rule
+      //fill load vector b using a quadrature rule (bulk term)
       for (const auto& ip : rule_q)
       {
         //quadrature point in local and global coordinates
@@ -104,16 +106,16 @@ private:
 
         const Scalar weight(ip.weight());
         const Scalar integrationElem(geo.integrationElement(qpLocal));
-        const Scalar qEvalution(problem_.q(qpGlobal));
+        const Scalar qEvaluation(problem_.q(qpGlobal));
 
         //quadrature for int_elem q*phi_elem,0 dV
-        b[elemIdxSLE] += weight * integrationElem * qEvalution;
+        b[elemIdxSLE] += weight * integrationElem * qEvaluation;
 
         //quadrature for int_elem q*phi_elem,i dV
         for (int i = 0; i < dim; i++)
         {
           b[elemIdxSLE + i + 1] +=
-            weight * integrationElem * qpGlobal[i] * qEvalution;
+            weight * integrationElem * qpGlobal[i] * qEvaluation;
         }
       }
 
@@ -206,13 +208,13 @@ private:
 
           const Scalar weight(ip.weight());
           const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
-          const Dune::FieldMatrix<Scalar, dim, dim> KEvalution =
+          const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
             problem_.K(qpGlobal);
 
           for (int i = 0; i < dim; i++)
           {
             KIntegrals[i] +=
-              weight * integrationElem * (KEvalution[i] * normal);
+              weight * integrationElem * (KEvaluation[i] * normal);
           }
         }
 
@@ -225,7 +227,7 @@ private:
 
           const Scalar weight(ip.weight());
           const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
-          const Dune::FieldMatrix<Scalar, dim, dim> KEvalution =
+          const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
             problem_.K(qpGlobal);
 
           for (int i = 0; i < dim; i++)
@@ -233,7 +235,7 @@ private:
             for (int j = 0; j < dim; j++)
             {
               KIntegralsX[i][j] += weight * integrationElem *
-                qpGlobal[j] * (KEvalution[i] * normal);
+                qpGlobal[j] * (KEvaluation[i] * normal);
             }
           }
         }
@@ -354,24 +356,52 @@ private:
         }
         else //boundary facet
         {
+          //fill load vector b using a quadrature rule (boundary terms)
+          for (const auto& ip : rule_boundary)
+          {
+            //quadrature point in local and global coordinates
+            const auto& qpLocal = ip.position();
+            const auto& qpGlobal = intersctGeo.global(qpLocal);
+
+            const Scalar weight(ip.weight());
+            const Scalar
+              integrationElem(intersctGeo.integrationElement(qpLocal));
+            const Scalar boundaryEvaluation(problem_.boundary(qpGlobal));
+            const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
+              problem_.K(qpGlobal);
+
+            //quadrature for int_intersct mu * g * phi_elem,0 ds
+            b[elemIdxSLE] += weight * integrationElem * mu * boundaryEvaluation;
+
+            //quadrature for
+            // int_intersct (mu * phi_elem,i - K[:,i] * normal) * g ds
+            //with K*grad(phi_elem,i) = K[:,i]
+            for (int i = 0; i < dim; i++)
+            { 
+              b[elemIdxSLE + i + 1] +=
+                weight * integrationElem * (mu * qpGlobal[i] -
+                  (KEvaluation[i] * normal)) * boundaryEvaluation;
+            }
+          }
+
           for (int i = 0; i < dim; i++)
-          { //we use the relations
-            // int_intersct mu * jump(phi_elem,0) * jump(phi_elem,i) ds
-            // = mu * int_intersct x[i] ds
-            //and for boundary facets
-            // int_intersct avg(K*grad(phi_elem,i)) * jump(phi_elem,0) ds
-            // = 0.5 * int_intersct K[:,i] * normal ds
+          {
+            //we evaluate the integrals
+            // int_intersct mu * phi_elem,0 * phi_elem,i ds
+            //and
+            // int_intersct phi_elem,0 * K*grad(phi_elem,i) * normal ds
+            //with K*grad(phi_elem,i) = K[:,i]
             A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
               mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
 
             for (int j = 0; j <= i; j++)
             {
-              //we use the relations
-              // int_intersct mu * jump(phi_elem,i) * jump(phi_elem,j) ds
-              // = mu * int_intersct x[i] * x[j] ds
-              //and for boundary facets
-              // int_intersct avg(K*grad(phi_elem,i)) * jump(phi_elem,j) ds
-              // = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
+              //we evaluate the integrals
+                // int_intersct mu * phi_elem,i * phi_elem,j ds
+                //and
+                // int_intersct (phi_elem,i * K[:,j]
+                //    + phi_elem,j * K[:,i]) * normal ds
+                //with K*grad(phi_elem,i) = K[:,i]
               A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
                 mu * quadraticIntregrals[i][j]
                 - 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh
index 942c4338083017f5eebbd2c91adcc00e5749db3b..64cc67fd5b9c163db702e1e184b8ccdbf0d13138 100644
--- a/dune/mmdg/problems/dgproblem.hh
+++ b/dune/mmdg/problems/dgproblem.hh
@@ -57,6 +57,13 @@ class DGProblem
       return 1;
     }
 
+    //returns the recommended quadrature order to compute an integral
+    //over x * boundary(x)
+    virtual int quadratureOrderBoundary () const
+    {
+      return 1;
+    }
+
     //returns the recommended quadrature order to compute an integral
     //over an entry K[i][j] of the permeability tensor
     virtual int quadratureOrder_K () const