diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
index 3ef04b1eb3a6f50650b74ff9ad1cf3380613f7b9..877bda21be9437ae223500a29d24fa4b25b7aec9 100644
--- a/examples/algebraicSolver.m
+++ b/examples/algebraicSolver.m
@@ -115,7 +115,7 @@ title('error estimator over number of total computational cost')
 %% local function for residual a posteriori error estimation
 % \eta(T)^2 = h_T^2 * || \Delta u ||_{L^2(T)}^2
 %               + h_T * || [[(Du - fvec) * n]] ||_{L^2(E) \cap \Omega}^2
-function indicators = estimate(blf, lf, u)
+function indicators = estimate(~, lf, u)
 p = u.fes.finiteElement.order;
 mesh =  u.fes.mesh;
 trafo = getAffineTransformation(mesh);
diff --git a/examples/convergenceRates.m b/examples/convergenceRates.m
index f6c1c35943b2fc953088281b9ff3b31faee86148..f8b6455f400b32a5dcf74fa926cfa62d4284fed5 100644
--- a/examples/convergenceRates.m
+++ b/examples/convergenceRates.m
@@ -87,9 +87,10 @@ end
 % \eta(T)^2 = h_T^2 * || \Delta u + f ||_{L^2(T)}^2
 %               + h_T * || [[Du * n]] ||_{L^2(E) \cap \Omega}^2
 %               + h_T * || Du * n - phi ||_{L^2(E) \cap \Gamma_N}^2
-function indicators = estimate(blf, lf, u)
-    p = u.fes.finiteElement.order;
-    mesh =  u.fes.mesh;
+function indicators = estimate(~, lf, u)
+    fes = u.fes;
+    p = fes.finiteElement.order;
+    mesh =  fes.mesh;
     
     % compute volume residual element-wise
     % For p=1, the diffusion term vanishes in the residual.
@@ -103,9 +104,11 @@ function indicators = estimate(blf, lf, u)
     
     % compute edge residual edge-wise
     qr = QuadratureRule.ofOrder(p, '1D');
+    dirichletBnd = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
+    neumannBnd = getCombinedBndEdges(mesh, fes.bnd.neumann);
     edgeRes = integrateNormalJump(Gradient(u), qr, ...
-        @(j) zeros(size(j)), {}, mesh.boundaries{1}, ...    % no jump on dirichlet edges
-        @(j,phi) j-phi, {lf.neumann}, mesh.boundaries{2},...% neumann jump on neumann edges
+        @(j) zeros(size(j)), {}, dirichletBnd, ...          % no jump on dirichlet edges
+        @(j,phi) j-phi, {lf.neumann}, neumannBnd,...        % neumann jump on neumann edges
         @(j) j.^2, {}, ':');                                % normal jump on inner edges
     
     % combine the resdiuals suitably
diff --git a/examples/goafemIterativeSolver.m b/examples/goafemIterativeSolver.m
index 50a7df3835d03592ed0f894e33081a62b8935177..357b256be0535e36cdfd97636a33d70296c106e2 100644
--- a/examples/goafemIterativeSolver.m
+++ b/examples/goafemIterativeSolver.m
@@ -116,7 +116,7 @@ title(['GOAFEM p=', num2str(p)])
 %% local function for residual a posteriori error estimation
 % \eta(T)^2 = h_T^2 * || \Delta u ||_{L^2(T)}^2
 %               + h_T * || [[(Du - fvec) * n]] ||_{L^2(E) \cap \Omega}^2
-function indicators = estimate(blf, lf, u)
+function indicators = estimate(~, lf, u)
     p = u.fes.finiteElement.order;
     mesh =  u.fes.mesh;
     trafo = getAffineTransformation(mesh);
diff --git a/examples/goalOrientedAFEM.m b/examples/goalOrientedAFEM.m
index 69297791ea0e34e073a1890ebafb6db2d40e3140..8cc4bc6b5fa244ee2aa2779e346adcd9beddcfa2 100644
--- a/examples/goalOrientedAFEM.m
+++ b/examples/goalOrientedAFEM.m
@@ -95,9 +95,10 @@ title(['goal error estimator over number of dofs'])
 %% local function for residual a posteriori error estimation
 % \eta(T)^2 = h_T^2 * || \Delta u ||_{L^2(T)}^2
 %               + h_T * || [[(Du - fvec) * n]] ||_{L^2(E) \cap \Omega}^2
-function indicators = estimate(blf, lf, u)
-    p = u.fes.finiteElement.order;
-    mesh =  u.fes.mesh;
+function indicators = estimate(~, lf, u)
+    fes = u.fes;
+    p = fes.finiteElement.order;
+    mesh =  fes.mesh;
     trafo = getAffineTransformation(mesh);
     
     % compute volume residual element-wise
@@ -113,9 +114,10 @@ function indicators = estimate(blf, lf, u)
     % compute edge residual edge-wise
     qr = QuadratureRule.ofOrder(max(p-1, 1), '1D');
     f = CompositeFunction(@(p,fvec) p-fvec, Gradient(u), lf.fvec);
+    dirichletBnd = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
     edgeResidual = integrateNormalJump(f, qr, ...
-        @(j) zeros(size(j)), {}, mesh.boundaries{1}, ...    % no jump on dirichlet edges
-        @(j) j.^2, {}, ':');                                % normal jump on inner edges
+        @(j) zeros(size(j)), {}, dirichletBnd, ...    % no jump on dirichlet edges
+        @(j) j.^2, {}, ':');                          % normal jump on inner edges
     
     % combine the resdiuals suitably
     hT = sqrt(trafo.area);