switch.py

""" boundary and skeleton integrals come in variants in sum factorization - implement the switch! """

import csv

from dune.codegen.generation import (construct_from_mixins,
                                     get_global_context_value,
                                     global_context,
                                     )
from dune.codegen.pdelab.geometry import world_dimension
from dune.codegen.pdelab.localoperator import generate_kernel, generate_kernels_per_integral
from dune.codegen.pdelab.signatures import (assembly_routine_args,
                                            assembly_routine_signature,
                                            kernel_name,
                                            )
from dune.codegen.options import get_form_option, get_option, form_option_context
from dune.codegen.cgen.clazz import ClassMember


def sumfact_generate_kernels_per_integral(integrals):
    dim = world_dimension()
    measure = get_global_context_value("integral_type")

    if measure == "cell":
        yield generate_kernel(integrals)

    if measure == "exterior_facet":
        # Maybe skip sum factorization on boundary integrals
        if not get_form_option("sumfact_on_boundary"):
            mixin = construct_from_mixins(mixins=[get_form_option("geometry_mixins")])()
            geometry = mixin.nonsumfact_fallback() or get_form_option("geometry_mixins")
            with form_option_context(sumfact=False, geometry_mixins=geometry):
                for k in generate_kernels_per_integral(integrals):
                    yield k
                return

        # Generate all necessary kernels
        for facedir in range(dim):
            for facemod in range(2):
                with global_context(facedir_s=facedir, facemod_s=facemod):
                    yield generate_kernel(integrals)

        # Generate switch statement
        yield generate_exterior_facet_switch()

    if measure == "interior_facet":
        # Generate all necessary kernels
        for facedir_s in range(dim):
            for facemod_s in range(2):
                for facedir_n in range(dim):
                    for facemod_n in range(2):
                        if decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
                            with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n):
                                yield generate_kernel(integrals)

        # Generate switch statement
        yield generate_interior_facet_switch()


def get_kernel_name(facedir_s=None, facemod_s=None, facedir_n=None, facemod_n=None):
    with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n):
        return kernel_name()


def decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
    # If we are not using YaspGrid, all variants need to be realized
    if get_form_option("geometry_mixins") == "sumfact_multilinear":
        # Reduce the variability according to grid info file
        if get_option("grid_info") is not None:
            filename = get_option("grid_info")
            with open(filename) as csv_file:
                csv_reader = csv.reader(csv_file, delimiter=" ")
                for row in csv_reader:
                    if (row[0] == 'skeleton' and
                            facedir_s == int(row[1])) and \
                            (facemod_s == int(row[2])) and \
                            (facedir_n == int(row[3])) and \
                            (facemod_n == int(row[4])):
                        return True
                return False
        else:
            return True

    # The PDELab machineries visit-once policy combined with Yasp avoids any visits
    # with facemod_s being True
    # NB: This is not true anymore for parallel computations, as we would like to
    #     skip computations on the overlap and that requires us to visit intersections
    #     on the right/upper part of the domain from within the domain.
    # if facemod_s:
    #     return False

    # A codim1 entity can never be on the upper resp. lower side of the ref element
    # in both inside and outside cell in a YaspGrid
    if facemod_n == facemod_s:
        return False

    # A codim1 entity has the same orientation in both the embedding in the inside
    # and outside cell for a YaspGrid
    if facedir_n != facedir_s:
        return False

    return True


def generate_exterior_facet_switch():
    # Extract the signature
    signature = assembly_routine_signature()
    args = ", ".join(tuple(a for c, a in assembly_routine_args()))
    dim = world_dimension()

    # Construct the switch statement
    block = []
    block.append("{")
    block.append("  size_t variant = ig.indexInInside();")
    block.append("  switch(variant)")
    block.append("  {")

    for facedir_s in range(dim):
        for facemod_s in range(2):
            block.append("    case {}: {}({}); break;".format(2 * facedir_s + facemod_s,
                                                              get_kernel_name(facedir_s=facedir_s,
                                                                              facemod_s=facemod_s,
                                                                              ),
                                                              args))

    block.append('    default: DUNE_THROW(Dune::Exception, "Variation not implemented.");')
    block.append("  }")
    block.append("}")

    return ClassMember(signature + block)


def generate_interior_facet_switch():
    # Extract the signature
    signature = assembly_routine_signature()
    args = ", ".join(tuple(a for c, a in assembly_routine_args()))
    dim = world_dimension()

    # Construct the switch statement
    block = []
    block.append("{")
    block.append("  size_t variant = ig.indexInOutside() + {} * ig.indexInInside();".format(2 * dim))
    block.append("  switch(variant)")
    block.append("  {")

    for facedir_s in range(dim):
        for facemod_s in range(2):
            for facedir_n in range(dim):
                for facemod_n in range(2):
                    if decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
                        block.append("    case {}: {}({}); break;".format((2 * facedir_s + facemod_s) * (2 * dim) + 2 * facedir_n + facemod_n,
                                                                          get_kernel_name(facedir_s=facedir_s,
                                                                                          facemod_s=facemod_s,
                                                                                          facedir_n=facedir_n,
                                                                                          facemod_n=facemod_n,
                                                                                          ),
                                                                          args))

    block.append('    default: DUNE_THROW(Dune::Exception, "Variation not implemented.");')
    block.append("  }")
    block.append("}")

    return ClassMember(signature + block)