diff --git a/cmake/modules/DuneCodegenMacros.cmake b/cmake/modules/DuneCodegenMacros.cmake
index 91e48d73f4c79ccc6a7258071dbbfeaf8a27a1c6..cff09c5ec49c0675680b5f49c1f24f9008c93cc5 100644
--- a/cmake/modules/DuneCodegenMacros.cmake
+++ b/cmake/modules/DuneCodegenMacros.cmake
@@ -178,10 +178,8 @@ function(dune_add_generated_executable)
   endif()
 
   # Parse a mapping of operators to build and their respective filenames
-  dune_execute_process(COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env python ${dune-codegen_path}/deplist.py ${GEN_INIFILE} ${GEN_TARGET}
-                       OUTPUT_VARIABLE depdata
-                       )
-  parse_python_data(PREFIX depdata INPUT ${depdata})
+  dune_execute_process(COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env python ${dune-codegen_path}/deplist.py ${GEN_INIFILE} ${GEN_TARGET} ${CMAKE_CURRENT_BINARY_DIR}/interface.log)
+  parse_python_data(PREFIX depdata FILE ${CMAKE_CURRENT_BINARY_DIR}/interface.log)
 
   if(DUNE_CODEGEN_PROFILING)
     # This is a bit silly, but cProfile only finds entry point scripts
diff --git a/cmake/modules/GeneratedSystemtests.cmake b/cmake/modules/GeneratedSystemtests.cmake
index 87af7161c8b9c8196dbddc3ae42e55bcffd3bd81..2a1907950a7f5bd68e85a5c536465a73634828cf 100644
--- a/cmake/modules/GeneratedSystemtests.cmake
+++ b/cmake/modules/GeneratedSystemtests.cmake
@@ -119,9 +119,14 @@ function(dune_add_formcompiler_system_test)
   configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${SYSTEMTEST_INIFILE} ${CMAKE_CURRENT_BINARY_DIR}/tmp_${SYSTEMTEST_INIFILE})
 
   # expand the given meta ini file into the build tree
-  execute_process(COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env dune_expand_metaini.py --cmake --ini ${CMAKE_CURRENT_SOURCE_DIR}/${SYSTEMTEST_INIFILE} --dir ${CMAKE_CURRENT_BINARY_DIR} --section formcompiler
-                  OUTPUT_VARIABLE output)
-  parse_python_data(PREFIX INIINFO INPUT "${output}")
+  execute_process(COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env dune_expand_metaini.py
+                          --cmake
+                          --ini ${CMAKE_CURRENT_SOURCE_DIR}/${SYSTEMTEST_INIFILE}
+                          --dir ${CMAKE_CURRENT_BINARY_DIR}
+                          --section formcompiler
+                          --file ${CMAKE_CURRENT_BINARY_DIR}/interface.log
+                          )
+  parse_python_data(PREFIX INIINFO FILE ${CMAKE_CURRENT_BINARY_DIR}/interface.log)
 
   foreach(inifile ${INIINFO_names})
     if(${INIINFO_${inifile}_suffix} STREQUAL "__empty")
@@ -147,10 +152,11 @@ function(dune_add_formcompiler_system_test)
     # just the way that dune-testtools does.
     dune_execute_process(COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env dune_extract_static.py
                                --ini ${inifile}
+                               --file ${CMAKE_CURRENT_BINARY_DIR}/interface.log
                          WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
                          OUTPUT_VARIABLE output
                          ERROR_MESSAGE "Error extracting static info from ${inifile}")
-    parse_python_data(PREFIX STAT INPUT "${output}")
+    parse_python_data(PREFIX STAT FILE ${CMAKE_CURRENT_BINARY_DIR}/interface.log)
 
     foreach(config ${STAT___CONFIGS})
       foreach(cd ${STAT___STATIC_DATA})
diff --git a/cmake/modules/deplist.py b/cmake/modules/deplist.py
index 9cb5d7d42cbfc712e18e37798dffc2d553416f8c..7b0afac3c7c5f216e9710cef394d175c0c6c9ee9 100755
--- a/cmake/modules/deplist.py
+++ b/cmake/modules/deplist.py
@@ -22,5 +22,5 @@ def get_filename(operator):
 result = {"__{}".format(o): get_filename(o) for o in operators}
 result["__operators"] = ";".join(operators)
 
-printForCMake(result)
+printForCMake(result, sys.argv[3])
 sys.exit(0)
diff --git a/dune/codegen/sumfact/horizontaladd.hh b/dune/codegen/sumfact/horizontaladd.hh
index fc62dc47dab77330a125bb5f91805c89808480b0..7dd122af4c877b7206e99d71ba8aa05fc2fc448d 100644
--- a/dune/codegen/sumfact/horizontaladd.hh
+++ b/dune/codegen/sumfact/horizontaladd.hh
@@ -4,6 +4,50 @@
 #include<dune/codegen/common/vectorclass.hh>
 
 
+// Only use our custom implementations if we have AVX2 or later!
+#if INSTRSET >= 8
+
+/** Implement a variant of horizontal_add(Vec2d) that avoids the haddpd
+ *  instruction and instead uses the shuffle port.
+ */
+static inline double permuting_horizontal_add (const Vec2d & a)
+{
+    return _mm_cvtsd_f64(_mm_add_pd(_mm_permute_pd(a,1),a));
+}
+
+/** Implement a variant of horizontal_add(Vec4d) that avoids the haddpd
+ *  instruction and instead uses the shuffle port.
+ */
+static inline double permuting_horizontal_add (const Vec4d& a)
+{
+    __m128d valupper = _mm256_extractf128_pd(a, 1);
+    __m128d vallower = _mm256_castpd256_pd128(a);
+    __m128d valval = _mm_add_pd(valupper, vallower);
+    __m128d res = _mm_add_pd(_mm_permute_pd(valval,1), valval);
+    return _mm_cvtsd_f64(res);
+}
+
+#if MAX_VECTOR_SIZE >= 512
+
+/** Implement a variant of horizontal_add(Vec8d) that avoids the haddpd
+ *  instruction and instead uses the shuffle port.
+ */
+static inline double permuting_horizontal_add(const Vec8d& a)
+{
+  return permuting_horizontal_add(a.get_low() + a.get_high());
+}
+
+#endif
+
+#else
+template<typename V>
+static inline double permuting_horizontal_add (const V& a)
+{
+    return horizontal_add(a);
+}
+
+#endif
+
 template<class V>
 typename base_floatingpoint<V>::value horizontal_add_lower(const V& x)
 {
@@ -16,4 +60,16 @@ typename base_floatingpoint<V>::value horizontal_add_upper(const V& x)
   return horizontal_add(x.get_high());
 }
 
