Alpha volume working, dim_tags bug in jacobian

python/dune/perftool/sumfact/__init__.py

@@ -4,6 +4,7 @@ from dune.perftool.sumfact.quadrature import (quadrature_inames,
 
 from dune.perftool.sumfact.basis import (lfs_inames,
                                          pymbolic_basis,
+                                         pymbolic_reference_gradient,
                                          pymbolic_trialfunction,
                                          pymbolic_trialfunction_gradient,
                                          )
@@ -18,6 +19,9 @@ class SumFactInterface(PDELabInterface):
     def pymbolic_basis(self, element, restriction, number):
         return pymbolic_basis(element, restriction, number)
 
+    def pymbolic_reference_gradient(self, element, restriction, number):
+        return pymbolic_reference_gradient(element, restriction, number)
+
     def pymbolic_trialfunction_gradient(self, element, restriction, component):
         return pymbolic_trialfunction_gradient(element, restriction, component)
 

python/dune/perftool/sumfact/basis.py

@@ -198,3 +198,58 @@ def pymbolic_basis(element, restriction, number):
     evaluate_basis(element, name, restriction)
 
     return prim.Variable(name)
+
+
+@backend(interface="evaluate_grad")
+@cached
+def evaluate_reference_gradient(element, name, restriction):
+    # from dune.perftool.pdelab.basis import name_leaf_lfs
+    # lfs = name_leaf_lfs(element, restriction)
+    # from dune.perftool.pdelab.spaces import name_lfs_bound
+    from dune.perftool.pdelab.geometry import name_dimension
+    temporary_variable(
+        name,
+        shape=(name_dimension(),))
+        # shape=())
+    quad_inames = quadrature_inames()
+    inames = lfs_inames(element, restriction)
+    assert(len(quad_inames) == len(inames))
+
+    theta = name_theta()
+    dtheta = name_dtheta()
+
+
+    # TODO WIP!
+    # Get geometric dimension
+    formdata = get_global_context_value('formdata')
+    dim = formdata.geometric_dimension
+
+    for i in range(dim):
+        calls = [prim.Call(ColMajorAccess(theta), (prim.Variable(i), prim.Variable(j)))
+                 for (i, j) in zip(quad_inames, inames)]
+        calls[i] = prim.Call(ColMajorAccess(dtheta), (prim.Variable(quad_inames[i]), prim.Variable(inames[i])))
+        calls = tuple(calls)
+
+        # assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(0)))
+        assignee = prim.Subscript(prim.Variable(name), (i,))
+        # assignee = prim.Variable(name)
+        expression = prim.Product(calls)
+
+        instruction(assignee=assignee,
+                    expression=expression,
+                    forced_iname_deps=frozenset(quad_inames + inames),
+                    forced_iname_deps_is_final=True,
+        )
+
+
+def pymbolic_reference_gradient(element, restriction, number):
+    assert number == 1
+    # TODO ?
+    #assert element.num_sum_elements() == 0
+
+    # TODO: Change name?
+    name = "js_{}".format(FEM_name_mangling(element))
+    name = restricted_name(name, restriction)
+    evaluate_reference_gradient(element, name, restriction)
+
+    return prim.Variable(name)

python/dune/perftool/sumfact/sumfact.py

@@ -79,6 +79,7 @@ def name_test_function_contribution(test):
     return restricted_name("contrib_{}phi_{}".format(grad, str(count)), test.restriction)
 
 
+# palpo TODO move somewhere else
 from pymbolic.mapper import IdentityMapper
 class IndexReplaceMapper(IdentityMapper):
     """Replace indices in pymbolic expression using a dictionary lookup"""

test/sumfact/poisson/poisson_dg_only_volume.mini

@@ -6,6 +6,7 @@ __sumfact_suffix = normal, sumfact | expand sumf
 
 cells = 1 1
 extension = 1. 1.
+printresidual = 1
 printmatrix = 1
 
 [wrapper.vtkcompare]