localoperator.py

from __future__ import absolute_import
from os.path import splitext

from dune.perftool.options import get_option
from dune.perftool.generation import (backend,
                                      base_class,
                                      class_basename,
                                      class_member,
                                      constructor_parameter,
                                      domain,
                                      dump_accumulate_timer,
                                      get_backend,
                                      get_global_context_value,
                                      global_context,
                                      iname,
                                      include_file,
                                      initializer_list,
                                      post_include,
                                      retrieve_cache_items,
                                      template_parameter,
                                      )
from dune.perftool.cgen.clazz import (AccessModifier,
                                      BaseClass,
                                      ClassMember,
                                      )
from dune.perftool.ufl.modified_terminals import Restriction
from pymbolic.primitives import Variable
from pytools import Record

import loopy as lp


def name_form(formdata, data):
    # Check wether the formdata has a name in UFL
    try:
        name = data.object_names[id(formdata.original_form)]
        return name
    except:
        for index, form in enumerate(data.forms):
            if formdata.preprocessed_form.equals(form):
                name = str(index)
                return name
    # If the form has no name and can not be found in data.forms something went wrong
    assert False


def name_localoperator_file(formdata, data):
    from dune.perftool.options import get_option
    if len(data.forms) == 1:
        filename = get_option("operator_file")
    else:
        suffix = '_' + name_form(formdata, data)
        basename, extension = splitext(get_option("operator_file"))
        filename = basename + suffix + extension
    return filename


@template_parameter(classtag="operator")
def lop_template_ansatz_gfs():
    return "GFSU"


@template_parameter(classtag="operator")
def lop_template_test_gfs():
    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
    gfs = lop_template_ansatz_gfs()
    return "using {} = typename {}::Traits::GridView::ctype;".format(name, gfs)


def name_domain_field():
    name = "DF"
    lop_domain_field(name)
    return name


def lop_template_gfs(ma):
    from ufl.classes import Argument, Coefficient
    if isinstance(ma.argexpr, Argument):
        if ma.argexpr.number() == 0:
            return lop_template_test_gfs()
        if ma.argexpr.number() == 1:
            return lop_template_ansatz_gfs()
    if isinstance(ma.argexpr, Coefficient):
        # Index 0 is reserved for trialfunction, index 1 is reserved for jacobian apply function
        assert ma.argexpr.count() < 2
        return lop_template_ansatz_gfs()
    assert False


def name_initree_constructor_param():
    return "iniParams"


@class_member(classtag="operator")
def define_initree(name):
    param_name = name_initree_constructor_param()
    include_file('dune/common/parametertree.hh', filetag="operatorfile")
    constructor_parameter("const Dune::ParameterTree&", param_name, classtag="operator")
    initializer_list(name, [param_name], classtag="operator")

    return "const Dune::ParameterTree& {};".format(name)


def ufl_measure_to_pdelab_measure(which):
    return {'cell': 'Volume',
            'exterior_facet': 'Boundary',
            'interior_facet': 'Skeleton',
            }.get(which)


@class_member(classtag="operator")
def _enum_pattern(which):
    return "enum {{ doPattern{} = true }};".format(which)


def enum_pattern():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    return _enum_pattern(ufl_measure_to_pdelab_measure(integral_type))


def _pattern_baseclass(measure):
    return base_class('Dune::PDELab::Full{}Pattern'.format(measure), classtag="operator")


def pattern_baseclass():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    return _pattern_baseclass(ufl_measure_to_pdelab_measure(integral_type))


@class_member(classtag="operator")
def _enum_alpha(which):
    return "enum {{ doAlpha{} = true }};".format(which)


def enum_alpha():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    return _enum_alpha(ufl_measure_to_pdelab_measure(integral_type))


def name_initree_member():
    define_initree("_iniParams")
    return "_iniParams"


@class_basename(classtag="operator")
def localoperator_basename(formdata, data):
    form_name = name_form(formdata, data)
    return "LocalOperator" + form_name.capitalize()


