uses nested indices for basis functions, i.e. phi[ix+2*iy]

python/dune/perftool/blockstructured/basis.py

@@ -21,6 +21,7 @@ from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
 from dune.perftool.pdelab.spaces import type_leaf_gfs
 from dune.perftool.pdelab.restriction import restricted_name
 from dune.perftool.blockstructured.spaces import lfs_inames
+from dune.perftool.blockstructured.tools import tensor_index_to_sequential_index
 
 from ufl import MixedElement
 
@@ -86,9 +87,9 @@ def pymbolic_basis(leaf_element, restriction, number, context=''):
     name = "phi_{}".format(FEM_name_mangling(leaf_element))
     name = restricted_name(name, restriction)
     evaluate_basis(leaf_element, name, restriction)
-    iname = lfs_inames(leaf_element, restriction, number, context=context)[0]
+    inames = lfs_inames(leaf_element, restriction, number, context=context)
 
-    return prim.Subscript(prim.Variable(name), (prim.Variable(iname), 0))
+    return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames,leaf_element.degree()+1),0))
 
 
 @kernel_cached
@@ -109,6 +110,6 @@ def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
     name = "js_{}".format(FEM_name_mangling(leaf_element))
     name = restricted_name(name, restriction)
     evaluate_reference_gradient(leaf_element, name, restriction)
-    iname = lfs_inames(leaf_element, restriction, number, context=context)[0]
+    inames = lfs_inames(leaf_element, restriction, number, context=context)
 
-    return prim.Subscript(prim.Variable(name), (prim.Variable(iname), 0))
+    return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames,leaf_element.degree()+1),0))

python/dune/perftool/blockstructured/spaces.py

@@ -14,7 +14,10 @@ def lfs_inames(element, restriction, count=None, context=''):
 
     lfs = name_leaf_lfs(element, restriction)
 
-    name = "micro_{}_{}_index".format(lfs, context)
-    domain(name, pow(element.degree() + 1, world_dimension()))
-
-    return (name, )
+    dim_names = ["x","y","z"] + [str(i) for i in range(4,world_dimension()+1)]
+    name = "micro_{}_{}_index_".format(lfs, context)
+    inames = tuple()
+    for i in range(world_dimension()):
+        inames = inames + (name+dim_names[i],)
+        domain(name+dim_names[i], element.degree() + 1)
+    return inames

python/dune/perftool/blockstructured/tools.py

@@ -96,16 +96,16 @@ def sub_facet_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(index * k ** i for i, index in enumerate(indices)))
+    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 to_tensor_index(iname, 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, iname):
+def micro_index_to_macro_index(element, inames):
     it = get_global_context_value("integral_type")
     if it == "cell":
         subelem_inames = sub_element_inames()
@@ -114,7 +114,5 @@ def micro_index_to_macro_index(element, iname):
 
     k = get_option("number_of_blocks")
     p = element.degree()
-    modified_index = prim.Sum((tensor_index_to_sequential_index(to_tensor_index(iname, p + 1), p * k + 1),
-                               tensor_index_to_sequential_index(tuple(prim.Variable(iname) * p for iname in subelem_inames), p * k + 1)))
-
-    return modified_index
+    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))))