diff --git a/examples/hyperelasticity/driver_hyperelasticity.py b/examples/hyperelasticity/driver_hyperelasticity.py
index 8ee6d40..e2192a7 100644
--- a/examples/hyperelasticity/driver_hyperelasticity.py
+++ b/examples/hyperelasticity/driver_hyperelasticity.py
@@ -42,10 +42,16 @@
 a :code:`ControlCostFunctional` to be compatible with :code:`scipy.optimize`
 
 This driver supports MPI to parallelize the sampling of the parameter field.
-e.g. 
 
-mpirun -n 4 python driver_hyperelasticity.py
+To run deterministic: 
+python driver_hyperelasticity.py
 
+
+To run with mean + variance risk measure:
+python driver_hyperelasticity.py -r mean_var
+
+To run with mean + variance risk measure and parllel sampling (e.g.):
+mpirun -n 4 python driver_hyperelasticity.py -r mean_var
 """
 
 
diff --git a/examples/hyperelasticity/setupHyperelasticityProblem.py b/examples/hyperelasticity/setupHyperelasticityProblem.py
index 85020e9..0135ecc 100644
--- a/examples/hyperelasticity/setupHyperelasticityProblem.py
+++ b/examples/hyperelasticity/setupHyperelasticityProblem.py
@@ -89,13 +89,13 @@ def setup_qoi(mesh, Vh, settings):
     """
     if settings["qoi_type"] == "all":
         u0 = dl.Function(Vh[soupy.STATE]).vector()
-        qoi = soupy.L2MisfitControlQoI(mesh, Vh, u0)
+        qoi = soupy.L2MisfitControlQoI(Vh, u0)
     elif settings["qoi_type"] == "stiffness":
         stiffness = Stiffnessvarf()
-        qoi = soupy.VariationalControlQoI(mesh, Vh, stiffness)
+        qoi = soupy.VariationalControlQoI(Vh, stiffness)
     elif settings["qoi_type"] == "point":
         local_displacement = LocalDisplacementVarf(settings["obs"]["location"], settings["obs"]["width"])
-        qoi = soupy.VariationalControlQoI(mesh, Vh, local_displacement)
+        qoi = soupy.VariationalControlQoI(Vh, local_displacement)
     else:
         raise ValueError("Settings qoi type not available")
 
diff --git a/examples/poisson/driver_poisson_mean.py b/examples/poisson/driver_poisson_mean.py
index 1ad0020..8678510 100644
--- a/examples/poisson/driver_poisson_mean.py
+++ b/examples/poisson/driver_poisson_mean.py
@@ -86,7 +86,7 @@ def boundary(x, on_boundary):
 u_target = dl.Expression("x[1] + sin(k*x[0]) * sin(k*x[1])", k=1.5*np.pi, degree=2, mpi_comm=comm_mesh)
 u_target_function = dl.interpolate(u_target, Vh_STATE)
 u_target_vector = u_target_function.vector()
-qoi = soupy.L2MisfitControlQoI(mesh, Vh, u_target_vector)
+qoi = soupy.L2MisfitControlQoI(Vh, u_target_vector)
 
 # 5. Define the ControlModel
 control_model = soupy.ControlModel(pde, qoi)
diff --git a/examples/semilinear_elliptic/driver_semilinear_cvar.py b/examples/semilinear_elliptic/driver_semilinear_cvar.py
index 2f0fec3..c034ed8 100644
--- a/examples/semilinear_elliptic/driver_semilinear_cvar.py
+++ b/examples/semilinear_elliptic/driver_semilinear_cvar.py
@@ -135,7 +135,7 @@ def print_on_root(print_str, mpi_comm=MPI.COMM_WORLD):
     print_on_root("Making target, QoI, and problem", comm_sampler)
     u_target_expr = get_target(args.target, args.param)
     u_target = dl.interpolate(u_target_expr, Vh[hp.STATE])
-    qoi = soupy.L2MisfitControlQoI(mesh, Vh, u_target.vector())
+    qoi = soupy.L2MisfitControlQoI(Vh, u_target.vector())
     control_model = soupy.ControlModel(pde, qoi)
 
     print_on_root("Making SAA risk measure and cost", comm_sampler)
