Skip to content
Snippets Groups Projects
geometry.py 9.92 KiB
from dune.perftool import Restriction
from dune.perftool.pdelab import restricted_name
from dune.perftool.generation import (cached,
                                      domain,
                                      get_global_context_value,
                                      iname,
                                      preamble,
                                      temporary_variable,
                                      )
from dune.perftool.pdelab.quadrature import (name_quadrature_position,
                                             name_quadrature_position_in_cell,
                                             quadrature_preamble,
                                             )
from ufl.algorithms import MultiFunction
from pymbolic.primitives import Variable


@iname
def _dimension_iname(context, count):
    if context:
        context = '_' + context
    name = 'idim{}{}'.format(context, str(count))
    formdata = get_global_context_value('formdata')
    dim = formdata.geometric_dimension
    domain(name, dim)
    return name


def dimension_iname(context='', count=0):
    return _dimension_iname(context, count)


def name_element_geometry_wrapper():
    return 'eg'


def type_element_geometry_wrapper():
    return 'EG'


def name_intersection_geometry_wrapper():
    return 'ig'


def type_intersection_geometry_wrapper():
    return 'IG'


def name_geometry_wrapper():
    """ Selects the 'native' geometry wrapper of the kernel """
    from dune.perftool.generation import get_global_context_value
    it = get_global_context_value("integral_type")

    if it == 'cell':
        return name_element_geometry_wrapper()
    if it == 'exterior_facet' or it == 'interior_facet':
        return name_intersection_geometry_wrapper()
    assert False


def type_geometry_wrapper():
    """ Selects the 'native' geometry wrapper of the kernel """
    from dune.perftool.generation import get_global_context_value
    it = get_global_context_value("integral_type")

    if it == 'cell':
        return type_element_geometry_wrapper()
    if it == 'exterior_facet' or it == 'interior_facet':
        return type_intersection_geometry_wrapper()
    assert False


@preamble
def define_restricted_cell(name, restriction):
    ig = name_intersection_geometry_wrapper()
    which = "inside" if restriction == Restriction.NEGATIVE else "outside"
    return "const auto& {} = {}.{}();".format(name,
                                              ig,
                                              which,
                                              )


def name_cell(restriction):
    if restriction == Restriction.NONE:
        eg = name_element_geometry_wrapper()
        return "{}.entity()".format(eg)
    else:
        which = "inside" if restriction == Restriction.NEGATIVE else "outside"
        name = "{}_cell".format(which)
        define_restricted_cell(name, restriction)
        return name


@preamble
def define_cell_geometry(name, restriction):
    cell = name_cell(restriction)
    return "auto {} = {}.geometry();".format(name,
                                             cell
                                             )


def name_cell_geometry(restriction):
    name = restricted_name("cell_geo", restriction)
    define_cell_geometry(name, restriction)
    return name


def type_cell_geometry(restriction):
    if restriction == Restriction.NONE:
        eg = type_element_geometry_wrapper()
        return "{}::Geometry".format(eg)
    else:
        ig = type_intersection_geometry_wrapper()
        return "{}::Entity::Geometry".format(ig)


@preamble
def define_intersection_geometry(name):
    ig = name_intersection_geometry_wrapper()
    return "auto {} = {}.geometry();".format(name,
                                             ig,
                                             )


def name_intersection_geometry():
    define_intersection_geometry("is_geo")
    return "is_geo"


def name_intersection():
    ig = name_intersection_geometry_wrapper()
    return "{}.intersection()".format(ig)


def name_geometry():
    """ Selects the 'native' geometry of the kernel """
    from dune.perftool.generation import get_global_context_value
    it = get_global_context_value("integral_type")

    if it == 'cell':
        return name_cell_geometry(Restriction.NONE)
    if it == 'exterior_facet' or it == 'interior_facet':
        return name_intersection_geometry()
    assert False


@preamble
def define_in_cell_geometry(restriction, name):
    ig = name_intersection_geometry_wrapper()
    which = "In" if restriction == Restriction.NEGATIVE else "Out"
    return "auto {} = {}.geometryIn{}side();".format(name,
                                                     ig,
                                                     which
                                                     )


