localoperator.py

        # 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):
                    from dune.codegen.pdelab.signatures import assembler_routine_name
                    with global_context(kernel=assembler_routine_name()):
                        kernel = [k for k in generate_kernels_per_integral(jac_apply_form.integrals_by_type(measure))]

                        if kernel:
                            enum_pattern()
                            pattern_baseclass()
                            enum_alpha()

                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:
                    with global_context(integral_type=it):
                        from dune.codegen.pdelab.signatures import assembly_routine_signature
                        operator_kernels[(it, 'jacobian_apply')] = [LoopyKernelMethod(assembly_routine_signature(), kernel=None)]
    if get_form_option("generate_jacobians"):
        with global_context(form_type="jacobian"):
            with form_option_context(conditional=get_form_option("sumfact") and get_form_option("sumfact_regular_jacobians"),
                                     geometry_mixins="generic",
                                     sumfact=False):
                for measure in set(i.integral_type() for i in jacform.integrals()):
                    if not measure_is_enabled(measure):
                        continue

                    logger.info("generate_jacobian_kernels: measure {}".format(measure))
                    with global_context(integral_type=measure):
                        from dune.codegen.pdelab.signatures import assembler_routine_name
                        with global_context(kernel=assembler_routine_name()):
                            kernel = [k for k in generate_kernels_per_integral(jacform.integrals_by_type(measure))]

                            if kernel:
                                enum_pattern()
                                pattern_baseclass()
                                enum_alpha()

                    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:
                    with global_context(integral_type=it):
                        from dune.codegen.pdelab.signatures import assembly_routine_signature
                        operator_kernels[(it, 'jacobian')] = [LoopyKernelMethod(assembly_routine_signature(), kernel=None)]

    return operator_kernels


def generate_control_kernels(forms):
    # All forms will we written in the residual method and
    # accumulation will be done in a class member instead of the
    # residual.
    logger = logging.getLogger(__name__)
    with global_context(form_type='residual'):
        operator_kernels = {}

        # Generate the necessary residual methods
        for measure in set(i.integral_type() for form in forms for i in form.integrals()):
            if not measure_is_enabled(measure):
                continue

            logger.info("generate_control_kernels: measure {}".format(measure))
            with global_context(integral_type=measure):
                from dune.codegen.pdelab.signatures import assembler_routine_name
                with global_context(kernel=assembler_routine_name()):
                    # TODO: Sumfactorization not yet implemented
                    assert not get_form_option('sumfact')

                    from dune.codegen.pdelab.adjoint import control_generate_kernels_per_integral
                    forms_measure = [form.integrals_by_type(measure) for form in forms]
                    kernel = [k for k in control_generate_kernels_per_integral(forms_measure)]

                    # The integrals might vanish due to unfulfillable boundary conditions.
                    # We only generate the local operator enums/base classes if they did not.
                    if kernel:
                        enum_pattern()
                        pattern_baseclass()
                        enum_alpha()

            operator_kernels[(measure, 'residual')] = kernel

        return operator_kernels


