From 3849f624bb0604b61c69d39becbe60ccf5ea0976 Mon Sep 17 00:00:00 2001
From: Marcel Koch <marcel.koch@uni-muenster.de>
Date: Wed, 12 Jul 2017 15:37:17 +0200
Subject: [PATCH] Adds transformations of JacobianInverse and
 JacobianDeterminant

---
 .../dune/perftool/blockstructured/__init__.py |  4 +-
 python/dune/perftool/blockstructured/basis.py | 46 +++++++++++++++++++
 .../perftool/blockstructured/quadrature.py    |  2 +-
 .../dune/perftool/blockstructured/spaces.py   | 17 +++++++
 python/dune/perftool/pdelab/__init__.py       |  7 ++-
 python/dune/perftool/pdelab/basis.py          | 32 ++++++-------
 python/dune/perftool/pdelab/localoperator.py  |  8 +++-
 python/dune/perftool/pdelab/spaces.py         |  4 +-
 .../ufl/transformations/blockstructured.py    | 39 ++++++++++++++++
 9 files changed, 137 insertions(+), 22 deletions(-)
 create mode 100644 python/dune/perftool/blockstructured/basis.py
 create mode 100644 python/dune/perftool/blockstructured/spaces.py
 create mode 100644 python/dune/perftool/ufl/transformations/blockstructured.py

diff --git a/python/dune/perftool/blockstructured/__init__.py b/python/dune/perftool/blockstructured/__init__.py
index 85836626..50214de6 100644
--- a/python/dune/perftool/blockstructured/__init__.py
+++ b/python/dune/perftool/blockstructured/__init__.py
@@ -1,3 +1,5 @@
 import dune.perftool.blockstructured.accumulation
 import dune.perftool.blockstructured.quadrature