+template<class V>
+typename base_floatingpoint<V>::value permuting_horizontal_add_lower(const V& x)
+{
+  return permuting_horizontal_add(x.get_low());
+}
+
+template<class V>
+typename base_floatingpoint<V>::value permuting_horizontal_add_upper(const V& x)
+{
+  return permuting_horizontal_add(x.get_high());
+}
+
 #endif
diff --git a/dune/codegen/sumfact/oc_horizontaladd.hh b/dune/codegen/sumfact/oc_horizontaladd.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d136aaf5ffc226e8e2463c5c40b3d64c8a476fe7
--- /dev/null
+++ b/dune/codegen/sumfact/oc_horizontaladd.hh
@@ -0,0 +1,25 @@
+#ifndef DUNE_CODEGEN_SUMFACT_OCHORIZONTALADD_HH
+#define DUNE_CODEGEN_SUMFACT_OCHORIZONTALADD_HH
+
+#include<immintrin.h>
+#include<dune/codegen/common/simdtraits.hh>
+
+template<class V>
+typename base_floatingpoint<V>::value permuting_horizontal_add_lower(const V& x)
+{
+  return horizontal_add(x.get_low());
+}
+
+template<class V>
+typename base_floatingpoint<V>::value permuting_horizontal_add_upper(const V& x)
+{
+  return horizontal_add(x.get_high());
+}
+
+template<class V>
+typename base_floatingpoint<V>::value permuting_horizontal_add(const V& x)
+{
+  return horizontal_add(x);
+}
+
+#endif
diff --git a/patches/apply_patches.sh b/patches/apply_patches.sh
index 7d1d45112d392a13a3dec4881ba43222cf60054b..2d0cdc6f543f2cb43d3a28bf55f562db0e64f13a 100755
--- a/patches/apply_patches.sh
+++ b/patches/apply_patches.sh
@@ -5,10 +5,6 @@ git apply ../../patches/loopy/Current.patch
 git apply ../../patches/loopy/0001-Disable-a-logging-statement-that-breaks.patch
 popd
 
-pushd dune/codegen/vectorclass
-git apply ../../../patches/vectorclass/0001-Better-implementation-of-horizontal_add.patch
-popd
-
 pushd python/ufl
 git apply ../../patches/ufl/0001-Remove-special-case-for-variable-in-ufl2dot.patch
 popd