def name_in_cell_geometry(restriction):
    assert restriction is not Restriction.NONE

    name = "geo_in_{}side".format("in" if restriction is Restriction.NEGATIVE else "out")
    define_in_cell_geometry(restriction, name)
    return name


def apply_in_cell_transformation(name, local, restriction):
    geo = name_in_cell_geometry(restriction)
    return quadrature_preamble("{} = {}.global({});".format(name,
                                                            geo,
                                                            local,
                                                            ),
                               assignees=frozenset({name}),
                               read_variables=frozenset({local}),
                               )


def name_in_cell_coordinates(local, basename, restriction):
    name = "{}_in_{}side".format(basename, "in" if restriction is Restriction.NEGATIVE else "out")
    temporary_variable(name, shape=(name_dimension(),), shape_impl=("fv",))
    apply_in_cell_transformation(name, local, restriction)
    return name


def to_cell_coordinates(local, basename, restriction):
    if restriction == Restriction.NONE:
        return local
    else:
        return name_in_cell_coordinates(local, basename, restriction)


@preamble
def define_dimension(name):
    geo = type_geometry_wrapper()
    return 'const int {} = {}::Entity::dimension;'.format(name, geo)


def name_dimension():
    define_dimension("dim")
    return "dim"


@preamble
def define_intersection_dimension(name):
    geo = type_intersection_geometry_wrapper()
    return "const int {} = {}::Geometry::mydimension;".format(name, geo)


def name_intersection_dimension():
    define_intersection_dimension("idim")
    return "idim"


def evaluate_unit_outer_normal(name):
    ig = name_intersection_geometry_wrapper()
    qp = name_quadrature_position()
    return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
                               assignees=frozenset({name}),
                               read_variables=frozenset({qp}),
                               )


@preamble
def declare_normal(name, shape, shape_impl):
    ig = name_intersection_geometry_wrapper()
    return "auto {} = {}.centerUnitOuterNormal();".format(name, ig)


def name_unit_outer_normal():
    name = "outer_normal"
    temporary_variable(name, shape=(name_dimension(),), decl_method=declare_normal)
    evaluate_unit_outer_normal(name)
    return "outer_normal"


def evaluate_unit_inner_normal(name):
    outer = name_unit_outer_normal()
    return quadrature_preamble("auto {} = -1. * {};".format(name, outer),
                               assignees=frozenset({name}),
                               read_variables=frozenset({outer}),
                               )


def name_unit_inner_normal():
    name = "inner_normal"
    temporary_variable(name, shape=(name_dimension(),), decl_method=declare_normal)
    evaluate_unit_inner_normal(name)
    return "inner_normal"


def type_jacobian_inverse_transposed(restriction):
    geo = type_cell_geometry(restriction)
    return "typename {}::JacobianInverseTransposed".format(geo)


@cached
def define_jacobian_inverse_transposed_temporary(restriction):
    @preamble
    def _define_jacobian_inverse_transposed_temporary(name, shape, shape_impl):
        t = type_jacobian_inverse_transposed(restriction)
        return "{} {};".format(t,
                               name,
                               )
    return _define_jacobian_inverse_transposed_temporary


def define_jacobian_inverse_transposed(name, restriction):
    dim = name_dimension()
    temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
    geo = name_cell_geometry(restriction)
    pos = name_quadrature_position_in_cell(restriction)
    return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
                                                                               geo,
                                                                               pos,
                                                                               ),
                               assignees=frozenset({name}),
                               read_variables=frozenset({pos}),
                               )


def name_jacobian_inverse_transposed(restriction):
    name = restricted_name("jit", restriction)
    define_jacobian_inverse_transposed(name, restriction)
    return name


def define_jacobian_determinant(name):
    temporary_variable(name, shape=())
    geo = name_geometry()
    pos = name_quadrature_position()
    code = "{} = {}.integrationElement({});".format(name,
                                                    geo,
                                                    pos,
                                                    )
    return quadrature_preamble(code,
                               assignees=frozenset({name}),
                               read_variables=frozenset({pos}),
                               )


def name_jacobian_determinant():
    name = 'detjac'
    define_jacobian_determinant(name)
    return name


def name_facet_jacobian_determinant():
    name = 'fdetjac'
    define_jacobian_determinant(name)
    return name