-import dune.perftool.blockstructured.argument
\ No newline at end of file
+import dune.perftool.blockstructured.argument
+import dune.perftool.blockstructured.spaces
+import dune.perftool.blockstructured.basis
diff --git a/python/dune/perftool/blockstructured/basis.py b/python/dune/perftool/blockstructured/basis.py
new file mode 100644
index 00000000..2f545523
--- /dev/null
+++ b/python/dune/perftool/blockstructured/basis.py
@@ -0,0 +1,46 @@
+from dune.perftool.generation import (backend,
+                                      kernel_cached,
+                                      get_backend,
+                                      instruction,
+                                      temporary_variable)
+from dune.perftool.tools import get_pymbolic_basename
+from dune.perftool.pdelab.driver import FEM_name_mangling
+from dune.perftool.pdelab.restriction import restricted_name
+from dune.perftool.pdelab.basis import (declare_cache_temporary,
+                                        name_localbasis_cache)
+from dune.perftool.pdelab.geometry import world_dimension
+from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
+from dune.perftool.blockstructured.spaces import lfs_inames
+import pymbolic.primitives as prim
+
+
+@backend(interface="evaluate_basis", name="blockstructured")
+@kernel_cached
+def evaluate_basis(leaf_element, name, restriction):
+    temporary_variable(name, shape=(4,), decl_method=declare_cache_temporary(leaf_element, restriction, 'Function'))
+    cache = name_localbasis_cache(leaf_element)
+    qp = pymbolic_quadrature_position_in_cell(restriction)
+    instruction(inames=get_backend("quad_inames")(),
+                code='{} = {}.evaluateFunction({}, lfs.finiteElement().localBasis());'.format(name,
+                                                                                             cache,
+                                                                                             str(qp),
+                                                                                             ),
+                assignees=frozenset({name}),
+                read_variables=frozenset({get_pymbolic_basename(qp)}),
+                )
+
+
+@backend(interface="evaluate_grad", name="blockstructured")
+@kernel_cached
+def evaluate_reference_gradient(leaf_element, name, restriction):
+    temporary_variable(name, shape=(4,1,world_dimension()), decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
+    cache = name_localbasis_cache(leaf_element)
+    qp = pymbolic_quadrature_position_in_cell(restriction)
+    instruction(inames=get_backend("quad_inames")(),
+                code='{} = {}.evaluateJacobian({}, lfs.finiteElement().localBasis());'.format(name,
+                                                                                             cache,
+                                                                                             str(qp),
+                                                                                             ),
+                assignees=frozenset({name}),
+                read_variables=frozenset({get_pymbolic_basename(qp)}),
+                )
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/quadrature.py b/python/dune/perftool/blockstructured/quadrature.py
index b960f8ef..a228c029 100644
--- a/python/dune/perftool/blockstructured/quadrature.py
+++ b/python/dune/perftool/blockstructured/quadrature.py
@@ -43,5 +43,5 @@ def pymbolic_quadrature_position():
 @backend(interface="qp_in_cell", name="blockstructured")
 def pymbolic_quadrature_position_in_cell(restriction):
     from dune.perftool.pdelab.geometry import to_cell_coordinates
-    quad_pos = get_backend("quad_pos", selector=lambda: "blockstructured")()
+    quad_pos = pymbolic_quadrature_position()
     return to_cell_coordinates(quad_pos, restriction)
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/spaces.py b/python/dune/perftool/blockstructured/spaces.py
new file mode 100644
index 00000000..6079d75c
--- /dev/null
+++ b/python/dune/perftool/blockstructured/spaces.py
@@ -0,0 +1,17 @@
+from dune.perftool.generation import (backend,
+                                      domain)
+
+
+@backend(interface="lfs_inames", name="blockstructured")
+def lfs_inames(element, restriction, count=None, context=''):
+    assert not ((context == '') and (count is None))
+    if count is not None:
+        if context != '':
+            context = "{}_{}".format(count, context)
+        else:
+            context = str(count)
+
+    name = "micro_{}_index".format(context)
+    domain(name, 4)
+
+    return name,
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 4d7c050d..32cec784 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -2,6 +2,8 @@
 
 # Trigger some imports that are needed to have all backend implementations visible
 # to the selection mechanisms
+from dune.perftool.generation import (get_backend)
+from dune.perftool.options import get_option
 from dune.perftool.pdelab.argument import (pymbolic_apply_function,
                                            pymbolic_apply_function_gradient,
                                            pymbolic_trialfunction,
@@ -56,7 +58,10 @@ class PDELabInterface(object):
     #
 
     def lfs_inames(self, element, restriction, number=None, context=''):
-        return lfs_inames(element, restriction, number, context)
+        if get_option("blockstructured"):
+            return get_backend("lfs_inames", selector=lambda: "blockstructured")(element, restriction, number, context)
+        else:
+            return get_backend("lfs_inames")(element, restriction, number, context)
 
     #
     # Test and trial function related generator functions
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index ccbb97f7..2475b6ed 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -2,18 +2,14 @@
 
 from dune.perftool.generation import (backend,
                                       class_member,
-                                      generator_factory,
                                       get_backend,
                                       include_file,
                                       instruction,
                                       kernel_cached,
-                                      preamble,
                                       temporary_variable,
                                       )
 from dune.perftool.options import get_option
 from dune.perftool.pdelab.spaces import (lfs_child,
-                                         lfs_iname,
-                                         lfs_inames,
                                          name_leaf_lfs,
                                          name_lfs,
                                          name_lfs_bound,
@@ -21,12 +17,7 @@ from dune.perftool.pdelab.spaces import (lfs_child,
                                          )
 from dune.perftool.pdelab.geometry import (component_iname,
                                            world_dimension,
-                                           name_jacobian_inverse_transposed,
-                                           to_cell_coordinates,
                                            )
-from dune.perftool.pdelab.localoperator import (lop_template_ansatz_gfs,
-                                                lop_template_test_gfs,
-                                                )
 from dune.perftool.tools import (get_pymbolic_basename,
                                  get_pymbolic_indices,
                                  )
@@ -88,8 +79,12 @@ def pymbolic_basis(leaf_element, restriction, number, context=''):
     assert leaf_element.num_sub_elements() == 0
     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)
+    if get_option("blockstructured"):
+        get_backend("evaluate_basis", selector=lambda: "blockstructured")(leaf_element, name, restriction)
+        iname, = get_backend("lfs_inames", selector=lambda: "blockstructured")(leaf_element, restriction, number, context=context)
+    else:
+        get_backend("evaluate_basis")(leaf_element, name, restriction)
+        iname, = get_backend("lfs_inames")(leaf_element, restriction, number, context=context)
 
     return Subscript(Variable(name), (Variable(iname),))
 
@@ -116,8 +111,12 @@ def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
     assert leaf_element.num_sub_elements() == 0
     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)
+    if get_option("blockstructured"):
+        get_backend("evaluate_grad", selector=lambda: "blockstructured")(leaf_element, name, restriction)
+        iname, = get_backend("lfs_inames", selector=lambda: "blockstructured")(leaf_element, restriction, number, context=context)
+    else:
+        get_backend("evaluate_grad")(leaf_element, name, restriction)
+        iname, = get_backend("lfs_inames")(leaf_element, restriction, number, context=context)
 
     return Subscript(Variable(name), (Variable(iname), 0))
 
@@ -155,8 +154,10 @@ def evaluate_coefficient(element, name, container, restriction, tree_path):
     if isinstance(sub_element, (VectorElement, TensorElement)):
         lfs = lfs_child(lfs, idims, shape=shape_as_pymbolic(shape), symmetry=element.symmetry())
 
-    from dune.perftool.pdelab.argument import pymbolic_coefficient
-    coeff = pymbolic_coefficient(container, lfs, index)
+    if get_option('blockstructured'):
+        coeff = get_backend("pymbolic_coefficient", selector=lambda: "blockstructured")(container, lfs, index)
+    else:
+        coeff = get_backend("pymbolic_coefficient")(container, lfs, index)
     assignee = Subscript(Variable(name), tuple(Variable(i) for i in idims))
 
     reduction_expr = Product((coeff, basis))
@@ -193,7 +194,6 @@ def evaluate_coefficient_gradient(element, name, container, restriction, tree_pa
     if isinstance(sub_element, (VectorElement, TensorElement)):
         lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]), symmetry=element.symmetry())
 
