From 28e6b457d93d8dd9eb111eb80e249059ca3860c9 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Fri, 6 Oct 2017 13:47:18 +0200
Subject: [PATCH] Also implement jacobian determinants on an operator level

---
 python/dune/perftool/pdelab/geometry.py | 52 +++++++++++++++++++------
 1 file changed, 41 insertions(+), 11 deletions(-)

diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 74a60a46..bb590f98 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -350,7 +350,7 @@ def define_jacobian_inverse_transposed(name, restriction):
                                )
 
 
-@backend(interface="define_jit", name="default")
+@backend(interface="name_jit", name="default")
 def name_jacobian_inverse_transposed(restriction):
     name = restricted_name("jit", restriction)
     define_jacobian_inverse_transposed(name, restriction)
@@ -364,18 +364,27 @@ def pymbolic_jacobian_inverse_transposed(i, j, restriction):
     return prim.Subscript(prim.Variable(name), (j, i))
 
 
+@preamble(kernel="operator")
+def define_constant_jacobian_determinant_eval(name):
+    from dune.perftool.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
+    gfs = name_ansatz_gfs_constructor_param()
+    rft = lop_template_range_field()
+
+    return "{} = {}.gridView().template begin<0>()->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
+
+
+@class_member(classtag="operator")
+def _define_constant_jacobian_determinant(name):
+    define_constant_jacobian_determinant_eval(name)
+    from dune.perftool.pdelab.localoperator import lop_template_ansatz_gfs
+    gfst = lop_template_ansatz_gfs()
+    return "typename {}::Traits::GridView::template Codim<0>::Geometry::ctype {};".format(gfst, name)
+
+
 @backend(interface="detjac", name="constant_transformation_matrix")
-@preamble
 def define_constant_jacobian_determinant(name):
-    geo = name_geometry()
-    pos = name_localcenter()
-
     valuearg(name, dtype=np.float64)
-
-    return "auto {} = {}.integrationElement({});".format(name,
-                                                         geo,
-                                                         pos,
-                                                         )
+    _define_constant_jacobian_determinant(name)
 
 
 @backend(interface="detjac", name="default")
@@ -394,6 +403,27 @@ def define_jacobian_determinant(name):
                                )
 
 
+@backend(interface="fdetjac", name="constant_transformation_matrix")
+@preamble
+def define_facet_jacobian_determinant(name):
+    # NB: This might be optimized to store *d* values on the operator level
+    #     We don't do that right now for laziness.
+    geo = name_geometry()
+    pos = name_localcenter()
+
+    valuearg(name, dtype=np.float64)
+
+    return "auto {} = {}.integrationElement({});".format(name,
+                                                         geo,
+                                                         pos,
+                                                         )
+
+
+@backend(interface="fdetjac", name="default")
+def define_facet_jacobian_determinant(name):
+    return define_jacobian_determinant(name)
+
+
 def name_jacobian_determinant():
     name = 'detjac'
     get_backend(interface="detjac", selector=option_switch("constant_transformation_matrix"))(name)
@@ -406,7 +436,7 @@ def pymbolic_jacobian_determinant():
 
 def name_facet_jacobian_determinant():
     name = 'fdetjac'
-    get_backend(interface="detjac", selector=option_switch("constant_transformation_matrix"))(name)
+    get_backend(interface="fdetjac", selector=option_switch("constant_transformation_matrix"))(name)
     return name
 
 
-- 
GitLab