def generate_localoperator_kernels(operator):
    logger = logging.getLogger(__name__)

    data = get_global_context_value("data")
    original_form = data.object_by_name[get_form_option("form")]
    from dune.codegen.ufl.preprocess import preprocess_form

    if get_form_option("adjoint"):
        # Generate adjoint operator
        #
        # The jacobian of the adjoint form is just the jacobian of the
        # original form with test and ansazt function swapped. A a
        # linear form you have to subtract the derivative of the
        # objective function w.r.t the ansatz function to get the
        # final residual formulation of the adjoint.
        #
        # Might not be true in all cases but works for the simple ones.
        assert get_form_option("objective_function") is not None
        assert get_form_option("control") is False

        from ufl import derivative, adjoint, action, replace
        from ufl.classes import Coefficient

        # Jacobian of the adjoint form
        jacform = derivative(original_form, original_form.coefficients()[0])
        adjoint_jacform = adjoint(jacform)

        # Derivative of objective function w.r.t. state
        objective = data.object_by_name[get_form_option("objective_function")]
        objective_jacobian = derivative(objective, objective.coefficients()[0])

        # Replace coefficient belonging to ansatz function with new coefficient
        element = objective.coefficients()[0].ufl_element()
        coeff = Coefficient(element, count=3)
        objective_jacobian = replace(objective_jacobian, {objective.coefficients()[0]: coeff})
        if len(adjoint_jacform.coefficients()) > 0:
            adjoint_jacform = replace(adjoint_jacform, {adjoint_jacform.coefficients()[0]: coeff})

        # Residual of the adjoint form
        adjoint_form = action(adjoint_jacform, original_form.coefficients()[0])
        adjoint_form = adjoint_form + objective_jacobian

        # Update form and original_form
        original_form = adjoint_form
        form = preprocess_form(adjoint_form).preprocessed_form

    elif get_form_option("control"):
        # Generate control operator
        #
        # This is the normal form derived w.r.t. the control
        # variable. We generate a form for every row of:
        #
        # \nabla  \hat{J}(m) = (\nabla R(z,m))^T \lambda + \nabla_m J(z,m)
        #
        # These forms will not depend on the test function anymore and
        # will need special treatment for the accumulation process.
        from ufl import action, diff
        from ufl.classes import Coefficient

        # Get control variables
        assert get_form_option("control_variable") is not None
        controls = [data.object_by_name[ctrl.strip()] for ctrl in get_form_option("control_variable").split(",")]

        # Transoform flat index to multiindex. Wrapper around numpy
        # unravel since we need to transform numpy ints to native
        # ints.
        def _unravel(flat_index, shape):
            multi_index = np.unravel_index(flat_index, shape)
            multi_index = tuple(int(i) for i in multi_index)
            return multi_index

        # Will be used to replace ansatz function with adjoint function
        element = original_form.coefficients()[0].ufl_element()
        coeff = Coefficient(element, count=3)

        # Store a form for every control
        forms = []
        for control in controls:
            shape = control.ufl_shape
            flat_length = int(np.prod(shape))
            for i in range(flat_length):
                c = control[_unravel(i, shape)]
                control_form = diff(original_form, c)
                control_form = action(control_form, coeff)
                objective = data.object_by_name[get_form_option("objective_function")]
                objective_gradient = diff(objective, c)
                control_form = control_form + objective_gradient
                forms.append(preprocess_form(control_form).preprocessed_form)

        # Used to create local operator default settings
        form = preprocess_form(original_form).preprocessed_form

    else:
        form = preprocess_form(original_form).preprocessed_form

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

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

    # Generate things needed for all local operator files
    local_operator_default_settings(operator, form)

    if get_form_option("control"):
        logger.info("generate_localoperator_kernels: create methods for control operator")
        operator_kernels.update(generate_control_kernels(forms))
    else:
        logger.info("generate_localoperator_kernels: create residual methods")
        operator_kernels.update(generate_residual_kernels(form, original_form))

        # Generate the necessary jacobian methods
        if not get_form_option("numerical_jacobian"):
            logger.info("generate_localoperator_kernels: create jacobian methods")
            operator_kernels.update(generate_jacobian_kernels(form, original_form))

    # Return the set of generated kernels
    return operator_kernels


def generate_localoperator_file(kernels, filename):
    logger = logging.getLogger(__name__)

    operator_methods = []
    for k in kernels.values():
        operator_methods.extend(k)

    # Generate all the realizations of sum factorization kernel objects needed in this operator
    sfkernels = [sf for sf in retrieve_cache_items("kernelimpl")]
    if sfkernels:
        logger.info("generate_localoperator_kernels: Create {} sumfact kernel realizations".format(len(sfkernels)))

    from dune.codegen.sumfact.realization import realize_sumfact_kernel_function
    for sf, qp in sfkernels:
        from dune.codegen.sumfact.tabulation import set_quadrature_points
        set_quadrature_points(qp)
        operator_methods.append(realize_sumfact_kernel_function(sf))

    if get_option('instrumentation_level') >= 3:
        include_file('dune/codegen/common/timer.hh', filetag='operatorfile')
        if get_option('use_likwid'):
            operator_methods.append(RegisterLikwidMethod())
        elif get_option('use_sde'):
            operator_methods.append(RegisterSSCMarksMethod())
        else:
            operator_methods.append(TimerMethod())
    elif get_option('opcounter'):
        include_file('dune/codegen/common/timer.hh', filetag='operatorfile')

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