diff --git a/patches/vectorclass/0001-Alternative-implementation-of-horizontal_add-on-AVX5.patch b/patches/vectorclass/0001-Alternative-implementation-of-horizontal_add-on-AVX5.patch
deleted file mode 100644
index c5ca6dc30e2135ab30a28c7373b94da344b8a7ac..0000000000000000000000000000000000000000
--- a/patches/vectorclass/0001-Alternative-implementation-of-horizontal_add-on-AVX5.patch
+++ /dev/null
@@ -1,44 +0,0 @@
-From a324181d74fd8cd81fb945a4f66e4502ffbc68a0 Mon Sep 17 00:00:00 2001
-From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
-Date: Thu, 30 Nov 2017 18:51:49 +0100
-Subject: [PATCH] Alternative implementation of horizontal_add on AVX512
-
----
- vectorf512.h | 19 +++++++++++++------
- 1 file changed, 13 insertions(+), 6 deletions(-)
-
-diff --git a/vectorf512.h b/vectorf512.h
-index 0845d12..6a15ac2 100644
---- a/vectorf512.h
-+++ b/vectorf512.h
-@@ -1339,14 +1339,21 @@ static inline Vec8d if_mul (Vec8db const & f, Vec8d const & a, Vec8d const & b)
- 
- 
- // General arithmetic functions, etc.
-+#if __GNUC__ < 7
-+extern __inline double
-+__attribute__ ((__gnu_inline__, __always_inline__, __artificial__))
-+_mm512_cvtsd_f64 (__m512d __A)
-+{
-+  return __A[0];
-+}
-+#endif
- 
- // Horizontal add: Calculates the sum of all vector elements.
--static inline double horizontal_add (Vec8d const & a) {
--#if defined(__INTEL_COMPILER)
--    return _mm512_reduce_add_pd(a);
--#else
--    return horizontal_add(a.get_low() + a.get_high());
--#endif
-+static inline double horizontal_add (Vec8d const & x) {
-+    __m512d intermediate = _mm512_add_pd(x, _mm512_castsi512_pd(_mm512_alignr_epi64(_mm512_castpd_si512(x), _mm512_castpd_si512(x), 1)));
-+    intermediate = _mm512_add_pd(intermediate, _mm512_castsi512_pd(_mm512_alignr_epi64(_mm512_castpd_si512(intermediate), _mm512_castpd_si512(intermediate), 2)));
-+    intermediate = _mm512_add_pd(intermediate, _mm512_castsi512_pd(_mm512_alignr_epi64(_mm512_castpd_si512(intermediate), _mm512_castpd_si512(intermediate), 4)));
-+    return _mm512_cvtsd_f64(intermediate);
- }
- 
- // function max: a > b ? a : b
--- 
-2.1.4
-
diff --git a/patches/vectorclass/0001-Better-implementation-of-horizontal_add.patch b/patches/vectorclass/0001-Better-implementation-of-horizontal_add.patch
deleted file mode 100644
index fee83d7ad7cedbacc588c530ad9581b49cfa3b54..0000000000000000000000000000000000000000
--- a/patches/vectorclass/0001-Better-implementation-of-horizontal_add.patch
+++ /dev/null
@@ -1,32 +0,0 @@
-From 69f4ea4dcd018eb74c39a076a60fc27c0496e1dd Mon Sep 17 00:00:00 2001
-From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
-Date: Mon, 19 Jun 2017 13:07:22 +0200
-Subject: [PATCH] Better implementation of horizontal_add
-
----
- vectorf256.h | 9 +++++----
- 1 file changed, 5 insertions(+), 4 deletions(-)
-
-diff --git a/vectorf256.h b/vectorf256.h
-index db509f8..2bbd9de 100644
---- a/vectorf256.h
-+++ b/vectorf256.h
-@@ -1692,10 +1692,11 @@ static inline Vec4d if_mul (Vec4db const & f, Vec4d const & a, Vec4d const & b)
- 
- // Horizontal add: Calculates the sum of all vector elements.
- static inline double horizontal_add (Vec4d const & a) {
--    __m256d t1 = _mm256_hadd_pd(a,a);
--    __m128d t2 = _mm256_extractf128_pd(t1,1);
--    __m128d t3 = _mm_add_sd(_mm256_castpd256_pd128(t1),t2);
--    return _mm_cvtsd_f64(t3);        
-+    const __m128d valupper = _mm256_extractf128_pd(a, 1);
-+    const __m128d vallower = _mm256_castpd256_pd128(a);
-+    const __m128d valval = _mm_add_pd(valupper, vallower);
-+    const __m128d res = _mm_add_pd(_mm_permute_pd(valval,1), valval);
-+    return _mm_cvtsd_f64(res);
- }
- 
- // function max: a > b ? a : b
--- 
-2.1.4
-
diff --git a/python/dune/codegen/loopy/vcl.py b/python/dune/codegen/loopy/vcl.py
index e0943a69ac35136abfe59a7f641e0f076681643c..2431275a5f18a3bc87272711ab5ed71c038ced0b 100644
--- a/python/dune/codegen/loopy/vcl.py
+++ b/python/dune/codegen/loopy/vcl.py
@@ -119,9 +119,17 @@ def vcl_function_mangler(knl, func, arg_dtypes):
         return lp.CallMangleInfo("select", (vcl,), (vcl, vcl, vcl))
 
     if func in ("horizontal_add", "horizontal_add_lower", "horizontal_add_upper"):
+        if get_option("permuting_horizontal_add"):
+            func = "permuting_{}".format(func)
+
         dtype = arg_dtypes[0]
         vcl = lp.types.NumpyType(get_vcl_type(dtype))
-        include_file("dune/codegen/sumfact/horizontaladd.hh", filetag="operatorfile")
+
+        if get_option("opcounter"):
+            include_file("dune/codegen/sumfact/oc_horizontaladd.hh", filetag="operatorfile")
+        else:
+            include_file("dune/codegen/sumfact/horizontaladd.hh", filetag="operatorfile")
+
         return lp.CallMangleInfo(func, (lp.types.NumpyType(dtype.dtype),), (vcl,))
 
     if isinstance(func, VCLPermute):
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 72697492907cbf4afdfe62d1a3789606a8e2c290..9407709a515ef497b9f9b7561f570cf0411dd4ed 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -37,6 +37,7 @@ class CodegenGlobalOptionsArray(ImmutableRecord):
     debug_cache_with_stack = CodegenOption(default=False, helpstr="Store stack along with cache objects. Makes debugging caching issues easier.")
     driver_file = CodegenOption(helpstr="The filename for the generated driver header")
     explicit_time_stepping = CodegenOption(default=False, helpstr="use explicit time stepping")
+    time_stepping_order = CodegenOption(default=1, helpstr="Order of the time stepping method")
     exact_solution_expression = CodegenOption(helpstr="name of the exact solution expression in the ufl file")
     compare_l2errorsquared = CodegenOption(helpstr="maximal allowed l2 error squared of difference between numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)")
     grid_info = CodegenOption(default=None, helpstr="Path to file with information about facedir and facemod variations. This can be used to limit the generation of skeleton kernels.")
@@ -60,6 +61,7 @@ class CodegenGlobalOptionsArray(ImmutableRecord):
     use_sde = CodegenOption(default=False, helpstr="Use sde instead of own performance measurements.")
     autotune_google_benchmark = CodegenOption(default=False, helpstr="Use google-benchmark library for autotuning (when autotuning is activated).")
     with_mpi = CodegenOption(default=True, helpstr="The module was configured with mpi")
+    permuting_horizontal_add = CodegenOption(default=True, helpstr="Whether SIMD horizontal_add should use a permuting implementation.")
 
     # Arguments that are mainly to be set by logic depending on other options
     max_vector_width = CodegenOption(default=256, helpstr=None)
diff --git a/python/dune/codegen/pdelab/driver/gridoperator.py b/python/dune/codegen/pdelab/driver/gridoperator.py
index abc878f0715dbcb4a0501d56247c4f03c98f20df..bdf8b2f6245acb63b2a67c518ab98fef597672af 100644
--- a/python/dune/codegen/pdelab/driver/gridoperator.py
+++ b/python/dune/codegen/pdelab/driver/gridoperator.py
@@ -77,8 +77,7 @@ def typedef_localoperator(name, form_ident):
     filename = get_form_option("filename", form_ident)
     include_file(filename, filetag="driver")
     lopname = localoperator_basename(form_ident)
-    range_type = type_range()
-    return "using {} = {}<{}, {}, {}>;".format(name, lopname, ugfs, vgfs, range_type)
+    return "using {} = {}<{}, {}>;".format(name, lopname, ugfs, vgfs)
 
 
 def type_localoperator(form_ident):
diff --git a/python/dune/codegen/pdelab/driver/instationary.py b/python/dune/codegen/pdelab/driver/instationary.py
index 5cf7170f1dca6f4274db4bf77b4b1627466e20cb..355fd0e743f6a59716319241018d3c5794c784bd 100644
--- a/python/dune/codegen/pdelab/driver/instationary.py
+++ b/python/dune/codegen/pdelab/driver/instationary.py
@@ -132,10 +132,25 @@ def name_time():
 def typedef_timesteppingmethod(name):
     r_type = type_range()
     explicit = get_option('explicit_time_stepping')
+    order = get_option('time_stepping_order')
     if explicit:
-        return "using {} = Dune::PDELab::ExplicitEulerParameter<{}>;".format(name, r_type)
+        if order == 1:
+            return "using {} = Dune::PDELab::ExplicitEulerParameter<{}>;".format(name, r_type)
+        elif order == 2:
+            return "using {} = Dune::PDELab::HeunParameter<{}>;".format(name, r_type)
+        elif order == 3:
+            return "using {} = Dune::PDELab::Shu3Parameter<{}>;".format(name, r_type)
+        elif order == 4:
+            return "using {} = Dune::PDELab::RK4Parameter<{}>;".format(name, r_type)
+        else:
+            raise NotImplementedError("Time stepping method not supported")
     else:
-        return "using {} = Dune::PDELab::OneStepThetaParameter<{}>;".format(name, r_type)
+        if order == 1:
+            return "using {} = Dune::PDELab::OneStepThetaParameter<{}>;".format(name, r_type)
+        elif order == 2:
+            return "using {} = Dune::PDELab::Alexander2Parameter<{}>;".format(name, r_type)
+        elif order == 3:
+            return "using {} = Dune::PDELab::Alexander3Parameter<{}>;".format(name, r_type)
 
 
 def type_timesteppingmethod():
@@ -150,8 +165,12 @@ def define_timesteppingmethod(name):
     if explicit:
         return "{} {};".format(tsm_type, name)
     else:
-        ini = name_initree()
-        return "{} {}({}.get<double>(\"instat.theta\",1.0));".format(tsm_type, name, ini)
+        order = get_option('time_stepping_order')
+        if order == 1:
+            ini = name_initree()
+            return "{} {}({}.get<double>(\"instat.theta\",1.0));".format(tsm_type, name, ini)
+        else:
+            return "{} {};".format(tsm_type, name)
 
 
 def name_timesteppingmethod():
diff --git a/python/dune/codegen/pdelab/geometry.py b/python/dune/codegen/pdelab/geometry.py
index b9d5b0d03e9c137c30fa49bfefee18d835e7f003..da5c7b395f151cb6e3f25c734243a6732c2a95fc 100644
--- a/python/dune/codegen/pdelab/geometry.py
+++ b/python/dune/codegen/pdelab/geometry.py
@@ -269,9 +269,9 @@ class EquidistantGeometryMixin(AxiparallelGeometryMixin):
 
     @preamble(kernel="operator")
     def _define_jacobian_determinant_eval(self, name):
-        from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
+        from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param
         gfs = name_ansatz_gfs_constructor_param()
-        rft = lop_template_range_field()
+        rft = type_floatingpoint()
         return "{} = {}.gridView().template begin<0>()->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
 
     def define_jacobian_inverse_transposed(self, name, restriction):
@@ -289,9 +289,9 @@ class EquidistantGeometryMixin(AxiparallelGeometryMixin):
 
     @preamble(kernel="operator")
     def _define_jacobian_inverse_transposed_eval(self, name):
-        from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
+        from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param
         gfs = name_ansatz_gfs_constructor_param()
-        rft = lop_template_range_field()
+        rft = type_floatingpoint()
         return "{} = {}.gridView().template begin<0>()->geometry().jacobianInverseTransposed(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
 
 
diff --git a/python/dune/codegen/pdelab/localoperator.py b/python/dune/codegen/pdelab/localoperator.py
index 0ca6bfb45ed80271af92800980df7405200dca8b..171f9af707fb744bddf97d9faae9cdf52875043a 100644
--- a/python/dune/codegen/pdelab/localoperator.py
+++ b/python/dune/codegen/pdelab/localoperator.py
@@ -38,6 +38,7 @@ from dune.codegen.cgen.clazz import (AccessModifier,
                                      BaseClass,
                                      ClassMember,
                                      )
+from dune.codegen.loopy.target import type_floatingpoint
 from dune.codegen.ufl.modified_terminals import Restriction
 
 import dune.codegen.loopy.mangler
@@ -73,11 +74,6 @@ def name_test_gfs_constructor_param():
     return "gfsv"
 
 
-@template_parameter(classtag="operator")
-def lop_template_range_field():
-    return "RF"
-
-
 @class_member(classtag="operator")
 def lop_domain_field(name):
     # TODO: Rethink for not Galerkin Method
@@ -501,7 +497,8 @@ def visit_integral(integral):
 
     # Start the visiting process!
     visitor = get_visitor(measure, subdomain_id)
-    visitor.accumulate(integrand)
+    with global_context(visitor=visitor):
+        visitor.accumulate(integrand)
 
     run_hook(name="after_visit",
              args=(visitor,),
@@ -857,7 +854,6 @@ def local_operator_default_settings(operator, form):
     localoperator_basename(operator)
     lop_template_ansatz_gfs()
     lop_template_test_gfs()
-    lop_template_range_field()
 
     # Make sure there is always the same constructor arguments, even if some of them are
     # not strictly needed. Also ensure the order.
@@ -871,7 +867,7 @@ def local_operator_default_settings(operator, form):
     base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")
     from dune.codegen.pdelab.driver import is_stationary
     if not is_stationary():
-        rf = lop_template_range_field()
+        rf = type_floatingpoint()
         base_class('Dune::PDELab::InstationaryLocalOperatorDefaultMethods<{}>'
                    .format(rf), classtag="operator")
 
diff --git a/python/dune/codegen/pdelab/quadrature.py b/python/dune/codegen/pdelab/quadrature.py
index ae6a7e2db7212b254e3857bc86813ea420ed0251..ccfb9f9c1a42888923eb07b753a013705190df01 100644
--- a/python/dune/codegen/pdelab/quadrature.py
+++ b/python/dune/codegen/pdelab/quadrature.py
@@ -12,7 +12,8 @@ from dune.codegen.generation import (class_member,
                                      temporary_variable,
                                      valuearg,
                                      )
-from dune.codegen.pdelab.localoperator import lop_template_range_field, name_ansatz_gfs_constructor_param
+from dune.codegen.loopy.target import type_floatingpoint
+from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param
 from dune.codegen.options import get_form_option
 
 from pymbolic.primitives import Variable, Subscript
@@ -51,7 +52,7 @@ class GenericQuadratureMixin(QuadratureMixinBase):
 
     @class_member(classtag="operator")
     def define_quadrature_weights(self, name):
-        rf = lop_template_range_field()
+        rf = type_floatingpoint()
         from dune.codegen.pdelab.geometry import local_dimension
         dim = local_dimension()
         self.eval_quadrature_weights(name)
@@ -89,7 +90,7 @@ class GenericQuadratureMixin(QuadratureMixinBase):
 
     @class_member(classtag="operator")
     def define_quadrature_points(self, name):
-        rf = lop_template_range_field()
+        rf = type_floatingpoint()
         from dune.codegen.pdelab.geometry import local_dimension
         dim = local_dimension()
         self.eval_quadrature_points(name)
diff --git a/python/dune/codegen/sumfact/accumulation.py b/python/dune/codegen/sumfact/accumulation.py
index 90fd01fd5b90b722c0a3ee194e1da82ed44f6eab..25be6f7eb80b365d8eead39c3cbb93c0e36c52bc 100644
--- a/python/dune/codegen/sumfact/accumulation.py
+++ b/python/dune/codegen/sumfact/accumulation.py
@@ -39,9 +39,6 @@ from dune.codegen.sumfact.permutation import (permute_backward,
 from dune.codegen.sumfact.tabulation import (basis_functions_per_direction,
                                              construct_basis_matrix_sequence,
                                              )
-from dune.codegen.sumfact.switch import (get_facedir,
-                                         get_facemod,
-                                         )
 from dune.codegen.sumfact.symbolic import SumfactKernel, SumfactKernelInterfaceBase
 from dune.codegen.ufl.modified_terminals import extract_modified_arguments
 from dune.codegen.tools import get_pymbolic_basename, get_leaf, ImmutableCuttingRecord
@@ -361,6 +358,22 @@ class SumfactAccumulationMixin(AccumulationMixinBase):
     def generate_accumulation_instruction(self, expr):
         return generate_accumulation_instruction(expr, self)
 
+    def get_facedir(self, restriction):
+        from dune.codegen.pdelab.restriction import Restriction
+        if restriction == Restriction.POSITIVE or get_global_context_value("integral_type") == "exterior_facet":
+            return get_global_context_value("facedir_s")
+        if restriction == Restriction.NEGATIVE:
+            return get_global_context_value("facedir_n")
+        return None
+
+    def get_facemod(self, restriction):
+        from dune.codegen.pdelab.restriction import Restriction
+        if restriction == Restriction.POSITIVE or get_global_context_value("integral_type") == "exterior_facet":
+            return get_global_context_value("facemod_s")
+        if restriction == Restriction.NEGATIVE:
+            return get_global_context_value("facemod_n")
+        return None
+
     def additional_matrix_sequence(self):
         return None
 
@@ -375,8 +388,8 @@ class SumfactPointDiagonalAccumulationMixin(SumfactAccumulationMixin):
         info = self.current_info[1]
         return construct_basis_matrix_sequence(transpose=True,
                                                derivative=info.grad_index,
-                                               facedir=get_facedir(info.restriction),
-                                               facemod=get_facemod(info.restriction),
+                                               facedir=self.get_facedir(info.restriction),
+                                               facemod=self.get_facemod(info.restriction),
                                                basis_size=get_basis_size(info),
                                                )
 
@@ -582,8 +595,8 @@ def generate_accumulation_instruction(expr, visitor):
     matrix_sequence = construct_basis_matrix_sequence(
         transpose=True,
         derivative=test_info.grad_index,
-        facedir=get_facedir(test_info.restriction),
-        facemod=get_facemod(test_info.restriction),
+        facedir=visitor.get_facedir(test_info.restriction),
+        facemod=visitor.get_facemod(test_info.restriction),
         basis_size=basis_size,
         additional_sequence=visitor.additional_matrix_sequence())
 
diff --git a/python/dune/codegen/sumfact/basis.py b/python/dune/codegen/sumfact/basis.py
index c49e052a90faac34ad6f64b539f93b3bb1d560cb..83e5206b5289719f3c4e85f4ca523976cb6feb79 100644
--- a/python/dune/codegen/sumfact/basis.py
+++ b/python/dune/codegen/sumfact/basis.py
@@ -31,9 +31,6 @@ from dune.codegen.sumfact.permutation import (permute_backward,
                                               sumfact_cost_permutation_strategy,
                                               sumfact_quadrature_permutation_strategy,
                                               )
-from dune.codegen.sumfact.switch import (get_facedir,
-                                         get_facemod,
-                                         )
 from dune.codegen.pdelab.argument import name_coefficientcontainer, name_applycontainer
 from dune.codegen.pdelab.basis import GenericBasisMixin
 from dune.codegen.pdelab.geometry import (local_dimension,
@@ -86,7 +83,7 @@ class SumfactBasisMixin(GenericBasisMixin):
         temporary_variable(name, shape=())
         quad_inames = self.quadrature_inames()
         inames = self.lfs_inames(element, restriction)
-        facedir = get_facedir(restriction)
+        facedir = self.get_facedir(restriction)
 
         # Collect the pairs of lfs/quad inames that are in use
         # On facets, the normal direction of the facet is excluded
@@ -106,7 +103,7 @@ class SumfactBasisMixin(GenericBasisMixin):
 
         # Add the missing direction on facedirs by evaluating at either 0 or 1
         if facedir is not None:
-            facemod = get_facemod(restriction)
+            facemod = self.get_facemod(restriction)
             prod = prod + (prim.Call(PolynomialLookup(name_polynomials(element.degree()), False),
                                      (prim.Variable(inames[facedir]), facemod)),)
 
@@ -141,7 +138,7 @@ class SumfactBasisMixin(GenericBasisMixin):
         temporary_variable(name, shape=())
         quad_inames = self.quadrature_inames()
         inames = self.lfs_inames(element, restriction)
-        facedir = get_facedir(restriction)
+        facedir = self.get_facedir(restriction)
 
         # Map the direction to a quadrature iname
         quadinamemapping = {}
@@ -161,7 +158,7 @@ class SumfactBasisMixin(GenericBasisMixin):
                 prod.append(tab.pymbolic((prim.Variable(quadinamemapping[i]), prim.Variable(inames[i]))))
 
         if facedir is not None:
-            facemod = get_facemod(restriction)
+            facemod = self.get_facemod(restriction)
             prod.append(prim.Call(PolynomialLookup(name_polynomials(element.degree()), index == facedir),
                                   (prim.Variable(inames[facedir]), facemod)),)
 
@@ -197,8 +194,8 @@ class SumfactBasisMixin(GenericBasisMixin):
 
         # Construct the matrix sequence for this sum factorization
         matrix_sequence = construct_basis_matrix_sequence(derivative=derivative,
-                                                          facedir=get_facedir(restriction),
-                                                          facemod=get_facemod(restriction),
+                                                          facedir=self.get_facedir(restriction),
+                                                          facemod=self.get_facemod(restriction),
                                                           basis_size=basis_size)
 
         inp = LFSSumfactKernelInput(matrix_sequence=matrix_sequence,
diff --git a/python/dune/codegen/sumfact/geometry.py b/python/dune/codegen/sumfact/geometry.py
index 8ac8aaa4faa1a4849c3b88fb54f87b0f87ca7b2b..811a201e9e11238bcdc4965ce2c7f75426fc3487 100644
--- a/python/dune/codegen/sumfact/geometry.py
+++ b/python/dune/codegen/sumfact/geometry.py
@@ -16,6 +16,7 @@ from dune.codegen.generation import (class_member,
                                      valuearg,
                                      )
 from dune.codegen.loopy.flatten import flatten_index
+from dune.codegen.loopy.target import type_floatingpoint
 from dune.codegen.options import get_option
 from dune.codegen.pdelab.geometry import (enforce_boundary_restriction,
                                           local_dimension,
@@ -28,7 +29,6 @@ from dune.codegen.pdelab.geometry import (enforce_boundary_restriction,
                                           )
 from dune.codegen.pdelab.localoperator import (name_ansatz_gfs_constructor_param,
                                                lop_template_ansatz_gfs,
-                                               lop_template_range_field,
                                                )
 from dune.codegen.pdelab.restriction import restricted_name
 from dune.codegen.sumfact.accumulation import basis_sf_kernels, sumfact_iname
@@ -39,7 +39,6 @@ from dune.codegen.sumfact.permutation import (permute_backward,
                                               sumfact_quadrature_permutation_strategy,
                                               )
 from dune.codegen.sumfact.quadrature import additional_inames
-from dune.codegen.sumfact.switch import get_facedir, get_facemod
 from dune.codegen.sumfact.symbolic import SumfactKernelInterfaceBase, SumfactKernel
 from dune.codegen.tools import get_pymbolic_basename, ImmutableCuttingRecord
 from dune.codegen.options import get_form_option, option_switch
@@ -162,8 +161,8 @@ class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
     def outer_normal(self):
         """ This is the *unnormalized* outer normal """
         name = "outer_normal"
-        facedir_s = get_facedir(Restriction.POSITIVE)
-        facemod_s = get_facemod(Restriction.POSITIVE)
+        facedir_s = self.get_facedir(Restriction.POSITIVE)
+        facemod_s = self.get_facemod(Restriction.POSITIVE)
 
         temporary_variable(name, shape=(world_dimension(),))
         for i in range(world_dimension()):
@@ -210,8 +209,8 @@ class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
         restriction = enforce_boundary_restriction(self)
 
         # Generate sum factorization kernel and add vectorization info
-        matrix_sequence = construct_basis_matrix_sequence(facedir=get_facedir(restriction),
-                                                          facemod=get_facemod(restriction),
+        matrix_sequence = construct_basis_matrix_sequence(facedir=self.get_facedir(restriction),
+                                                          facemod=self.get_facemod(restriction),
                                                           basis_size=(2,) * world_dimension())
         inp = GeoCornersInput(matrix_sequence=matrix_sequence,
                               direction=self.indices[0],
@@ -253,8 +252,8 @@ class SumfactAxiParallelGeometryMixin(AxiparallelGeometryMixin):
         assert isinstance(i, int)
 
         # Use facemod_s and facedir_s
-        if i == get_facedir(Restriction.POSITIVE):
-            if get_facemod(Restriction.POSITIVE):
+        if i == self.get_facedir(Restriction.POSITIVE):
+            if self.get_facemod(Restriction.POSITIVE):
                 return 1
             else:
                 return -1
@@ -270,7 +269,7 @@ class SumfactEqudistantGeometryMixin(EquidistantGeometryMixin, SumfactAxiParalle
     def facet_jacobian_determinant(self, o):
         name = "fdetjac"
         self.define_facet_jacobian_determinant(name)
-        facedir = get_facedir(Restriction.POSITIVE)
+        facedir = self.get_facedir(Restriction.POSITIVE)
         globalarg(name, shape=(world_dimension(),))
         return prim.Subscript(prim.Variable(name), (facedir,))
 
@@ -283,7 +282,7 @@ class SumfactEqudistantGeometryMixin(EquidistantGeometryMixin, SumfactAxiParalle
     @kernel_cached(kernel="operator")
     def _define_facet_jacobian_determinant_eval(self, name):
         gfs = name_ansatz_gfs_constructor_param()
-        rft = lop_template_range_field()
+        rft = type_floatingpoint()
         code = ["{",
                 "  auto e = *({}.gridView().template begin<0>());".format(gfs),
                 "  int dir=0;",
@@ -304,8 +303,7 @@ class SumfactEqudistantGeometryMixin(EquidistantGeometryMixin, SumfactAxiParalle
         restriction = Restriction.NONE
         if self.measure == "interior_facet":
             restriction = Restriction.POSITIVE
-        from dune.codegen.sumfact.switch import get_facedir
-        face = get_facedir(restriction)
+        face = self.get_facedir(restriction)
 
         lowcorner = name_lowerleft_corner()
         meshwidth = name_meshwidth()
@@ -483,8 +481,7 @@ def define_corner(name, low):
 
 @class_member(classtag="operator")
 def define_mesh_width(name):
-    from dune.codegen.pdelab.localoperator import lop_template_range_field
-    rft = lop_template_range_field()
+    rft = type_floatingpoint()
     define_mesh_width_eval(name)
     return "Dune::FieldVector<{}, {}> {};".format(rft, world_dimension(), name)
 
@@ -528,8 +525,8 @@ def _name_jacobian(i, j, restriction, visitor):
     """
     # Create matrix sequence with derivative in j direction
     matrix_sequence = construct_basis_matrix_sequence(derivative=j,
-                                                      facedir=get_facedir(restriction),
-                                                      facemod=get_facemod(restriction),
+                                                      facedir=visitor.get_facedir(restriction),
+                                                      facemod=visitor.get_facemod(restriction),
                                                       basis_size=(2,) * world_dimension())
 
     # Sum factorization input for the i'th component of the geometry mapping
diff --git a/python/dune/codegen/sumfact/permutation.py b/python/dune/codegen/sumfact/permutation.py
index 7f37dfeae795031ec81b292b96e65649f8fc78cc..6e9fdaaad5dc9021d0ec076df7cc743806eab512 100644
--- a/python/dune/codegen/sumfact/permutation.py
+++ b/python/dune/codegen/sumfact/permutation.py
@@ -3,7 +3,6 @@
 import itertools
 
 from dune.codegen.options import get_option
-from dune.codegen.sumfact.switch import get_facedir, get_facemod
 from dune.codegen.ufl.modified_terminals import Restriction
 
 
@@ -125,8 +124,9 @@ def sumfact_quadrature_permutation_strategy(dim, restriction):
         # all others can be derived by rotating the cube and matching edge
         # directions.
         def _order_on_self(restriction):
-            facedir = get_facedir(restriction)
-            facemod = get_facemod(restriction)
+            from dune.codegen.sumfact.accumulation import SumfactAccumulationMixin
+            facedir = SumfactAccumulationMixin.get_facedir(None, restriction)
+            facemod = SumfactAccumulationMixin.get_facemod(None, restriction)
 
             quadrature_order = {
                 (0, 0): (0, 1, 2),
diff --git a/python/dune/codegen/sumfact/quadrature.py b/python/dune/codegen/sumfact/quadrature.py
index 91e99c4c3c0cd333d6e64353ed06ac0c129a7e05..ae387426ebda79a8435811c548e91761f2de1b3c 100644
--- a/python/dune/codegen/sumfact/quadrature.py
+++ b/python/dune/codegen/sumfact/quadrature.py
@@ -10,7 +10,6 @@ from dune.codegen.generation import (domain,
                                      quadrature_mixin,
                                      temporary_variable,
                                      )
-from dune.codegen.sumfact.switch import get_facedir
 from dune.codegen.sumfact.tabulation import (quadrature_points_per_direction,
                                              local_quadrature_points_per_direction,
                                              name_oned_quadrature_points,
@@ -22,7 +21,6 @@ from dune.codegen.pdelab.geometry import (local_dimension,
                                           )
 from dune.codegen.pdelab.quadrature import GenericQuadratureMixin
 from dune.codegen.options import get_form_option
-from dune.codegen.sumfact.switch import get_facedir
 from dune.codegen.loopy.target import dtype_floatingpoint
 
 from loopy import CallMangleInfo
diff --git a/python/dune/codegen/sumfact/switch.py b/python/dune/codegen/sumfact/switch.py
index a031790dd5321841eb8c100195cde675651fe1d3..65d35c17e4a8d2753e680c8c440fee6e9f4c6adf 100644
--- a/python/dune/codegen/sumfact/switch.py
+++ b/python/dune/codegen/sumfact/switch.py
@@ -169,25 +169,3 @@ def generate_interior_facet_switch():
     block.append("}")
 
     return ClassMember(signature + block)
-
-
-def get_facedir(restriction):
-    from dune.codegen.pdelab.restriction import Restriction
-    if restriction == Restriction.POSITIVE or get_global_context_value("integral_type") == "exterior_facet":
-        return get_global_context_value("facedir_s")
-    if restriction == Restriction.NEGATIVE:
-        return get_global_context_value("facedir_n")
-    if restriction == Restriction.NONE:
-        return None
-    assert False
-
-
-def get_facemod(restriction):
-    from dune.codegen.pdelab.restriction import Restriction
-    if restriction == Restriction.POSITIVE or get_global_context_value("integral_type") == "exterior_facet":
-        return get_global_context_value("facemod_s")
-    if restriction == Restriction.NEGATIVE:
-        return get_global_context_value("facemod_n")
-    if restriction == Restriction.NONE:
-        return None
-    assert False
diff --git a/python/dune/codegen/sumfact/tabulation.py b/python/dune/codegen/sumfact/tabulation.py
index e2e4a771ec9918f960e5ea88821aa7acb885a0c2..0170c59463eed570fcce679636899e0d633d32fe 100644
--- a/python/dune/codegen/sumfact/tabulation.py
+++ b/python/dune/codegen/sumfact/tabulation.py
@@ -18,12 +18,10 @@ from dune.codegen.generation import (class_member,
                                      transform,
                                      valuearg
                                      )
-from dune.codegen.loopy.target import dtype_floatingpoint
+from dune.codegen.loopy.target import dtype_floatingpoint, type_floatingpoint
 from dune.codegen.loopy.vcl import ExplicitVCLCast, get_vcl_type_size
 from dune.codegen.options import get_option
-from dune.codegen.pdelab.localoperator import (name_domain_field,
-                                               lop_template_range_field,
-                                               )
+from dune.codegen.pdelab.localoperator import name_domain_field
 from dune.codegen.pdelab.quadrature import quadrature_order
 from dune.codegen.tools import maybe_wrap_subscript, ceildiv
 from loopy import CallMangleInfo
@@ -345,7 +343,7 @@ def name_oned_quadrature_points(bound):
 
 @class_member(classtag="operator")
 def typedef_polynomials(name, degree):
-    range_field = lop_template_range_field()
+    range_field = type_floatingpoint()
     domain_field = name_domain_field()
 
     include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="operatorfile")
@@ -380,7 +378,7 @@ def name_polynomials(degree):
 
 
 def sort_quadrature_points_weights(qp, qw, bound):
-    range_field = lop_template_range_field()
+    range_field = type_floatingpoint()
     domain_field = name_domain_field()
     include_file("dune/codegen/sumfact/onedquadrature.hh", filetag="operatorfile")
     return frozenset({instruction(code="onedQuadraturePointsWeights<{}, {}, {}>({}, {});"
diff --git a/test/adjoint/poisson_mc_driver.hh b/test/adjoint/poisson_mc_driver.hh
index 5bd222b800baf36c09bd172c48d706a135912c80..a3ac578afc8d9a2039dd9437ef48fe971a37fcb6 100644
--- a/test/adjoint/poisson_mc_driver.hh
+++ b/test/adjoint/poisson_mc_driver.hh
@@ -60,7 +60,7 @@ bool driver(int argc, char** argv){
   Dune::PDELab::constraints(p1_bctype, p1_dirichlet_gfs_, p1_dirichlet_gfs__cc);
 
   // Set up grid grid operators...
-  using LOP_R = ROperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType>;
+  using LOP_R = ROperator<P1_dirichlet_GFS, P1_dirichlet_GFS>;
   using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
   using GO_r = Dune::PDELab::GridOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, P1_dirichlet_GFS_CC, P1_dirichlet_GFS_CC>;
   LOP_R lop_r(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree);
@@ -122,7 +122,7 @@ bool driver(int argc, char** argv){
   GF_X x_gf(p1_dirichlet_gfs_, x_r);
 
   // Local operator for adjoint problem
-  using LOP_Adjoint = RAdjointOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType, GF_X>;
+  using LOP_Adjoint = RAdjointOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, GF_X>;
   LOP_Adjoint lop_adjoint(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree, x_gf);
 
   // Grid operator for adjoint problem
@@ -156,7 +156,7 @@ bool driver(int argc, char** argv){
   DJDM dJdm(7,0.0);
 
   // Local operator for control problem
-  using LOP_Control = RControlOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType, GF_Adjoint, DJDM>;
+  using LOP_Control = RControlOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, GF_Adjoint, DJDM>;
   LOP_Control lop_control(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree, gf_adjoint, dJdm);
 
   // Grid operator for control problem
diff --git a/test/coeffeval/coeffeval_poisson.cc b/test/coeffeval/coeffeval_poisson.cc
index 4f9e8861b5139dacb8183327ad053b3cb2373b8b..f0d645ffb5a4521f80c82684de736412d49f737e 100644
--- a/test/coeffeval/coeffeval_poisson.cc
+++ b/test/coeffeval/coeffeval_poisson.cc
@@ -74,7 +74,7 @@ int main(int argc, char** argv)
   GF c_gf(p2_gfs, c);
 
   // Local Operator
-  using LOP_R = PoissonLocalOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType, GF>;
+  using LOP_R = PoissonLocalOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, GF>;
   LOP_R lop_r(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree, c_gf);
 
   // Constraints stuff
diff --git a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation_driver.cc
index 2f8cc6c8614e0cdb954a52fa5de808dbff466b07..d9fa57bf0d692186a5761667508f9e7c1b645715 100644
--- a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation_driver.cc
+++ b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation_driver.cc
@@ -109,7 +109,7 @@ int main(int argc, char** argv){
     Dune::PDELab::constraints(dg1_gfs_, dg1_gfs__cc);
 
     // Set up grid grid operators...
-    using LOP_R = rOperator<DG1_GFS, DG1_GFS, RangeType>;
+    using LOP_R = rOperator<DG1_GFS, DG1_GFS>;
     using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
     using GO_r = Dune::PDELab::GridOperator<DG1_GFS, DG1_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, DG1_GFS_CC, DG1_GFS_CC>;
     LOP_R lop_r(dg1_gfs_, dg1_gfs_, initree);
diff --git a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc
index 4b57deaa50bc6e435b40dec8060d55e927797b2a..cc4c1d01d3d6d5bb67f123e98696e7945002dc45 100644
--- a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc
+++ b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc
@@ -251,7 +251,7 @@ int main(int argc, char** argv){
     Dune::PDELab::constraints(dg1_gfs_, dg1_gfs__cc);
 
     // Set up grid grid operators...
-    using LOP_R = CLASSNAME<DG1_GFS, DG1_GFS, RangeType>;
+    using LOP_R = CLASSNAME<DG1_GFS, DG1_GFS>;
     using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
     using GO_r = Dune::PDELab::GridOperator<DG1_GFS, DG1_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, DG1_GFS_CC, DG1_GFS_CC>;
     LOP_R lop_r(dg1_gfs_, dg1_gfs_, initree);
diff --git a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation_driver.cc
index 90a5c381b538f8c41403161f96ba78a960fe491b..ad75cfc4ed894a45879fd2e078f01fa53aff3b09 100644
--- a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation_driver.cc
+++ b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation_driver.cc
@@ -225,7 +225,7 @@ int main(int argc, char** argv){
     Dune::PDELab::constraints(dg1_gfs_, dg1_gfs__cc);
 
     // Set up grid grid operators...
-    using LOP_R = CLASSNAME<DG1_GFS, DG1_GFS, RangeType>;
+    using LOP_R = CLASSNAME<DG1_GFS, DG1_GFS>;
     using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
     using GO_r = Dune::PDELab::FastDGGridOperator<DG1_GFS, DG1_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, DG1_GFS_CC, DG1_GFS_CC>;
     LOP_R lop_r(dg1_gfs_, dg1_gfs_, initree);
diff --git a/test/sumfact/preconditioner/test_preconditioner_2d.cc b/test/sumfact/preconditioner/test_preconditioner_2d.cc
index 2781b32399f800d1d6334565aefbede360fdbfa9..3e8fcfc117288feba4d5bd76abb5ca548e18f649 100644
--- a/test/sumfact/preconditioner/test_preconditioner_2d.cc
+++ b/test/sumfact/preconditioner/test_preconditioner_2d.cc
@@ -56,7 +56,7 @@ int main(int argc, char** argv){
     Dune::PDELab::constraints(dg2_gfs_, dg2_gfs__cc);
     
     // Set up grid grid operators...    
-    using FullLOP = FullOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using FullLOP = FullOperator<DG2_GFS, DG2_GFS>;
     using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
     using FullGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, FullLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     FullLOP fulllop(dg2_gfs_, dg2_gfs_, initree);
@@ -65,17 +65,17 @@ int main(int argc, char** argv){
     FullGO fullgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, fulllop, mb);
 
     // Additional grid operators for preconditioner
-    using BDLOP = BlockDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BDLOP = BlockDiagonalOperator<DG2_GFS, DG2_GFS>;
     using BDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     BDLOP bdlop(dg2_gfs_, dg2_gfs_, initree);
     BDGO bdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bdlop, mb);
 
-    using BODLOP = BlockOffDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BODLOP = BlockOffDiagonalOperator<DG2_GFS, DG2_GFS>;
     using BODGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BODLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     BODLOP bodlop(dg2_gfs_, dg2_gfs_, initree);
     BODGO bodgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bodlop, mb);
 
-    using PDLOP = PointDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using PDLOP = PointDiagonalOperator<DG2_GFS, DG2_GFS>;
     using PDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, PDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     PDLOP pdlop(dg2_gfs_, dg2_gfs_, initree);
     PDGO pdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, pdlop, mb);
diff --git a/test/sumfact/preconditioner/test_preconditioner_3d.cc b/test/sumfact/preconditioner/test_preconditioner_3d.cc
index ec84a1efbf7e19f46d8531622374f2902a3139eb..90ea7aee7def84e9fc2395a798ace1bc39018f7d 100644
--- a/test/sumfact/preconditioner/test_preconditioner_3d.cc
+++ b/test/sumfact/preconditioner/test_preconditioner_3d.cc
@@ -56,7 +56,7 @@ int main(int argc, char** argv){
     Dune::PDELab::constraints(dg2_gfs_, dg2_gfs__cc);
     
     // Set up grid grid operators...    
-    using FullLOP = FullOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using FullLOP = FullOperator<DG2_GFS, DG2_GFS>;
     using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
     using FullGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, FullLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     FullLOP fulllop(dg2_gfs_, dg2_gfs_, initree);
@@ -65,17 +65,17 @@ int main(int argc, char** argv){
     FullGO fullgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, fulllop, mb);
 
     // Additional grid operators for preconditioner
-    using BDLOP = BlockDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BDLOP = BlockDiagonalOperator<DG2_GFS, DG2_GFS>;
     using BDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     BDLOP bdlop(dg2_gfs_, dg2_gfs_, initree);
     BDGO bdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bdlop, mb);
 
-    using BODLOP = BlockOffDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using BODLOP = BlockOffDiagonalOperator<DG2_GFS, DG2_GFS>;
     using BODGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, BODLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     BODLOP bodlop(dg2_gfs_, dg2_gfs_, initree);
     BODGO bodgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, bodlop, mb);
 
-    using PDLOP = PointDiagonalOperator<DG2_GFS, DG2_GFS, RangeType>;
+    using PDLOP = PointDiagonalOperator<DG2_GFS, DG2_GFS>;
     using PDGO = Dune::PDELab::FastDGGridOperator<DG2_GFS, DG2_GFS, PDLOP, MatrixBackend, DF, RangeType, RangeType, DG2_GFS_CC, DG2_GFS_CC>;
     PDLOP pdlop(dg2_gfs_, dg2_gfs_, initree);
     PDGO pdgo(dg2_gfs_, dg2_gfs__cc, dg2_gfs_, dg2_gfs__cc, pdlop, mb);