diff --git a/test/tests/indicators/laplacian_jump_indicator/biharmonic.i b/test/tests/indicators/laplacian_jump_indicator/biharmonic.i
new file mode 100644
index 000000000000..858f08822c9b
--- /dev/null
+++ b/test/tests/indicators/laplacian_jump_indicator/biharmonic.i
@@ -0,0 +1,123 @@
+[GlobalParams]
+  # Parameters used by Functions.
+  vars   = 'c'
+  vals   = '50'
+[]
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  xmin = -.5
+  xmax = .5
+  ymin = -.5
+  ymax = .5
+  nx = 10
+  ny = 10
+[]
+
+[Variables]
+  [./u]
+    order = THIRD
+    family = HERMITE
+  [../]
+[]
+
+[Kernels]
+  [./biharmonic]
+    type = Biharmonic
+    variable = u
+  [../]
+  [./body_force]
+    type = BodyForce
+    variable = u
+    function = forcing_func
+  [../]
+[]
+
+[BCs]
+  active = 'all_value all_flux'
+  [./all_value]
+    type = FunctionPenaltyDirichletBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+  [./all_flux]
+    type = FunctionPenaltyFluxBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+  [./all_laplacian]
+    type = BiharmonicLapBC
+    variable = u
+    boundary = 'left right top bottom'
+    laplacian_function = lapu_func
+  [../]
+[]
+
+[Adaptivity]
+  [./Indicators]
+    [./error]
+      type = LaplacianJumpIndicator
+      variable = u
+      scale_by_flux_faces = true
+    [../]
+  [../]
+[]
+
+[Executioner]
+  type = Steady
+
+  # Note: the unusually tight tolerances here are due to the penalty
+  # BCs (currently the only way of accurately Dirichlet boundary
+  # conditions on Hermite elements in MOOSE).
+  nl_rel_tol = 1.e-15
+  l_tol = 1.e-15
+
+  # We have exact Jacobians
+  solve_type = 'NEWTON'
+
+  # Use 6x6 quadrature to ensure the forcing function is integrated
+  # accurately.
+  [./Quadrature]
+    type = GAUSS
+    order = ELEVENTH
+  [../]
+[]
+
+[Functions]
+  [./u_func]
+    type   = ParsedGradFunction
+    value  = 'exp(-c*(x^2+y^2))'
+    grad_x = '-2*c*exp(-c*(x^2+y^2))*x'
+    grad_y = '-2*c*exp(-c*(x^2+y^2))*y'
+  [../]
+  [./lapu_func]
+    type   = ParsedFunction
+    value  = '4*c*(c*(x^2+y^2) - 1)*exp(-c*(x^2+y^2))'
+  [../]
+  [./forcing_func]
+    type   = ParsedFunction
+    value  = '16*c^2*(c^2*(x^2+y^2)^2 - 4*c*(x^2+y^2) + 2)*exp(-c*(x^2+y^2))'
+  [../]
+[]
+
+[Postprocessors]
+  [./l2_error]
+    type = ElementL2Error
+    variable = u
+    function = u_func
+  [../]
+  [./h1_error]
+    type = ElementH1Error
+    variable = u
+    function = u_func
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/test/tests/indicators/laplacian_jump_indicator/biharmonic.py b/test/tests/indicators/laplacian_jump_indicator/biharmonic.py
new file mode 100755
index 000000000000..7be12e137c39
--- /dev/null
+++ b/test/tests/indicators/laplacian_jump_indicator/biharmonic.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+from sympy import *
+from fractions import Fraction
+
+r, theta, t, c = sympify('r, theta, t, c')
+
+# A simple Gaussian function (all derivatives wrt theta are 0).
+# This script should print the following:
+#
+# Steady version: f = du/dt + Biharmonic(u)
+# lapu = 4*c*(c*r**2 - 1)*exp(-c*r**2)
+# bihu = 16*c**2*(c**2*r**4 - 4*c*r**2 + 2)*exp(-c*r**2)
+#
+# Time dependent version: f = du/dt + Biharmonic(u)
+u = exp(-t) * exp(-c * r**2)
+
+ut = diff(u,t)
+ur = diff(u,r)
+urr = diff(ur,r)
+u_theta = diff(u,theta)
+u_theta_theta = diff(u_theta,theta)
+
+# Gradient component
+gradu_r=ur
+gradu_t=(1/r)*u_theta
+# print('gradu_r = {}'.format(gradu_r))
+# print('gradu_t = {}'.format(gradu_t))
+
+lapu = simplify(urr + ur/r + u_theta_theta/r**2)
+print('lapu = {}'.format(lapu))
+
+# The biharmonic operator in polar coodinates has a bunch of nested
+# derivatives in r. Here we are ignoring the theta derivatives because
+# our Gaussian does not have any theta dependence.
+# https://en.wikipedia.org/wiki/Biharmonic_equation
+bihu = simplify((1/r) * diff(r * diff((1/r) * diff(r*diff(u,r),r),r),r))
+print('bihu = {}'.format(bihu))
+
+# 16*c**2*(c**2*r**4 - 4*c*r**2 + 2)*exp(-c*r**2 - t) - exp(-t)*exp(-c*r**2)
+f = ut + bihu
+print('Time-dependent forcing function = {}'.format(f))
diff --git a/test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_dirichlet.csv b/test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_dirichlet.csv
new file mode 100644
index 000000000000..b2a7e3184f73
--- /dev/null
+++ b/test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_dirichlet.csv
@@ -0,0 +1,4 @@
+.1,8.411260e-04,3.207318e-02
+.05,5.885662e-05,4.515238e-03
+.025,3.721949e-06,5.738347e-04
+.0125,2.332737e-07,7.202340e-05
diff --git a/test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_weak_bc.csv b/test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_weak_bc.csv
new file mode 100644
index 000000000000..88eb7b4adb92
--- /dev/null
+++ b/test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_weak_bc.csv
@@ -0,0 +1,4 @@
+.1,8.411280e-04,3.207327e-02
+.05,5.886178e-05,4.515250e-03
+.025,3.722349e-06,5.738352e-04
+.0125,2.332996e-07,7.202341e-05
diff --git a/test/tests/indicators/laplacian_jump_indicator/biharmonic_transient.i b/test/tests/indicators/laplacian_jump_indicator/biharmonic_transient.i
new file mode 100644
index 000000000000..44c523ba328e
--- /dev/null
+++ b/test/tests/indicators/laplacian_jump_indicator/biharmonic_transient.i
@@ -0,0 +1,122 @@
+[GlobalParams]
+  # Parameters used by Functions.
+  vars   = 'c'
+  vals   = '50'
+[]
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  xmin = -.5
+  xmax = .5
+  ymin = -.5
+  ymax = .5
+  nx = 10
+  ny = 10
+[]
+
+[Variables]
+  [./u]
+    order = THIRD
+    family = HERMITE
+  [../]
+[]
+
+[Kernels]
+  [./biharmonic]
+    type = Biharmonic
+    variable = u
+  [../]
+  [./body_force]
+    type = BodyForce
+    variable = u
+    function = forcing_func
+  [../]
+[]
+
+[BCs]
+  [./all_value]
+    type = FunctionPenaltyDirichletBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+  [./all_flux]
+    type = FunctionPenaltyFluxBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+[]
+
+[Adaptivity]
+  [./Indicators]
+    [./error]
+      type = LaplacianJumpIndicator
+      variable = u
+      scale_by_flux_faces = true
+    [../]
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 4
+  dt = 0.1
+
+  # Note: the unusually tight tolerances here are due to the penalty
+  # BCs (currently the only way of accurately Dirichlet boundary
+  # conditions on Hermite elements in MOOSE).
+  nl_rel_tol = 1.e-15
+  l_tol = 1.e-15
+
+  # We have exact Jacobians
+  solve_type = 'NEWTON'
+
+  # Use 6x6 quadrature to ensure the forcing function is integrated
+  # accurately.
+  [./Quadrature]
+    type = GAUSS
+    order = ELEVENTH
+  [../]
+[]
+
+[Functions]
+  [./u_func]
+    type   = ParsedGradFunction
+    value  = 'exp(-c*(x^2+y^2))*exp(-t)'
+    grad_x = '-2*c*exp(-c*(x^2+y^2))*x*exp(-t)'
+    grad_y = '-2*c*exp(-c*(x^2+y^2))*y*exp(-t)'
+  [../]
+  [./forcing_func]
+    type   = ParsedFunction
+    value  = '16*c^2*(c^2*(x^2+y^2)^2 - 4*c*(x^2+y^2) + 2)*exp(-c*(x^2+y^2))*exp(-t)'
+  [../]
+[]
+
+[ICs]
+  [./u_ic]
+    type = FunctionIC
+    function = u_func
+    variable = u
+  [../]
+[]
+
+[Postprocessors]
+  [./l2_error]
+    type = ElementL2Error
+    variable = u
+    function = u_func
+  [../]
+  [./h1_error]
+    type = ElementH1Error
+    variable = u
+    function = u_func
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_out.e b/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_out.e
new file mode 100644
index 000000000000..6073b917e60a
Binary files /dev/null and b/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_out.e differ
diff --git a/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_transient_out.e b/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_transient_out.e
new file mode 100644
index 000000000000..26fd4008d033
Binary files /dev/null and b/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_transient_out.e differ
diff --git a/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_weak_bc_out.e b/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_weak_bc_out.e
new file mode 100644
index 000000000000..c990cfbad23d
Binary files /dev/null and b/test/tests/indicators/laplacian_jump_indicator/gold/biharmonic_weak_bc_out.e differ
diff --git a/test/tests/indicators/laplacian_jump_indicator/gold/laplacian_jump_indicator_test_out.e b/test/tests/indicators/laplacian_jump_indicator/gold/laplacian_jump_indicator_test_out.e
deleted file mode 100644
index 84ecbadd4705..000000000000
Binary files a/test/tests/indicators/laplacian_jump_indicator/gold/laplacian_jump_indicator_test_out.e and /dev/null differ
diff --git a/test/tests/indicators/laplacian_jump_indicator/gold/laplacian_jump_transient_out.e b/test/tests/indicators/laplacian_jump_indicator/gold/laplacian_jump_transient_out.e
deleted file mode 100644
index 685656bcd939..000000000000
Binary files a/test/tests/indicators/laplacian_jump_indicator/gold/laplacian_jump_transient_out.e and /dev/null differ
diff --git a/test/tests/indicators/laplacian_jump_indicator/laplacian_jump_indicator_test.i b/test/tests/indicators/laplacian_jump_indicator/laplacian_jump_indicator_test.i
deleted file mode 100644
index c61a20134e1a..000000000000
--- a/test/tests/indicators/laplacian_jump_indicator/laplacian_jump_indicator_test.i
+++ /dev/null
@@ -1,68 +0,0 @@
-[Mesh]
-  type = GeneratedMesh
-  dim = 2
-  nx = 10
-  ny = 10
-  nz = 10
-[]
-
-[Variables]
-  [./u]
-    order = THIRD
-    family = HERMITE
-  [../]
-[]
-
-[Kernels]
-  [./diff]
-    type = Diffusion
-    variable = u
-  [../]
-  [./conv]
-    type = Convection
-    variable = u
-    velocity = '1 0 0'
-  [../]
-[]
-
-[BCs]
-  [./left]
-    type = DirichletBC
-    variable = u
-    boundary = left
-    value = 0
-  [../]
-  [./right]
-    type = DirichletBC
-    variable = u
-    boundary = right
-    value = 1
-  [../]
-  [./divbc]
-    type = DivergenceBC
-    variable = u
-    boundary = 'left right'
-  [../]
-[]
-
-[Adaptivity]
-  [./Indicators]
-    [./error]
-      type = LaplacianJumpIndicator
-      variable = u
-      scale_by_flux_faces = true
-    [../]
-  [../]
-[]
-
-[Executioner]
-  type = Steady
-  nl_rel_tol = 1.e-11
-
-  # Preconditioned JFNK (default)
-  solve_type = 'PJFNK'
-[]
-
-[Outputs]
-  exodus = true
-[]
diff --git a/test/tests/indicators/laplacian_jump_indicator/laplacian_jump_transient.i b/test/tests/indicators/laplacian_jump_indicator/laplacian_jump_transient.i
deleted file mode 100644
index 1cbabb569ea1..000000000000
--- a/test/tests/indicators/laplacian_jump_indicator/laplacian_jump_transient.i
+++ /dev/null
@@ -1,70 +0,0 @@
-[Mesh]
-  type = GeneratedMesh
-  dim = 2
-  nx = 4
-  ny = 4
-  nz = 0
-[]
-
-[Variables]
-  [./u]
-    order = THIRD
-    family = HERMITE
-  [../]
-[]
-
-[Kernels]
-  [./diff]
-    type = Diffusion
-    variable = u
-  [../]
-  [./conv]
-    type = Convection
-    variable = u
-    velocity = '1 0 0'
-  [../]
-  [./time]
-    type = TimeDerivative
-    variable = u
-  [../]
-[]
-
-[BCs]
-  [./left]
-    type = DirichletBC
-    variable = u
-    boundary = left
-    value = 0
-  [../]
-  [./right]
-    type = DirichletBC
-    variable = u
-    boundary = right
-    value = 10
-  [../]
-  [./divbc]
-    type = DivergenceBC
-    variable = u
-    boundary = 'left right'
-  [../]
-[]
-
-[Adaptivity]
-  [./Indicators]
-    [./error]
-      type = LaplacianJumpIndicator
-      variable = u
-      scale_by_flux_faces = true
-    [../]
-  [../]
-[]
-
-[Executioner]
-  type = Transient
-  num_steps = 4
-  dt = 0.1
-[]
-
-[Outputs]
-  exodus = true
-[]
diff --git a/test/tests/indicators/laplacian_jump_indicator/plot.py b/test/tests/indicators/laplacian_jump_indicator/plot.py
new file mode 100755
index 000000000000..1ce9470fb32b
--- /dev/null
+++ b/test/tests/indicators/laplacian_jump_indicator/plot.py
@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+import matplotlib.pyplot as plt
+import numpy as np
+
+filenames = ['biharmonic_centered_dirichlet.csv',
+             'biharmonic_centered_weak_bc.csv']
+
+for filename in filenames:
+    fig = plt.figure()
+    ax1 = fig.add_subplot(111)
+    data = np.genfromtxt(filename, delimiter=',')
+
+    # Extract vectors
+    logh = np.log10(data[:,0])
+    log_L2 = np.log10(data[:,1])
+    log_H1 = np.log10(data[:,2])
+
+    # Compute linear fits
+    l2_fit = np.polyfit(logh, log_L2, 1)
+    h1_fit = np.polyfit(logh, log_H1, 1)
+
+    # Make log-log plots
+    ax1.plot(logh, log_L2, linewidth=2, marker='o', label=r'$L^2$ error')
+    ax1.plot(logh, np.poly1d(l2_fit)(logh), linestyle='--', color='black')
+    ax1.plot(logh, log_H1, linewidth=2, marker='o', label=r'$H^1$ error')
+    ax1.plot(logh, np.poly1d(h1_fit)(logh), linestyle='--', color='black')
+
+    # Put text labels on plot
+    ax1.text(-1.6, -6, '{:.2f}'.format(l2_fit[0]))
+    ax1.text(-1.6, -3.75, '{:.2f}'.format(h1_fit[0]))
+
+    # Add axis labels, legend, and print
+    ax1.set_xlabel('log(h)')
+    ax1.legend(loc='upper left')
+    plt.savefig(filename.rsplit( ".", 1)[0] + '.pdf')
diff --git a/test/tests/indicators/laplacian_jump_indicator/tests b/test/tests/indicators/laplacian_jump_indicator/tests
index f62846e9d33b..c001c5dc4cac 100644
--- a/test/tests/indicators/laplacian_jump_indicator/tests
+++ b/test/tests/indicators/laplacian_jump_indicator/tests
@@ -1,14 +1,18 @@
 [Tests]
-  [./test]
+  [./test_biharmonic]
     type = 'Exodiff'
-    input = 'laplacian_jump_indicator_test.i'
-    exodiff = 'laplacian_jump_indicator_test_out.e'
-    scale_refine = 2
+    input = 'biharmonic.i'
+    exodiff = 'biharmonic_out.e'
   [../]
-
-  [./test_transient]
+  [./test_biharmonic_weak_bc]
     type = 'Exodiff'
-    input = 'laplacian_jump_transient.i'
-    exodiff = 'laplacian_jump_transient_out.e'
+    input = 'biharmonic.i'
+    exodiff = 'biharmonic_weak_bc_out.e'
+    cli_args = "BCs/active='all_value all_laplacian' Outputs/file_base=biharmonic_weak_bc_out"
+  [../]
+  [./test_biharmonic_transient]
+    type = 'Exodiff'
+    input = 'biharmonic_transient.i'
+    exodiff = 'biharmonic_transient_out.e'
   [../]
 []