def class_type_from_cache(classtag):
    from dune.perftool.generation import retrieve_cache_items

    # get the basename
    basename = [i for i in retrieve_cache_items(condition="{} and basename".format(classtag))]
    assert len(basename) == 1
    basename = basename[0]

    # get the template parameters
    tparams = [i for i in retrieve_cache_items(condition="{} and template_param".format(classtag))]
    tparam_str = ''
    if len(tparams) > 0:
        tparam_str = '<{}>'.format(', '.join(t for t in tparams))

    return basename, basename + tparam_str


def assembler_routine_name():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    form_type = get_global_context_value("form_type")

    part1 = {"residual": "alpha"}.get(form_type, form_type)
    part2 = ufl_measure_to_pdelab_measure(integral_type).lower()

    return "{}_{}".format(part1, part2)


def assembly_routine_signature(formdata):
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    form_type = get_global_context_value("form_type")

    # Check if form is linear
    from dune.perftool.pdelab.driver import is_linear
    linear = is_linear(formdata.original_form)

    if form_type == 'residual':
        if integral_type == 'cell':
            from dune.perftool.pdelab.signatures import alpha_volume_signature
            return alpha_volume_signature()
        if integral_type == 'exterior_facet':
            from dune.perftool.pdelab.signatures import alpha_boundary_signature
            return alpha_boundary_signature()
        if integral_type == 'interior_facet':
            from dune.perftool.pdelab.signatures import alpha_skeleton_signature
            return alpha_skeleton_signature()

    if form_type == 'jacobian':
        if integral_type == 'cell':
            from dune.perftool.pdelab.signatures import jacobian_volume_signature
            return jacobian_volume_signature()
        if integral_type == 'exterior_facet':
            from dune.perftool.pdelab.signatures import jacobian_boundary_signature
            return jacobian_boundary_signature()
        if integral_type == 'interior_facet':
            from dune.perftool.pdelab.signatures import jacobian_skeleton_signature
            return jacobian_skeleton_signature()

    if form_type == 'jacobian_apply':
        if linear:
            if integral_type == 'cell':
                from dune.perftool.pdelab.signatures import jacobian_apply_volume_signature
                return jacobian_apply_volume_signature()
            if integral_type == 'exterior_facet':
                from dune.perftool.pdelab.signatures import jacobian_apply_boundary_signature
                return jacobian_apply_boundary_signature()
            if integral_type == 'interior_facet':
                from dune.perftool.pdelab.signatures import jacobian_apply_skeleton_signature
                return jacobian_apply_skeleton_signature()
        else:
            if integral_type == 'cell':
                from dune.perftool.pdelab.signatures import nonlinear_jacobian_apply_volume_signature
                return nonlinear_jacobian_apply_volume_signature()
            if integral_type == 'exterior_facet':
                from dune.perftool.pdelab.signatures import nonlinear_jacobian_apply_boundary_signature
                return nonlinear_jacobian_apply_boundary_signature()
            if integral_type == 'interior_facet':
                from dune.perftool.pdelab.signatures import nonlinear_jacobian_apply_skeleton_signature
                return nonlinear_jacobian_apply_skeleton_signature()

    assert False


class AccumulationSpace(Record):
    def __init__(self,
                 lfs=None,
                 index=None,
                 restriction=None,
                 element=None,
                 ):
        Record.__init__(self,
                        lfs=lfs,
                        index=index,
                        restriction=restriction,
                        element=element,
                        )

    def get_args(self):
        if self.lfs is None:
            return ()
        else:
            return (self.lfs, self.index)

    def get_restriction(self):
        if self.restriction is None:
            return ()
        else:
            return (self.restriction,)


