Overhaul of quadrature and range based for loop generator

python/dune/perftool/loopy/target.py

@@ -67,3 +67,45 @@ class DuneTarget(TargetBase):
             body.append(gen_code.ast)
 
         return str(body), gen_code.implemented_domains
+
+    def get_value_arg_decl(self, name, shape, dtype, is_written):
+        assert shape == ()
+
+        return "blubb"
+
+    def get_global_arg_decl(self, name, shape, dtype, is_written):
+        return "blubb"
+
+    def emit_sequential_loop(self, codegen_state, iname, iname_dtype,
+            static_lbound, static_ubound, inner):
+
+        from loopy.codegen import wrap_in
+
+        # Some of our loops are special as they use PDELab specific concepts.
+        # Fortunately those loops are tied to specific inames.
+        from dune.perftool.pdelab.quadrature import quadrature_iname
+        if iname == quadrature_iname():
+            from cgen import CustomLoop
+            from dune.perftool.pdelab.quadrature import quadrature_loop_statement
+            loop_stmt = quadrature_loop_statement()
+            return wrap_in(CustomLoop,
+                           "for({})".format(loop_stmt),
+                           inner)
+
+
+        # From here on it is the default implementation taken from loopys CTarget.
+        ecm = codegen_state.expression_to_code_mapper
+
+        from loopy.symbolic import aff_to_expr
+
+        from pymbolic.mapper.stringifier import PREC_NONE
+        from cgen import For
+
+        return wrap_in(For,
+                "%s %s = %s"
+                % (self.dtype_to_typename(iname_dtype),
+                    iname, ecm(aff_to_expr(static_lbound), PREC_NONE, "i")),
+                "%s <= %s" % (
+                    iname, ecm(aff_to_expr(static_ubound), PREC_NONE, "i")),
+                "++%s" % iname,
+                inner)

python/dune/perftool/pdelab/quadrature.py

@@ -18,25 +18,31 @@ def quadrature_preamble(code, **kw):
 
 
 @symbol
-def name_quadrature_rule():
-    return "rule"
+def name_quadrature_point():
+    return "qp"
 
 
 @symbol
-def name_quadrature_point():
-    rule = name_quadrature_rule()
-    return "{}.position()".format(rule)
+def name_quadrature_position():
+    qp = name_quadrature_point()
+    return "{}.position()".format(qp)
+
+
+@symbol
+def name_quadrature_weight():
+    qp = name_quadrature_point()
+    return "{}.weight()".format(qp)
 
 
 def define_quadrature_factor(fac):
-    rule = name_quadrature_rule()
+    weight = name_quadrature_weight()
     geo = name_geometry()
-    qp = name_quadrature_point()
-    code = "{} = {}->weight()*{}.integrationElement({});".format(fac,
-                                                                 rule,
-                                                                 geo,
-                                                                 qp,
-                                                                 )
+    pos = name_quadrature_position()
+    code = "{} = {}*{}.integrationElement({});".format(fac,
+                                                       weight,
+                                                       geo,
+                                                       pos,
+                                                       )
     return quadrature_preamble(code, assignees=fac)
 
 
@@ -45,3 +51,22 @@ def name_factor():
     temporary_variable("fac", shape=())
     define_quadrature_factor("fac")
     return "fac"
+
+
+@symbol
+def name_order():
+    # TODO how to determine integration order in a generic fashion here?
+    # Or maybe just read it from the ini file? Or define some mapping from the
+    # finite elements?
+    return "2"
+
+
+@cached
+def quadrature_loop_statement():
+    qp = name_quadrature_point()
+    geo = name_geometry()
+    order = name_order()
+    return "const auto& {} : quadratureRule({}, {})".format(qp,
+                                                            geo,
+                                                            order,
+                                                            )