Implement facet jacobian determinants

This time only for sumfact, as we need direction dependency.

python/dune/perftool/sumfact/__init__.py

@@ -67,6 +67,10 @@ class SumFactInterface(PDELabInterface):
     def pymbolic_quadrature_weight(self):
         return quadrature_weight(self.visitor)
 
+    #
+    # Geometry related generator functions
+    #
+
     def pymbolic_spatial_coordinate(self):
         import dune.perftool.sumfact.geometry
         ret = get_backend(interface="spatial_coordinate", selector=option_switch("diagonal_transformation_matrix"))(self.visitor.do_predicates, self.visitor)
@@ -83,3 +87,7 @@ class SumFactInterface(PDELabInterface):
         ret, indices = pymbolic_unit_inner_normal(self.visitor.indices)
         self.visitor.indices = indices
         return ret
+
+    def pymbolic_facet_jacobian_determinant(self):
+        from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
+        return pymbolic_facet_jacobian_determinant()

python/dune/perftool/sumfact/geometry.py

 """ Sum factorized geometry evaluations """
 
 from dune.perftool.generation import (backend,
+                                      class_member,
                                       domain,
                                       get_backend,
                                       get_counted_variable,
@@ -16,6 +17,7 @@ from dune.perftool.pdelab.geometry import (local_dimension,
                                            world_dimension,
                                            name_geometry,
                                            )
+from dune.perftool.sumfact.switch import get_facedir
 from dune.perftool.sumfact.symbolic import SumfactKernelInputBase
 from dune.perftool.sumfact.vectorization import attach_vectorization_info
 from dune.perftool.options import get_option, option_switch
@@ -194,3 +196,46 @@ def pymbolic_unit_inner_normal(visitor_indices):
     else:
         from dune.perftool.pdelab.geometry import pymbolic_unit_inner_normal as _norm
         return _norm(), visitor_indices
+
+
+def pymbolic_facet_jacobian_determinant():
+    if get_option("constant_transformation_matrix"):
+        return pymbolic_constant_facet_jacobian_determinant()
+    else:
+        from dune.perftool.pdelab.geometry import pymbolic_facet_jacobian_determinant as _norm
+        return _norm()
+
+
+def pymbolic_constant_facet_jacobian_determinant():
+    facedir = get_facedir(Restriction.NEGATIVE)
+    assert isinstance(facedir, int)
+
+    name = "fdetjac"
+    define_constant_facet_jacobian_determinant(name)
+    globalarg(name, dtype=np.float64, shape=(world_dimension(),))
+    return prim.Subscript(prim.Variable(name), (facedir,))
+
+
+@class_member(classtag="operator")
+def define_constant_facet_jacobian_determinant(name):
+    define_constant_facet_jacobian_determinant_eval(name)
+    from dune.perftool.pdelab.localoperator import lop_template_ansatz_gfs
+    gfst = lop_template_ansatz_gfs()
+    return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
+
+
+@kernel_cached(kernel="operator")
+def define_constant_facet_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()
+    code = ["{",
+            "  auto e = *({}.gridView().template begin<0>());".format(gfs),
+            "  int dir=0;",
+            "  for(auto is = {0}.gridView().ibegin(e); is != {0}.gridView().iend(e); ++(++is), ++dir)".format(gfs),
+            "    {}[dir] = is->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, rft, world_dimension() - 1),
+            "}",
+            ]
+    instruction(code="\n".join(code),
+                kernel="operator",
+                )