def determine_accumulation_space(expr, number, measure):
    from dune.perftool.ufl.modified_terminals import extract_modified_arguments
    args = extract_modified_arguments(expr, argnumber=number)

    if measure == 'exterior_facet':
        for ma in args:
            ma.restriction = Restriction.NEGATIVE

    # If this is a residual term we return a dummy object
    if len(args) == 0:
        return AccumulationSpace()

    # There should be but one modified argument, as the splitting eliminated all others.
    assert(len(args) == 1)
    ma, = args

    # Extract information on the finite element
    from ufl.functionview import select_subelement
    subel = select_subelement(ma.argexpr.ufl_element(), ma.component)

    # And generate a local function space for it!
    from dune.perftool.pdelab.spaces import name_lfs, name_lfs_bound, lfs_child, lfs_iname
    lfs = name_lfs(ma.argexpr.ufl_element(), ma.restriction, ma.component)
    from dune.perftool.generation import valuearg
    from loopy.types import NumpyType
    valuearg(lfs, dtype=NumpyType("str"))

    if len(subel.value_shape()) != 0:
        from dune.perftool.pdelab.geometry import dimension_iname
        idims = tuple(dimension_iname(context='arg', count=i) for i in range(len(subel.value_shape())))
        lfs = lfs_child(lfs, idims, shape=subel.value_shape(), symmetry=subel.symmetry())
        subel = subel.sub_elements()[0]

    lfsi = Variable(lfs_iname(subel, ma.restriction, count=number))

    # If the LFS is not yet a pymbolic expression, make it one
    from pymbolic.primitives import Expression
    if not isinstance(lfs, Expression):
        lfs = Variable(lfs)

    return AccumulationSpace(lfs=lfs,
                             index=lfsi,
                             restriction=ma.restriction,
                             )


def boundary_predicates(expr, measure, subdomain_id):
    predicates = frozenset([])

    if subdomain_id not in ['everywhere', 'otherwise']:
        # We need to reconstruct the subdomain_data parameter of the measure
        # I am *totally* confused as to why this information is not at hand anyway,
        # but conversation with Martin pointed me to dolfin.fem.assembly where this
        # is done in preprocessing with the limitation of only one possible type of
        # modified measure per integral type.

        # Get the original form and inspect the present measures
        from dune.perftool.generation import get_global_context_value
        original_form = get_global_context_value("formdata").original_form

        sd = original_form.subdomain_data()
        assert len(sd) == 1
        subdomains, = list(sd.values())
        domain, = list(sd.keys())
        for k in list(subdomains.keys()):
            if subdomains[k] is None:
                    del subdomains[k]

        # Finally extract the original subdomain_data (which needs to be present!)
        assert measure in subdomains
        subdomain_data = subdomains[measure]

        # Determine the name of the parameter function
        name = get_global_context_value("data").object_names[id(subdomain_data)]

        # Trigger the generation of code for this thing in the parameter class
        from ufl.checks import is_cellwise_constant
        cellwise_constant = is_cellwise_constant(expr)
        from dune.perftool.pdelab.parameter import intersection_parameter_function
        intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int32')

        predicates = predicates.union(['{} == {}'.format(name, subdomain_id)])

    return predicates


@iname
def grad_iname(index, dim):
    from dune.perftool.pdelab.index import name_index
    name = name_index(index)
    domain(name, dim)
    return name


@backend(interface="accum_insn")
def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
    # First we do the tree traversal to get a pymbolic expression representing this expression
    pymbolic_expr = visitor(accterm.term)

    # If this is a gradient, we generate an iname
    additional_inames = frozenset({})
    if accterm.argument.index:
        from ufl.domain import find_geometric_dimension
        dim = find_geometric_dimension(accterm.argument.expr)
        for i in accterm.argument.index._indices:
            if i not in visitor.dimension_indices:
                additional_inames = additional_inames.union(frozenset({grad_iname(i, dim)}))

    # It may happen that an entire accumulation term vanishes. We do nothing in that case
    if pymbolic_expr == 0:
        return

    # We also traverse the test function to get its pymbolic equivalent
    test_expr = visitor(accterm.argument.expr)

    # Combine expression and test function
    from pymbolic.primitives import Product
    pymbolic_expr = Product((pymbolic_expr, test_expr))

    # Collect the lfs and lfs indices for the accumulate call
    test_lfs = determine_accumulation_space(accterm.argument.expr, 0, measure)
    # In the jacobian case, also determine the space for the ansatz space
    ansatz_lfs = determine_accumulation_space(accterm.term, 1, measure)

    from dune.perftool.pdelab.argument import name_accumulation_variable
    accumvar = name_accumulation_variable((ansatz_lfs.get_restriction() + test_lfs.get_restriction()))

    predicates = boundary_predicates(accterm.term, measure, subdomain_id)

    rank = 1 if ansatz_lfs.lfs is None else 2

    from dune.perftool.pdelab.argument import PDELabAccumulationFunction
    from pymbolic.primitives import Call
    expr = Call(PDELabAccumulationFunction(accumvar, rank),
                (ansatz_lfs.get_args() + test_lfs.get_args() + (pymbolic_expr,))
                )

    from dune.perftool.generation import instruction
    from dune.perftool.options import option_switch
    quad_inames = visitor.interface.quadrature_inames()

    instruction(assignees=(),
                expression=expr,
                forced_iname_deps=additional_inames.union(frozenset(visitor.inames).union(frozenset(quad_inames))),
                forced_iname_deps_is_final=True,
                predicates=predicates
                )