diff --git a/soupy/modeling/__init__.py b/soupy/modeling/__init__.py
index 8060a14..8eb04f6 100644
--- a/soupy/modeling/__init__.py
+++ b/soupy/modeling/__init__.py
@@ -29,7 +29,7 @@
 
 from .meanVarRiskMeasureSAA import meanVarRiskMeasureSAASettings, MeanVarRiskMeasureSAA 
 
-from .penalization import Penalization, L2Penalization, WeightedL2Penalization, MultiPenalization
+from .penalization import Penalization, L2Penalization, WeightedL2Penalization, MultiPenalization, VariationalPenalization
 
 from .riskMeasure import RiskMeasure
 
diff --git a/soupy/modeling/controlQoI.py b/soupy/modeling/controlQoI.py
index 8dfba0e..3b326a6 100644
--- a/soupy/modeling/controlQoI.py
+++ b/soupy/modeling/controlQoI.py
@@ -100,11 +100,10 @@ class VariationalControlQoI(ControlQoI):
     """
     Class for a QoI defined by its variational form 
     """
-    def __init__(self, mesh, Vh, form_handler):
+    def __init__(self, Vh, form_handler):
         """
         Constructor
 
-        :param mesh: The mesh object
         :param Vh: List of function spaces for the state, parameter, 
             adjoint, and optimization variables
         :type Vh: list of :py:class:`dolfin.FunctionSpace`
@@ -112,7 +111,6 @@ def __init__(self, mesh, Vh, form_handler):
             :code:`__call__` method that takes as input the state, parameter, and control variables
             as functions and returns the variational form 
         """
-        self.mesh = mesh
         self.Vh = Vh
         self.x = [dl.Function(Vh[STATE]).vector(), dl.Function(Vh[PARAMETER]).vector(),
                   dl.Function(Vh[ADJOINT]).vector(), dl.Function(Vh[CONTROL]).vector()]
@@ -254,18 +252,16 @@ class L2MisfitControlQoI(ControlQoI):
 
     where :math:`u_d` is the reference state.
     """
-    def __init__(self, mesh, Vh, ud):
+    def __init__(self, Vh, ud):
         """
         Constructor
 
-        :param mesh: The mesh object
         :param Vh: List of function spaces for the state, parameter, 
             adjoint, and optimization variables
         :type Vh: list of :py:class:`dolfin.FunctionSpace`
         :param ud: The reference state as a vector
         :type ud: :py:class:`dolfin.Vector`
         """
-        self.mesh = mesh
         self.Vh = Vh
         self.x = [dl.Function(Vh[STATE]).vector(), dl.Function(Vh[PARAMETER]).vector(),
                   dl.Function(Vh[ADJOINT]).vector(), dl.Function(Vh[CONTROL]).vector()]
diff --git a/soupy/modeling/penalization.py b/soupy/modeling/penalization.py
index 6991c26..492bd1e 100644
--- a/soupy/modeling/penalization.py
+++ b/soupy/modeling/penalization.py
@@ -33,6 +33,80 @@ def hessian(self, z, zhat, out):
         raise NotImplementedError("Child class should implement hessian")
 
 
+
+class VariationalPenalization(Penalization):
+    """
+    Penalization given by a variational form in terms of the control variable :math:`z`
+    """
+    def __init__(self, Vh, form_handler):
+        """
+        Constructor: 
+
+        :param Vh: List of function spaces for the state, parameter, 
+            adjoint, and optimization variables
+        :type Vh: list of :py:class:`dolfin.FunctionSpace`
+        :param form_handler: The form handler for the penalization with a
+            :code:`__call__` method that takes the control variable :code:`z` as a :py:class:`dolfin.Function` and returns the variational form. 
+        """
+
+        self.Vh = Vh 
+        self.form_handler = form_handler
+
+        self.z_fun = dl.Function(self.Vh[CONTROL])
+        self.z = self.z_fun.vector()
+
+        self.dz_fun = dl.Function(self.Vh[CONTROL])
+        self.dz = self.dz_fun.vector()
+
+        self.z_trial = dl.TrialFunction(self.Vh[CONTROL])
+        self.z_test = dl.TestFunction(self.Vh[CONTROL])
+
+
+        self.penalization_form = self.form_handler(self.z_fun)
+        self.grad_form = dl.derivative(self.penalization_form, self.z_fun, self.z_test)
+
+    def init_vector(self, z):
+        if isinstance(z, AugmentedVector):
+            pass 
+        else:
+            z.init(self.z_fun.vector().local_range())
+
+    def cost(self, z):
+        if isinstance(z, AugmentedVector):
+            z = z.get_vector()
+
+        self.z.zero()
+        self.z.axpy(1.0, z)
+        return dl.assemble(self.penalization_form)
+
+    def grad(self, z, out):
+        if isinstance(z, AugmentedVector):
+            out.set_scalar(0)
+            z = z.get_vector()
+            out = out.get_vector()
+
+        self.z.zero()
+        self.z.axpy(1.0, z)
+        dl.assemble(self.grad_form, tensor=out)
+
+    def hessian(self, z, zhat, out):
+        if isinstance(z, AugmentedVector):
+            out.set_scalar(0)
+            z = z.get_vector()
+            zhat = zhat.get_vector()
+            out = out.get_vector()
+
+        self.z.zero()
+        self.z.axpy(1.0, z)
+
+        self.dz.zero()
+        self.dz.axpy(1.0, zhat)
+
+        hessian_action = dl.derivative(self.grad_form, self.z_fun, self.dz_fun)
+        dl.assemble(hessian_action, tensor=out)
+
+
+
 class MultiPenalization(Penalization):
     """
     Penalization term for the sum of individual penalties
