From 846f46fb24fa02630dbaabce9171109a0b8b029c Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Wed, 20 Apr 2016 15:38:57 +0200
Subject: [PATCH] Laplacian runs!

---
 python/dune/perftool/generation/loopy.py |  5 +--
 python/dune/perftool/loopy/target.py     | 10 +++++-
 python/dune/perftool/pdelab/basis.py     | 39 +++++++++++++++---------
 python/dune/perftool/pdelab/geometry.py  | 14 ++++++---
 4 files changed, 47 insertions(+), 21 deletions(-)

diff --git a/python/dune/perftool/generation/loopy.py b/python/dune/perftool/generation/loopy.py
index ee43debc..8214ec76 100644
--- a/python/dune/perftool/generation/loopy.py
+++ b/python/dune/perftool/generation/loopy.py
@@ -64,11 +64,12 @@ def temporary_variable(name, **kwargs):
     if 'dtype' not in kwargs:
         kwargs['dtype'] = numpy.float64
 
-    decl_method = kwargs.pop('decl_method', default_declaration)
     shape = kwargs.get('shape', ())
     shape_impl = kwargs.pop('shape_impl', ('arr',) * len(shape))
+    decl_method = kwargs.pop('decl_method', default_declaration)
 
-    decl_method(name, shape, shape_impl)
+    if decl_method is not None:
+        decl_method(name, shape, shape_impl)
 
     return loopy.TemporaryVariable(name, **kwargs)
 
diff --git a/python/dune/perftool/loopy/target.py b/python/dune/perftool/loopy/target.py
index 7395bbe2..2553dcf9 100644
--- a/python/dune/perftool/loopy/target.py
+++ b/python/dune/perftool/loopy/target.py
@@ -18,7 +18,15 @@ _registry = {'float32': 'float',
 
 
 class MyMapper(LoopyCCodeMapper):
-    var_subst_map = {}
+    def map_subscript(self, expr, enclosing_prec, type_context):
+        ret = str(expr.aggregate)
+        from pymbolic.primitives import Variable
+        if isinstance(expr.index, Variable):
+            ret = ret + '[{}]'.format(str(expr.index))
+        else:
+            for i in expr.index:
+                ret = ret + '[{}]'.format(str(i))
+        return ret
 
 
 class DuneTarget(TargetBase):
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 00042f14..06b0a797 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -11,7 +11,7 @@ from dune.perftool.generation import (cached,
                                       symbol,
                                       temporary_variable,
                                       )
-from dune.perftool.pdelab.quadrature import (name_quadrature_point,
+from dune.perftool.pdelab.quadrature import (name_quadrature_position,
                                              quadrature_iname,
                                              )
 from dune.perftool.pdelab.geometry import (name_dimension,
@@ -201,10 +201,10 @@ def name_basis(element):
 
 @cached
 def evaluate_reference_gradient(element, name):
-    temporary_variable(name, shape=(name_lfs_bound(element), name_dimension()), shape_impl=('vec', 'fv'))
+    temporary_variable(name, shape=(name_lfs_bound(element), name_dimension()), decl_method=None)
     cache = name_localbasis_cache(element)
     lfs = name_lfs(element)
-    qp = name_quadrature_point()
+    qp = name_quadrature_position()
     instruction(inames=(quadrature_iname(),
                         ),
                 code='auto& {} = {}.evaluateJacobian({}, {}.finiteElement().localBasis());'.format(name,
@@ -232,11 +232,12 @@ def evaluate_basis_gradient(element, name):
     instruction(inames=(index,
                         quadrature_iname(),
                         ),
-                code='{}.mv({}[{}][0], {});'.format(jac,
-                                                    reference_gradients,
-                                                    index,
-                                                    name,
-                                                    ),
+                code='{}.mv({}[{}][0], {}[{}]);'.format(jac,
+                                                        reference_gradients,
+                                                        index,
+                                                        name,
+                                                        index,
+                                                        ),
                 assignees=name,
                 read_variables=frozenset({reference_gradients}),
                 )
@@ -265,27 +266,37 @@ def evaluate_trialfunction(element, name):
                                                       ),
                 assignees=frozenset({name}),
                 read_variables=frozenset({basis}),
-
                 )
 
 
+def reset_trialfunction_gradient(name):
+    return instruction(inames=(quadrature_iname(),
+                               ),
+                       code='{} = 0.0;'.format(name),
+                       assignees=frozenset({name}),
+                       )
+
+
 @cached
 def evaluate_trialfunction_gradient(element, name):
     # TODO this is of course not yet correct
     temporary_variable(name, shape=(name_dimension(),), shape_impl=('fv',))
+    reset = reset_trialfunction_gradient(name)
     lfs = name_lfs(element)
     index = lfs_iname(element, context='trialgrad')
     basis = name_basis_gradient(element)
     instruction(inames=(quadrature_iname(),
                         index,
                         ),
-                code='{}.axpy(x({}, {}), {});'.format(name,
-                                                      lfs,
-                                                      index,
-                                                      basis
-                                                      ),
+                code='{}.axpy(x({}, {}), {}[{}]);'.format(name,
+                                                          lfs,
+                                                          index,
+                                                          basis,
+                                                          index,
+                                                          ),
                 assignees=frozenset({name}),
                 read_variables=frozenset({basis}),
                 forced_iname_deps=frozenset({quadrature_iname(), index}),
                 forced_iname_deps_is_final=True,
+                depends_on=frozenset({reset}),
                 )
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index a5e24dd0..e22c87bd 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -35,12 +35,18 @@ def name_dimension():
     return "dim"
 
 
+@symbol
+def type_jacobian_inverse_transposed():
+    geo = type_geometry()
+    return "typename {}::Geometry::JacobianInverseTransposed".format(geo)
+
+
 @preamble
 def define_jacobian_inverse_transposed_temporary(name, shape, shape_impl):
-    geo = name_geometry()
-    return "auto {} = {}.jacobianInverseTransposed({{{{ 0.0 }}}});".format(name,
-                                                                           geo,
-                                                                           )
+    t = type_jacobian_inverse_transposed()
+    return "{} {};".format(t,
+                           name,
+                           )
 
 
 def define_jacobian_inverse_transposed(name):
-- 
GitLab