def generate_kernel(integrals):
    for integral in integrals:
        integrand = integral.integrand()
        measure = integral.integral_type()
        subdomain_id = integral.subdomain_id()
        subdomain_data = integral.subdomain_data()

        # Maybe make the jacobian inverse diagonal!
        if get_option('diagonal_transformation_matrix'):
            from dune.perftool.ufl.transformations.axiparallel import diagonal_jacobian
            integrand = diagonal_jacobian(integrand)

        # Gather dimension indices
        from dune.perftool.ufl.dimensionindex import dimension_index_mapping
        dimension_indices = dimension_index_mapping(integrand)

        # Generate code for the LFS trees present in the form
        from dune.perftool.ufl.modified_terminals import extract_modified_arguments
        test_ma = extract_modified_arguments(integrand, argnumber=0)
        trial_ma = extract_modified_arguments(integrand, coeffcount=0)
        apply_ma = extract_modified_arguments(integrand, coeffcount=1)

        import itertools
        for ma in itertools.chain(test_ma, trial_ma, apply_ma):
            if measure == 'exterior_facet':
                ma.restriction = Restriction.NEGATIVE

            from dune.perftool.pdelab.spaces import traverse_lfs_tree
            traverse_lfs_tree(ma)

        # Now split the given integrand into accumulation expressions
        from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms
        accterms = split_into_accumulation_terms(integrand)

        # Iterate over the terms and generate a kernel
        for term in accterms:
            # Adjust the index map for the visitor
            from copy import deepcopy
            indexmap = deepcopy(dimension_indices)
            for i, j in term.indexmap.items():
                if i in indexmap:
                    indexmap[j] = indexmap[i]

            # Get a transformer instance for this kernel
            if get_option('sumfact'):
                from dune.perftool.sumfact import SumFactInterface
                interface = SumFactInterface()
            else:
                from dune.perftool.pdelab import PDELabInterface
                interface = PDELabInterface()
            from dune.perftool.ufl.visitor import UFL2LoopyVisitor
            visitor = UFL2LoopyVisitor(interface, measure, indexmap)
            get_backend(interface="accum_insn")(visitor, term, measure, subdomain_id)

    knl = extract_kernel_from_cache("kernel_default")

    # All items with the kernel tags can be destroyed once a kernel has been generated
    from dune.perftool.generation import delete_cache_items
    delete_cache_items("kernel_default")

    return knl


