diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 2795d33094715fa5976099214d93e386701a23be..592e606b1baf174845cf4ed741e6f76b73efdf8c 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -433,33 +433,43 @@ def generate_localoperator_kernels(formdata, data):
 
                 # Maybe add numerical differentiation
                 if get_option("numerical_jacobian"):
+                    # Include headers for numerical methods
                     include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
 
+                    # Numerical jacobian base class
                     _, loptype = class_type_from_cache("operator")
                     which = ufl_measure_to_pdelab_measure(measure)
-
-                    # Numerical jacobian methods
                     base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")
 
-                    # Add the initializer list for that base class
+                    # Numerical jacobian initializer list
                     ini = name_initree_member()
                     ini_constructor = name_initree_constructor_param()
-
-                    # Numerical jacobian methods
                     initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
                                      ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
                                      classtag="operator",
                                      )
 
+                    # In the case of matrix free operator evaluation we need jacobian apply methods
                     if get_option("matrix_free"):
-                        # Numeical jacobian apply base class
-                        base_class("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype), classtag="operator")
-
-                        # Numerical jacobian apply initializer list
-                        initializer_list("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype),
-                                         ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
-                                         classtag="operator",
-                                         )
+                        from dune.perftool.pdelab.driver import is_linear
+                        if is_linear(form):
+                            # Numeical jacobian apply base class
+                            base_class("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype), classtag="operator")
+
+                            # Numerical jacobian apply initializer list
+                            initializer_list("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype),
+                                             ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
+                                              classtag="operator",
+                                             )
+                        else:
+                            # Numerical nonlinear jacobian apply base class
+                            base_class("Dune::PDELab::NumericalNonlinearJacobianApply{}<{}>".format(which, loptype), classtag="operator")
+
+                            # Numerical nonlinear jacobian apply initializer list
+                            initializer_list("Dune::PDELab::NumericalNonlinearJacobianApply{}<{}>".format(which, loptype),
+                                             ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
+                                              classtag="operator",
+                                             )
 
                 operator_kernels[(measure, 'residual')] = kernel