Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
1 result
Blame Permalink
tools.py 4.33 KiB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107 






from dune.codegen.tools import get_pymbolic_basename
from dune.codegen.generation import (iname,
                                     domain,
                                     temporary_variable,
                                     instruction, globalarg, preamble)
from dune.codegen.pdelab.geometry import world_dimension
from dune.codegen.options import get_form_option
import pymbolic.primitives as prim


# add inames over the micro elements in tensor representation,
# i.e. each element has (i_1,i_2,...,i_d) indices
@iname
def sub_element_inames():
    name = "subel"
    dim = world_dimension()
    dim_names = ["x", "y", "z"] + [str(i) for i in range(4, dim + 1)]
    inames = tuple()
    for i in range(dim):
        inames = inames + ("subel_" + dim_names[i],)
        domain("subel_" + dim_names[i], get_form_option("number_of_blocks"))
    return inames


# compute sequential index for given tensor index, the same as index in base-k to base-10
def tensor_index_to_sequential_index(indices, k):
    return prim.Sum(tuple(prim.Variable(index) * k ** i for i, index in enumerate(indices)))


# compute tensor index for given sequential index, the same as index in base-10 to base-k
def sequential_index_to_tensor_index(iname, k):
    return tuple(prim.Remainder(prim.Variable(iname) / (k**i), k) for i in range(world_dimension()))


# compute index for higher order FEM for a given Q1 index
def micro_index_to_macro_index(element, inames):
    subelem_inames = sub_element_inames()

    k = get_form_option("number_of_blocks")
    p = element.degree()
    return prim.Sum(tuple((p * prim.Variable(si) + prim.Variable(bi)) * (p * k + 1) ** i
                          for i, (si, bi) in enumerate(zip(subelem_inames, inames))))


# translate a point in the micro element into macro coordinates
def define_point_in_macro(name, point_in_micro, visitor):
    # TODO this won't work for 2d mannifolds
    dim = world_dimension()
    if get_form_option('vectorization_blockstructured'):
        temporary_variable(name, shape=(dim,), managed=True)
    else:
        temporary_variable(name, shape=(dim,), shape_impl=('fv',))

    # point_macro = (point_micro + index_micro) / number_of_blocks
    # where index_micro = tensor index of the micro element
    subelem_inames = sub_element_inames()
    for i in range(dim):
        if isinstance(point_in_micro, prim.Subscript):
            expr = prim.Subscript(point_in_micro.aggregate, point_in_micro.index + (i,))
        else:
            expr = prim.Subscript(point_in_micro, (i,))
        expr = prim.Sum((expr, prim.Variable(subelem_inames[i]),))
        expr = prim.Quotient(expr, get_form_option('number_of_blocks'))
        # TODO relax within inames
        instruction(assignee=prim.Subscript(prim.Variable(name), (i,)),
                    expression=expr,
                    within_inames=frozenset(subelem_inames + visitor.quadrature_inames()),
                    tags=frozenset({subelem_inames[i]})
                    )


def name_point_in_macro(point_in_micro, visitor):
    assert isinstance(point_in_micro, prim.Expression)
    name = get_pymbolic_basename(point_in_micro) + "_macro"
    define_point_in_macro(name, point_in_micro, visitor)
    return name


@preamble
def define_container_alias(name, container, lfs, element, is_const):
    k = get_form_option("number_of_blocks")
    p = element.degree()
    dim = world_dimension()
    element_stride = tuple(p * (p * k + 1)**i for i in range(0, dim))
    index_stride = tuple((p * k + 1)**i for i in range(0, dim))
    globalarg(name, shape=(k,) * dim + (p + 1,) * dim, strides=element_stride + index_stride, managed=True)
    if is_const:
        return "const auto {} = &{}({},0);".format(name, container, lfs.name)
    else:
        return "auto {} = &{}.container()({},0);".format(name, container, lfs.name)


def name_accumulation_alias(container, accumspace):
    name = container + "_" + accumspace.lfs.name + "_alias"
    name_tail = container + "_" + accumspace.lfs.name + "_alias_tail"

    define_container_alias(name, container, accumspace.lfs, accumspace.element, is_const=False)
    define_container_alias(name_tail, container, accumspace.lfs, accumspace.element, is_const=False)
    return name


def name_container_alias(container, lfs, element):
    name = container + "_" + lfs.name + "_alias"

    define_container_alias(name, container, lfs, element, is_const=True)
    return name

Datenschutz und Impressum