def extract_kernel_from_cache(tag):
    # Extract the information, which is needed to create a loopy kernel.
    # First extracting it, might be useful to alter it before kernel generation.
    from dune.perftool.generation import retrieve_cache_functions, retrieve_cache_items
    from dune.perftool.loopy.target import DuneTarget
    domains = [i for i in retrieve_cache_items("{} and domain".format(tag))]

    if not domains:
        domains = ["{[stupid] : 0<=stupid<1}"]

    instructions = [i for i in retrieve_cache_items("{} and instruction".format(tag))]
    temporaries = {i.name: i for i in retrieve_cache_items("{} and temporary".format(tag))}
    arguments = [i for i in retrieve_cache_items("{} and argument".format(tag))]
    silenced = [l for l in retrieve_cache_items("{} and silenced_warning".format(tag))]
    transformations = [t for t in retrieve_cache_items("{} and transformation".format(tag))]

    # Construct an options object
    from loopy import Options
    opt = Options(ignore_boostable_into=True)

    # Create the kernel
    from loopy import make_kernel, preprocess_kernel
    kernel = make_kernel(domains,
                         instructions,
                         arguments,
                         temporary_variables=temporaries,
                         target=DuneTarget(),
                         options=opt,
                         silenced_warnings=silenced,
                         )

    from loopy import make_reduction_inames_unique
    kernel = make_reduction_inames_unique(kernel)

    # Apply the transformations that were gathered during tree traversals
    for trafo in transformations:
        kernel = trafo[0](kernel, *trafo[1])

    from dune.perftool.loopy import heuristic_duplication
    kernel = heuristic_duplication(kernel)

    # Maybe apply vectorization strategies
    if get_option("vectorize"):
        if get_option("sumfact"):
            # Vectorization of the quadrature loop
            insns = [i.id for i in lp.find_instructions(kernel, lp.match.Tagged("quadvec"))]
            if insns:
                from dune.perftool.sumfact.quadrature import quadrature_inames
                inames = quadrature_inames()

                from dune.perftool.loopy.transformations.collect_rotate import collect_vector_data_rotate
                kernel = collect_vector_data_rotate(kernel, insns, inames)
        else:
            raise NotImplementedError("Only vectorizing sumfactoized code right now!")

    # Now add the preambles to the kernel
    preambles = [(i, p) for i, p in enumerate(retrieve_cache_items("{} and preamble".format(tag)))]
    kernel = kernel.copy(preambles=preambles)

    # Do the loopy preprocessing!
    kernel = preprocess_kernel(kernel)

    return kernel


def name_time_dumper_os():
    return "os"


def name_time_dumper_reset():
    return "reset"


def name_time_dumper_t():
    return "t"


def name_time_dumper_counter():
    return "counter"


def name_time_dumper_exec():
    return "exec"


class TimerMethod(ClassMember):
    def __init__(self):
        os = name_time_dumper_os()
        reset = name_time_dumper_reset()
        t = name_time_dumper_t()
        ex = name_time_dumper_exec()

        # TODO: operator counting only works if alpha_volume_kernel exists
        content = ["template <typename Stream>",
                   "void dump_timers(Stream& {}, char* {}, bool {})".format(os, ex, reset),
                   "{",
                   "  double {} = 0.0;".format(t),
                   "#ifdef ENABLE_COUNTER",
                   "  auto counter = HP_TIMER_OPCOUNTERS(alpha_volume_kernel);",
                   "  counter.reset();",
                   "#endif",
                   ""]
        dump_timers = [i for i in retrieve_cache_items(condition='dump_timers')]
        content.extend(map(lambda x: '  ' + x, dump_timers))
        content.extend(["#ifdef ENABLE_COUNTERS",
                        "  counter.reportOperations({});".format(os),
                        "#endif"])
        content.append("}")
        ClassMember.__init__(self, content)


class LoopyKernelMethod(ClassMember):
    def __init__(self, signature, kernel, add_timings=True, initializer_list=[]):
        from loopy import generate_body
        from cgen import LiteralLines, Block
        content = signature

        # Add initializer list if this is a constructor
        if initializer_list:
            content[-1] = content[-1] + " :"
            for init in initializer_list[:-1]:
                content.append(" " * 4 + init + ",")
            content.append(" " * 4 + initializer_list[-1])

        content.append('{')
        if kernel is not None:
            # Add kernel preamble
            for i, p in kernel.preambles:
                content.append('  ' + p)

            # Start timer
            if add_timings and get_option('timer'):
                timer_name = assembler_routine_name() + '_kernel'
                post_include('HP_DECLARE_TIMER({});'.format(timer_name), filetag='operatorfile')
                content.append('  ' + 'HP_TIMER_START({});'.format(timer_name))
                dump_accumulate_timer(timer_name)

            # Add kernel body
            content.extend(l for l in generate_body(kernel).split('\n')[1:-1])

            # Stop timer
            if add_timings and get_option('timer'):
                content.append('  ' + 'HP_TIMER_STOP({});'.format(timer_name))

        content.append('}')
        ClassMember.__init__(self, content)