@@ -61,6 +135,8 @@ def __init__(self, Vh, penalization_list, alpha_list=None):
         else:
             self.alpha_list = [1.0] * len(self.penalization_list)
 
+    def init_vector(self, z):
+        self.penalization_list[0].init_vector(z)
     
     def cost(self, z):
         """
diff --git a/soupy/test/ptest_meanVarRiskMeasureSAA.py b/soupy/test/ptest_meanVarRiskMeasureSAA.py
index ce5fc3f..e0a3f94 100644
--- a/soupy/test/ptest_meanVarRiskMeasureSAA.py
+++ b/soupy/test/ptest_meanVarRiskMeasureSAA.py
@@ -75,7 +75,7 @@ def testCostValue(self):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(1.0))**2*dl.dx 
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         risk = MeanVarRiskMeasureSAA(model, prior, comm_sampler=self.comm_sampler)
@@ -96,7 +96,7 @@ def _setup_pde_and_distributions(self, is_fwd_linear):
         return pde, prior, control_dist
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 
diff --git a/soupy/test/ptest_scipyCostWrapper.py b/soupy/test/ptest_scipyCostWrapper.py
index fbe12d6..cf2dbb2 100644
--- a/soupy/test/ptest_scipyCostWrapper.py
+++ b/soupy/test/ptest_scipyCostWrapper.py
@@ -93,7 +93,7 @@ def _sample_control(self, z):
         self.collective.bcast(z, root=0)
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 
diff --git a/soupy/test/ptest_superquantileSAA.py b/soupy/test/ptest_superquantileSAA.py
index fce940d..d1c53b4 100644
--- a/soupy/test/ptest_superquantileSAA.py
+++ b/soupy/test/ptest_superquantileSAA.py
@@ -103,7 +103,7 @@ def _setup_pde_and_distributions(self, is_fwd_linear):
         return pde, prior, control_dist
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 
@@ -279,7 +279,7 @@ def computeSuperquantileValue(self, sample_size, beta):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         rm_settings = superquantileRiskMeasureSAASettings()
diff --git a/soupy/test/test_SAACostFunctional.py b/soupy/test/test_SAACostFunctional.py
index e266b35..62ef2fa 100644
--- a/soupy/test/test_SAACostFunctional.py
+++ b/soupy/test/test_SAACostFunctional.py
@@ -68,7 +68,7 @@ def _setup_pde_and_distributions(self, is_fwd_linear):
         return pde, prior, control_dist
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 
diff --git a/soupy/test/test_controlCostFunctional.py b/soupy/test/test_controlCostFunctional.py
index d490633..b8af44b 100644
--- a/soupy/test/test_controlCostFunctional.py
+++ b/soupy/test/test_controlCostFunctional.py
@@ -76,7 +76,7 @@ def testCostValue(self):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(2.0))**2*dl.dx 
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
         cost = DeterministicControlCostFunctional(model, prior, penalization=None)
         z = model.generate_vector(CONTROL)
@@ -102,7 +102,7 @@ def finiteDifferenceCheck(self, qoi_varf, is_fwd_linear=True):
         settings['LINEAR'] = is_fwd_linear
 
         pde, prior, control_dist = setupPoissonPDEProblem(self.Vh, settings) 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         cost = DeterministicControlCostFunctional(model, prior, penalization=None)
 
@@ -196,7 +196,7 @@ def testSavedSolution(self):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
         cost = DeterministicControlCostFunctional(model, prior, penalization=None)
 
diff --git a/soupy/test/test_controlQoI.py b/soupy/test/test_controlQoI.py
index 1a0f1ee..e8562f6 100644
--- a/soupy/test/test_controlQoI.py
+++ b/soupy/test/test_controlQoI.py
@@ -58,7 +58,7 @@ def finiteDifferenceGradient(self, nonlinear_qoi, i):
 
         # ud = dl.Function(self.Vh[STATE])
         # ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=2))
-        # nonlinear_qoi = VariationalControlQoI(self.mesh, self.Vh, varf)
+        # nonlinear_qoi = VariationalControlQoI(self.Vh, varf)
 
 
         u_fun = dl.Function(self.Vh[STATE])
@@ -137,7 +137,7 @@ def testFiniteDifference(self):
         ud = dl.Function(self.Vh[STATE])
         ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=2))
         for varf in varfs:
-            nonlinear_qoi = VariationalControlQoI(self.mesh, self.Vh, varf)
+            nonlinear_qoi = VariationalControlQoI(self.Vh, varf)
             for i in [STATE, PARAMETER, CONTROL]:
                 self.finiteDifferenceGradient(nonlinear_qoi, i)
                 for j in [STATE, PARAMETER, CONTROL]:
@@ -150,7 +150,7 @@ def testVariationalControlQoI(self):
         ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=5))
         chi = dl.Expression("0.5", degree=5)
         l2varf = L2MisfitVarfHandler(ud, chi=chi)
-        l2qoi = VariationalControlQoI(self.mesh, self.Vh, l2varf)
+        l2qoi = VariationalControlQoI(self.Vh, l2varf)
         u = dl.Function(self.Vh[STATE])
         m = dl.Function(self.Vh[PARAMETER])
         z = dl.Function(self.Vh[CONTROL])
@@ -166,7 +166,7 @@ def testVariationalControlQoI(self):
     def testL2MisfitQoI(self):
         ud = dl.Function(self.Vh[STATE])
         ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=2))
-        nonlinear_qoi = L2MisfitControlQoI(self.mesh, self.Vh, ud.vector())
+        nonlinear_qoi = L2MisfitControlQoI(self.Vh, ud.vector())
         for i in [STATE, PARAMETER, CONTROL]:
             self.finiteDifferenceGradient(nonlinear_qoi, i)
             for j in [STATE, PARAMETER, CONTROL]:
diff --git a/soupy/test/test_meanVarRiskMeasureStochastic.py b/soupy/test/test_meanVarRiskMeasureStochastic.py
index 3e74c21..4653c8f 100644
--- a/soupy/test/test_meanVarRiskMeasureStochastic.py
+++ b/soupy/test/test_meanVarRiskMeasureStochastic.py
@@ -63,7 +63,7 @@ def testCostValue(self):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(1.0))**2*dl.dx 
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         risk = MeanVarRiskMeasureStochastic(model, prior)
@@ -91,7 +91,7 @@ def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         rm_settings = meanVarRiskMeasureStochasticSettings()
diff --git a/soupy/test/test_scipyCostWrapper.py b/soupy/test/test_scipyCostWrapper.py
index c27069b..894940d 100644
--- a/soupy/test/test_scipyCostWrapper.py
+++ b/soupy/test/test_scipyCostWrapper.py
@@ -84,7 +84,7 @@ def _sample_control(self, z):
         self.control_prior.sample(self.control_noise, z)
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 
diff --git a/soupy/test/test_stochasticCostFunctional.py b/soupy/test/test_stochasticCostFunctional.py
index 7451d73..191927c 100644
--- a/soupy/test/test_stochasticCostFunctional.py
+++ b/soupy/test/test_stochasticCostFunctional.py
@@ -65,7 +65,7 @@ def costValueCheck(self, sample_size, use_penalization):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(1.0))**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
         risk = MeanVarRiskMeasureStochastic(model, prior)
         alpha = 2.0
@@ -124,7 +124,7 @@ def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         rm_settings = meanVarRiskMeasureStochasticSettings()
diff --git a/soupy/test/test_variationalPenalization.py b/soupy/test/test_variationalPenalization.py
new file mode 100644
index 0000000..f5f4855
--- /dev/null
+++ b/soupy/test/test_variationalPenalization.py
@@ -0,0 +1,187 @@
+# Copyright (c) 2023, The University of Texas at Austin 
+# & Georgia Institute of Technology
+#
+# All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the SOUPy package. For more information see
+# https://github.com/hippylib/soupy/
+#
+# SOUPy is free software; you can redistribute it and/or modify it under the
+# terms of the GNU General Public License (as published by the Free
+# Software Foundation) version 3.0 dated June 2007.
+
+import unittest
+
+import numpy as np
+import dolfin as dl
+
+import sys, os 
+sys.path.append(os.environ.get('HIPPYLIB_PATH'))
+import hippylib as hp
+
+import logging
+logging.getLogger('FFC').setLevel(logging.WARNING)
+
+sys.path.append('../../')
+from soupy import STATE, PARAMETER, ADJOINT, CONTROL, L2Penalization, VariationalPenalization, AugmentedVector
+
+class TestPenalization(unittest.TestCase):
+    def setUp(self):
+        self.mesh = dl.UnitSquareMesh(32, 32)
+
+    def finiteDifferenceGradient(self, penalty, z, dz, delta=1e-5, grad_tol=1e-3):
+        """
+        Finite difference check in direction dz 
+        """
+        grad = dl.Vector()
+        penalty.init_vector(grad)
+        if isinstance(z, AugmentedVector):
+            grad = AugmentedVector(grad)
+
+        p0 = penalty.cost(z)
+        penalty.grad(z, grad)
+
+        z.axpy(delta, dz)
+        p1 = penalty.cost(z)
+
+        finite_diff = (p1 - p0)/delta
+        symbolic_derivative = grad.inner(dz)
+        print("Finite difference derivative: %g" %(finite_diff))
+        print("Implemented derivative: %g" %(symbolic_derivative))
+        self.assertTrue(abs(finite_diff - symbolic_derivative) < grad_tol)
+
+        # Make sure that for augmented vectors, the derivative component in t is always zero
+        if isinstance(z, AugmentedVector):
+            self.assertTrue(grad.get_scalar() < 1e-12)
+
+
+    def finiteDifferenceHessian(self, penalty, z, zhat, delta=1e-5, hess_tol=1e-3):
+        """
+        Finite difference hessian check in direction dz 
+        """
+        grad0 = dl.Vector()
+        grad1 = dl.Vector()
+        Hzhat = dl.Vector()
+
+        penalty.init_vector(grad0)
+        penalty.init_vector(grad1)
+        penalty.init_vector(Hzhat)
+
+        if isinstance(z, AugmentedVector):
+            grad0 = AugmentedVector(grad0)
+            grad1 = AugmentedVector(grad1)
+            Hzhat = AugmentedVector(Hzhat)
+
+        penalty.grad(z, grad0)
+        penalty.hessian(z, zhat, Hzhat)
+
+        z.axpy(delta, zhat)
+        penalty.grad(z, grad1)
+
+        Hzhat_fd = (grad1.get_local() - grad0.get_local())/delta
+        hess_diff = Hzhat_fd - Hzhat.get_local()
+
+        if np.linalg.norm(Hzhat.get_local()) > 0:
+            hess_err = np.linalg.norm(hess_diff)/np.linalg.norm(Hzhat.get_local())
+        else:
+            hess_err = np.linalg.norm(hess_diff)
+
+        print("Hessian difference: %g" %(hess_err))
+        self.assertTrue(hess_err < hess_tol)
+
+        # Make sure that for augmented vectors, the derivative component in t is always zero
+        if isinstance(z, AugmentedVector):
+            self.assertTrue(Hzhat.get_scalar() < 1e-12)
+
+
+    def compareFunctionValue(self, penalty, z, p_true, tol=1e-3):
+        p = penalty.cost(z)
+        p_err = abs(p - p_true)
+        print("Evaluation: %g" %(p))
+        print("Truth: %g" %(p_true))
+        self.assertTrue(p_err < tol)
+
+
+    def testL2Penalization(self):
+        """
+        Compare the VariationalPenalization against the L2Penalization class \
+            Tests value for both :py:class:`dolfin.Vector` and \
+            :py:class:`soupy.AugmentedVector` as inputs
+
+        """
+
+        print("Test function L2 Penalization")
+        Vh_CONTROL = dl.FunctionSpace(self.mesh, "CG", 2)
+        Vh = [None, None, None, Vh_CONTROL]
+
+        alpha = 2.0
+        n_test = 3
+        penalty_L2 = L2Penalization(Vh, alpha)
+
+        def l2_varf(z):
+            return dl.Constant(alpha) * z**2 * dl.dx 
+
+        variational_l2 = VariationalPenalization(Vh, l2_varf)
+
+        # Define the points 
+        z_fun = dl.Function(Vh_CONTROL)
+        z = z_fun.vector()
+        dz_fun = dl.Function(Vh_CONTROL)
+        dz = dz_fun.vector()
+
+        for i in range(n_test):
+            hp.parRandom.normal(1.0, z)
+            p_true = penalty_L2.cost(z)
+            self.compareFunctionValue(variational_l2, z, p_true, tol=1e-8)
+
+            # Also check if input is an augmented vector 
+            zt = AugmentedVector(z)
+            self.compareFunctionValue(variational_l2, zt, p_true, tol=1e-8)
+
+
+
+    def testFiniteDifference(self):
+        print("Test finite difference using total variation")
+        Vh_CONTROL = dl.FunctionSpace(self.mesh, "CG", 2)
+        Vh = [None, None, None, Vh_CONTROL]
+
+        penalty_weight = 2.0
+        tv_smoothing = 1e-3
+        delta = 1e-5
+        tol = 1e-3
+
+        def total_variation(z):
+            return dl.Constant(penalty_weight) * dl.sqrt(dl.inner(dl.grad(z), dl.grad(z)) + dl.Constant(tv_smoothing))*dl.dx 
+        
+        tv_penalty = VariationalPenalization(Vh, total_variation)
+
+        # Define the points 
+        z_fun = dl.Function(Vh_CONTROL)
+        z = z_fun.vector()
+        dz_fun = dl.Function(Vh_CONTROL)
+        dz = dz_fun.vector()
+
+        z_expression = dl.Expression("cos(k*x[0])*x[1] + exp(-pow(x[0]-0.5, 2) - pow(x[1]-0.5, 2))", 
+            k=1, degree=2)
+        dz_expression = dl.Expression("cos(k*x[0])*x[1] + exp(-pow(x[0]-0.5, 2) - pow(x[1]-0.5, 2))", 
+            k=2.88, degree=5)
+
+        n_test = 3
+        # check for different functions 
+        for i in range(n_test):
+            z_expression.k = i + 1
+            dz_expression.k = i + 2.88
+            z_fun.interpolate(z_expression)
+            dz_fun.interpolate(dz_expression)
+
+            self.finiteDifferenceGradient(tv_penalty, z, dz, delta=delta, grad_tol=tol)
+            self.finiteDifferenceHessian(tv_penalty, z, dz, delta=delta, hess_tol=tol)
+
+            self.finiteDifferenceGradient(tv_penalty, AugmentedVector(z), AugmentedVector(dz), delta=delta, grad_tol=tol)
+            self.finiteDifferenceHessian(tv_penalty, AugmentedVector(z), AugmentedVector(dz), delta=delta, hess_tol=tol)
+
+if __name__ == "__main__":
+    unittest.main()
+
+