-    from dune.perftool.pdelab.argument import pymbolic_coefficient
     if get_option('blockstructured'):
         coeff = get_backend("pymbolic_coefficient", selector=lambda: "blockstructured")(container, lfs, index)
     else:
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 40ae1230..e0659c6a 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -229,7 +229,7 @@ def determine_accumulation_space(expr, number, measure, idims=None):
     subel = select_subelement(ma.argexpr.ufl_element(), ma.tree_path)
 
     # And generate a local function space for it!
-    from dune.perftool.pdelab.spaces import name_lfs, name_lfs_bound, lfs_child, lfs_iname
+    from dune.perftool.pdelab.spaces import name_lfs, name_lfs_bound, lfs_child, lfs_inames
     lfs = name_lfs(ma.argexpr.ufl_element(), ma.restriction, ma.tree_path)
     from dune.perftool.generation import valuearg
     from loopy.types import NumpyType
@@ -242,7 +242,11 @@ def determine_accumulation_space(expr, number, measure, idims=None):
         lfs = lfs_child(lfs, idims, shape=subel.value_shape(), symmetry=subel.symmetry())
         subel = subel.sub_elements()[0]
 
-    lfsi = Variable(lfs_iname(subel, ma.restriction, count=number))
+    if get_option("blockstructured"):
+        iname, = get_backend("lfs_inames", selector=lambda: "blockstructured")(subel, ma.restriction, count=number)
+    else:
+        iname, = get_backend("lfs_inames")(subel, ma.restriction, count=number)
+    lfsi = Variable(iname)
 
     # If the LFS is not yet a pymbolic expression, make it one
     from pymbolic.primitives import Expression
diff --git a/python/dune/perftool/pdelab/spaces.py b/python/dune/perftool/pdelab/spaces.py
index 5a7e2852..e698d68c 100644
--- a/python/dune/perftool/pdelab/spaces.py
+++ b/python/dune/perftool/pdelab/spaces.py
@@ -5,6 +5,7 @@ from dune.perftool.generation import (domain,
                                       generator_factory,
                                       include_file,
                                       preamble,
+                                      backend,
                                       )
 from dune.perftool.pdelab.restriction import restricted_name
 
@@ -41,7 +42,7 @@ def define_lfs_bound(lfs, bound):
 
 
 def name_lfs_bound(lfs):
-    # LFS might either be given by an UFL element or by a string describig its name
+    # LFS might either be given by an UFL element or by a string describing its name
     from ufl import FiniteElementBase
     if isinstance(lfs, FiniteElementBase):
         return name_lfs_bound(name_lfs(lfs))
@@ -217,6 +218,7 @@ def lfs_iname(element, restriction, count=None, context=''):
     return _lfs_iname(element, restriction, context)
 
 
+@backend(interface="lfs_inames")
 def lfs_inames(element, restriction, count=None, context=''):
     return (lfs_iname(element, restriction, count, context),)
 
diff --git a/python/dune/perftool/ufl/transformations/blockstructured.py b/python/dune/perftool/ufl/transformations/blockstructured.py
new file mode 100644
index 00000000..d25c95e1
--- /dev/null
+++ b/python/dune/perftool/ufl/transformations/blockstructured.py
@@ -0,0 +1,39 @@
+from dune.perftool.options import get_option
+from dune.perftool.ufl.transformations import ufl_transformation
+from dune.perftool.ufl.transformations.replace import ReplaceExpression
+from ufl.algorithms import MultiFunction
+from ufl import as_ufl
+from ufl.classes import JacobianInverse, JacobianDeterminant, Product, Division, Indexed
+
+
+class ReplaceReferenceTransformation(MultiFunction):
+    def __init__(self, k):
+        MultiFunction.__init__(self)
+        self.k = k
+        self.visited_jit = False
+
+    def expr(self, o):
+        return self.reuse_if_untouched(o, *tuple(self(op) for op in o.ufl_operands))
+
+#TODO abs uses c abs -> only works for ints!!!
+    def abs(self,o):
+        if isinstance(o.ufl_operands[0], JacobianDeterminant):
+            return Division(o, as_ufl(self.k**2))
+        else:
+            return self.reuse_if_untouched(o, *tuple(self(op) for op in o.ufl_operands))
+
+    def jacobian_determinant(self,o):
+        return Division(o, as_ufl(self.k**2))
+
+    def indexed(self, o):
+        expr = o.ufl_operands[0]
+        multiindex = o.ufl_operands[1]
+        if isinstance(expr, JacobianInverse):
+            return Product(as_ufl(self.k), Indexed(expr, multiindex))
+        else:
+            return self.reuse_if_untouched(o, *tuple(self(op) for op in o.ufl_operands))
+
+
+@ufl_transformation(name="blockstructured")
+def blockstructured(expr):
+    return ReplaceReferenceTransformation(get_option("number_of_blocks"))(expr)
\ No newline at end of file
-- 
GitLab