def cgen_class_from_cache(tag, members=[]):
    from dune.perftool.generation import retrieve_cache_items

    # Generate the name by concatenating basename and template parameters
    basename, fullname = class_type_from_cache(tag)

    base_classes = [bc for bc in retrieve_cache_items('{} and baseclass'.format(tag))]
    constructor_params = [bc for bc in retrieve_cache_items('{} and constructor_param'.format(tag))]
    il = [i for i in retrieve_cache_items('{} and initializer'.format(tag))]
    pm = [m for m in retrieve_cache_items('{} and member'.format(tag))]
    tparams = [i for i in retrieve_cache_items('{} and template_param'.format(tag))]

    # Construct the constructor
    constructor_knl = extract_kernel_from_cache(tag)
    from dune.perftool.loopy.target import DuneTarget
    constructor_knl = constructor_knl.copy(target=DuneTarget(declare_temporaries=False))
    signature = "{}({})".format(basename, ", ".join(next(iter(p.generate(with_semicolon=False))) for p in constructor_params))
    constructor = LoopyKernelMethod([signature], constructor_knl, add_timings=False, initializer_list=il)

    # Take any temporary declarations from the kernel and make them class members
    target = DuneTarget()
    from loopy.codegen import CodeGenerationState
    codegen_state = CodeGenerationState(kernel=constructor_knl,
                                        implemented_data_info=None,
                                        implemented_domain=None,
                                        implemented_predicates=frozenset(),
                                        seen_dtypes=frozenset(),
                                        seen_functions=frozenset(),
                                        seen_atomic_dtypes=frozenset(),
                                        var_subst_map={},
                                        allow_complex=False,
                                        is_generating_device_code=True,
                                        )
    decls = target.get_device_ast_builder().get_temporary_decls(codegen_state, 0)

    from dune.perftool.cgen import Class
    return Class(basename, base_classes=base_classes, members=[constructor] + members + pm + decls, tparam_decls=tparams)


