compute integration element

python/dune/perftool/blockstructured/geometry.py

 from dune.perftool.generation import (get_backend,
                                       temporary_variable,
-                                      instruction,
-                                      preamble)
+                                      instruction)
 from dune.perftool.tools import get_pymbolic_basename
 from dune.perftool.options import (get_form_option,
                                    option_switch)
-from dune.perftool.pdelab.geometry import (name_jacobian_determinant,
-                                           world_dimension,
+from dune.perftool.pdelab.geometry import (world_dimension,
                                            local_dimension,
                                            apply_to_global_transformation,
                                            name_facet_jacobian_determinant,
-                                           name_element_corners)
+                                           name_element_corners,
+                                           component_iname)
 from dune.perftool.blockstructured.tools import sub_element_inames
 import pymbolic.primitives as prim
 from loopy.match import Writes
@@ -18,14 +17,49 @@ from loopy.match import Writes
 
 # TODO warum muss within_inames_is_final=True gesetzt werden?
 def compute_jacobian_2d(name):
-    pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+    pymbolic_pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
     corners = name_element_corners()
-    instruction(expression=prim.Subscript(prim.Variable(corners), (0,0)),
-                assignee=prim.Subscript(prim.Variable(name), (0,0)),
-                within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()),
-                within_inames_is_final=True,
-                depends_on=frozenset({Writes(corners), Writes(get_pymbolic_basename(pos))})
-                )
+
+    jac_inames = tuple(component_iname(context="jac", count=i) for i in range(world_dimension()))
+
+    def _corner_diff(first, second, additional_terms=tuple()):
+        simplified_names = tuple("cd" + n.split('_')[2] + n.split('_')[3] for n in additional_terms)
+        name = "corner_diff_"+"_".join((str(first), str(second))+simplified_names)
+        temporary_variable(name, shape_impl=('fv',), shape=(world_dimension(),))
+        diminame = component_iname(context='corner')
+
+        if additional_terms:
+            xs_sum = prim.Sum(tuple(prim.Subscript(prim.Variable(x), (prim.Variable(diminame),)) for x in additional_terms))
+        else:
+            xs_sum = 0
+
+        instruction(expression=prim.Sum((prim.Subscript(prim.Variable(corners), (first, prim.Variable(diminame))),
+                                         -1 * prim.Subscript(prim.Variable(corners), (second, prim.Variable(diminame))),
+                                         -1 * xs_sum)),
+                    assignee=prim.Subscript(prim.Variable(name), prim.Variable(diminame)),
+                    within_inames=frozenset(), within_inames_is_final=True,
+                    depends_on=frozenset({Writes(corners)}))
+        return name
+
+    c = _corner_diff(2, 0)
+    b = _corner_diff(1, 0)
+    a = _corner_diff(3, 0, (b, c))
+
+    expr_jac = [None, None]
+    expr_jac[0] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (1,)),
+                                          prim.Subscript(prim.Variable(a), (prim.Variable(jac_inames[0]),)))),
+                            prim.Subscript(prim.Variable(b), (prim.Variable(jac_inames[0]),))))
+    expr_jac[1] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (0,)),
+                                          prim.Subscript(prim.Variable(a), (prim.Variable(jac_inames[1]),)))),
+                            prim.Subscript(prim.Variable(c), (prim.Variable(jac_inames[1]),))))
+
+    for i, expr in enumerate(expr_jac):
+        instruction(expression=expr,
+                    assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_inames[i]),i)),
+                    within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()),
+                    within_inames_is_final=True,
+                    depends_on=frozenset({Writes(b), Writes(get_pymbolic_basename(pymbolic_pos))})
+                    )
 
 
 
@@ -44,14 +78,18 @@ def name_jacobian_matrix():
 
 
 def pymbolic_compute_determinant(matrix):
-    pass
+    dim = world_dimension()
+    matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
+    expr_determinant = prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])),
+                                 -1*prim.Product((matrix_entry[1][0], matrix_entry[0][1]))))
+    return expr_determinant
 
 
 def define_jacobian_determinant(name):
     temporary_variable(name, shape=())
     jacobian = name_jacobian_matrix()
-    determinant = pymbolic_compute_determinant(jacobian)
-    return instruction(expression=prim.Subscript(prim.Variable(jacobian), (0,0)),
+    expr_determinant = pymbolic_compute_determinant(jacobian)
+    return instruction(expression=prim.Call(prim.Variable("abs"), (expr_determinant,)),
                        assignee=prim.Variable(name),
                        within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()),
                        within_inames_is_final=True,
@@ -67,7 +105,13 @@ def name_jacobian_determinant():
 
 # scale determinant according to the order of the blockstructure
 def pymbolic_jacobian_determinant():
-    return prim.Quotient(prim.Variable(name_jacobian_determinant()),
+    if get_form_option("constant_transformation_matrix"):
+        from dune.perftool.pdelab.geometry import name_jacobian_determinant
+        n_det = name_jacobian_determinant()
+    else:
+        from dune.perftool.blockstructured.geometry import name_jacobian_determinant
+        n_det = name_jacobian_determinant()
+    return prim.Quotient(prim.Variable(n_det),
                          prim.Power(get_form_option("number_of_blocks"), local_dimension()))