def generate_localoperator_kernels(formdata, data):
    # Extract the relevant attributes of the form data
    form = formdata.preprocessed_form

    # Reset the generation cache
    from dune.perftool.generation import delete_cache_items
    delete_cache_items()

    # Manage includes and base classes that we always need
    include_file('dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh', filetag="operatorfile")
    include_file('dune/pdelab/localoperator/idefault.hh', filetag="operatorfile")
    include_file('dune/pdelab/localoperator/flags.hh', filetag="operatorfile")
    include_file('dune/pdelab/localoperator/pattern.hh', filetag="operatorfile")

    # Trigger this one once early on to assure that template
    # parameters are set in the right order
    localoperator_basename(formdata, data)
    lop_template_ansatz_gfs()
    lop_template_test_gfs()
    lop_template_range_field()
    from dune.perftool.pdelab.parameter import parameterclass_basename
    parameterclass_basename(formdata, data)

    # Make sure there is always the same constructor arguments (even if parameter class is empty)
    from dune.perftool.pdelab.localoperator import name_initree_member
    name_initree_member()
    from dune.perftool.pdelab.parameter import name_paramclass
    name_paramclass()

    # Add right base classes for stationary/instationary operators
    base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")
    from dune.perftool.pdelab.driver import is_stationary
    if not is_stationary():
        rf = lop_template_range_field()
        base_class('Dune::PDELab::InstationaryLocalOperatorDefaultMethods<{}>'
                   .format(rf), classtag="operator")

        # Create set time method in parameter class
        from dune.perftool.pdelab.parameter import define_set_time_method
        define_set_time_method()

    # Have a data structure collect the generated kernels
    operator_kernels = {}

    with global_context(form_type='residual'):
        # Generate the necessary residual methods
        for measure in set(i.integral_type() for i in form.integrals()):
            with global_context(integral_type=measure):
                enum_pattern()
                pattern_baseclass()
                enum_alpha()
                with global_context(kernel=assembler_routine_name()):
                    kernel = generate_kernel(form.integrals_by_type(measure))

                # 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)
                    base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")

                    # Numerical jacobian initializer list
                    ini = name_initree_member()
                    ini_constructor = name_initree_constructor_param()
                    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"):
                        from dune.perftool.pdelab.driver import is_linear
                        if is_linear(formdata.original_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

    # Generate the necessary jacobian methods
    if not get_option("numerical_jacobian"):
        from ufl import derivative
        jacform = derivative(formdata.original_form, formdata.original_form.coefficients()[0])

        from dune.perftool.ufl.preprocess import preprocess_form
        jacform = preprocess_form(jacform).preprocessed_form

        with global_context(form_type="jacobian"):
            for measure in set(i.integral_type() for i in jacform.integrals()):
                with global_context(integral_type=measure):
                    with global_context(kernel=assembler_routine_name()):
                        kernel = generate_kernel(jacform.integrals_by_type(measure))
                operator_kernels[(measure, 'jacobian')] = kernel

            # Generate dummy functions for those kernels, that vanished in the differentiation process
            # We *could* solve this problem by using lambda_* terms but we do not really want that, so
            # we use empty jacobian assembly methods instead
            alpha_measures = set(i.integral_type() for i in form.integrals())
            jacobian_measures = set(i.integral_type() for i in jacform.integrals())
            for it in alpha_measures - jacobian_measures:
                operator_kernels[(it, 'jacobian')] = None

        # Jacobian apply methods for matrix-free computations
        if get_option("matrix_free"):
            # The apply vector has reserved index 1 so we directly use Coefficient class from ufl
            from ufl import Coefficient
            apply_coefficient = Coefficient(form.coefficients()[0].ufl_element(), 1)

            # Create application of jacobian on vector
            from ufl import action
            jac_apply_form = action(jacform, apply_coefficient)

            # Create kernel for jacobian application
            with global_context(form_type="jacobian_apply"):
                for measure in set(i.integral_type() for i in jac_apply_form.integrals()):
                    with global_context(integral_type=measure):
                        with global_context(kernel=assembler_routine_name()):
                            kernel = generate_kernel(jac_apply_form.integrals_by_type(measure))
                    operator_kernels[(measure, 'jacobian_apply')] = kernel

                    # Generate dummy functions for those kernels, that vanished in the differentiation process
                    # We *could* solve this problem by using lambda_* terms but we do not really want that, so
                    # we use empty jacobian assembly methods instead
                    alpha_measures = set(i.integral_type() for i in form.integrals())
                    jacobian_apply_measures = set(i.integral_type() for i in jac_apply_form.integrals())
                    for it in alpha_measures - jacobian_apply_measures:
                        operator_kernels[(it, 'jacobian_apply')] = None

    # Return the set of generated kernels
    return operator_kernels


def generate_localoperator_file(formdata, kernels, filename):
    operator_methods = []

    # Make generables from the given kernels
    for method, kernel in kernels.items():
        it, ft = method
        with global_context(integral_type=it, form_type=ft):
            signature = assembly_routine_signature(formdata)
            operator_methods.append(LoopyKernelMethod(signature, kernel))

    if get_option('timer'):
        include_file('dune/perftool/common/timer.hh', filetag='operatorfile')
        operator_methods.append(TimerMethod())

    # Write the file!
    from dune.perftool.file import generate_file
    param = cgen_class_from_cache("parameterclass")
    # TODO take the name of this thing from the UFL file
    lop = cgen_class_from_cache("operator", members=operator_methods)
    generate_file(filename, "operatorfile", [param, lop])


def generate_localoperator_basefile(formdatas, data):
    filename = get_option("operator_file")
    for formdata in formdatas:
        lop_filename = name_localoperator_file(formdata, data)
        include_file(lop_filename, filetag="operatorbasefile")

    from dune.perftool.file import generate_file
    generate_file(filename, "operatorbasefile", [])