diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index dfcaeba86bdd21bff5b956aa330f016243e2d7c6..7c0b169a108c2b1063d2c0bbc2fcc2ebea90c899 100644
--- a/python/dune/perftool/compile.py
+++ b/python/dune/perftool/compile.py
@@ -36,29 +36,13 @@ def read_ufl(uflfile):
formdatas = []
forms = data.forms
for index, form in enumerate(forms):
- from dune.perftool.ufl.execution import (JacobianDeterminant,
- JacobianInverse)
- formdatas.append(compute_form_data(form,
- do_apply_function_pullbacks=True,
- do_apply_geometry_lowering=True,
- do_apply_integral_scaling=True,
- preserve_geometry_types=(JacobianDeterminant, JacobianInverse)
- )
- )
+ from dune.perftool.ufl.preprocess import preprocess_form
+ formdatas.append(preprocess_form(form))
forms[index] = formdatas[index].preprocessed_form
# We expect at least one form
assert len(data.forms) >= 1
- # apply some transformations unconditionally!
- from dune.perftool.ufl.transformations import transform_form
- from dune.perftool.ufl.transformations.indexpushdown import pushdown_indexed
- from dune.perftool.ufl.transformations.reindexing import reindexing
- for i in range(len(forms)):
- forms[i] = transform_form(forms[i], pushdown_indexed)
- forms[i] = transform_form(forms[i], reindexing)
- formdatas[i].preprocessed_form = forms[i]
-
return formdatas, data
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index bb9efa516b971c03007322216aa59e5cff34523b..ae3396eb0f011d6688c09173683abe2424534941 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -204,7 +204,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
if self.reference_grad:
from dune.perftool.pdelab.basis import name_reference_gradient
- return Subscript(Variable(name_reference_gradient(leaf_element, restriction)), (Variable(iname),))
+ return Subscript(Variable(name_reference_gradient(leaf_element, restriction)), (Variable(iname), 0))
else:
from dune.perftool.pdelab.basis import name_basis
return Subscript(Variable(name_basis(leaf_element, restriction)), (Variable(iname),))
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 416e1cbb29b000a3c1815445faee7b9e854bec49..88b5eed964a2d9c5463dbdead355987ab9f7bdc6 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -365,7 +365,7 @@ def evaluate_coefficient_gradient(element, name, container, restriction, compone
coeff = pymbolic_coefficient(container, lfs, index)
assignee = Subscript(Variable(name), tuple(Variable(i) for i in idims))
- reduction_expr = Product((coeff, Subscript(Variable(basis), (Variable(index), Variable(idims[-1])))))
+ reduction_expr = Product((coeff, Subscript(Variable(basis), (Variable(index), 0, Variable(idims[-1])))))
instruction(expression=Reduction("sum", index, reduction_expr, allow_simultaneous=True),
assignee=assignee,
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 143a1356214859963345328a0529c07e9800156d..1834f4a4fd2be7f08f04e58bcf80d7d5c4564ee1 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -50,12 +50,19 @@ class GeometryMapper(MultiFunction):
return Variable(name_jacobian_determinant())
def jacobian_inverse(self, o):
+ restriction = self.restriction
+ if self.measure == 'exterior_facet':
+ restriction = Restriction.NEGATIVE
+
self.transpose_necessary = True
- return Variable(name_jacobian_inverse_transposed(self.restriction))
+ return Variable(name_jacobian_inverse_transposed(restriction))
def jacobian(self, o):
raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
+ def facet_jacobian_determinant(self, o):
+ return Variable(name_facet_jacobian_determinant())
+
@iname
def _dimension_iname(context, count):
@@ -341,3 +348,9 @@ def name_jacobian_determinant():
name = 'detjac'
define_jacobian_determinant(name)
return name
+
+
+def name_facet_jacobian_determinant():
+ name = 'fdetjac'
+ define_jacobian_determinant(name)
+ return name
\ No newline at end of file
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 1b65c0a33653917aa1065d86789fdbe7e86cb799..612eb77b1e4e4b23469ae3be02298ec5ed07c915 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -485,7 +485,7 @@ def generate_localoperator_kernels(formdata, data):
# In the case of matrix free operator evaluation we need jacobian apply methods
if get_option("matrix_free"):
from dune.perftool.pdelab.driver import is_linear
- if is_linear(form):
+ if is_linear(formdata.original_form):
# Numeical jacobian apply base class
base_class("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype), classtag="operator")
@@ -509,17 +509,10 @@ def generate_localoperator_kernels(formdata, data):
# Generate the necessary jacobian methods
if not get_option("numerical_jacobian"):
from ufl import derivative
- from ufl.algorithms import compute_form_data
- from dune.perftool.ufl.execution import (JacobianDeterminant,
- JacobianInverse,
- )
jacform = derivative(formdata.original_form, formdata.original_form.coefficients()[0])
- jacform = compute_form_data(jacform,
- do_apply_function_pullbacks=True,
- do_apply_geometry_lowering=True,
- do_apply_integral_scaling=True,
- preserve_geometry_types=(JacobianDeterminant, JacobianInverse),
- ).preprocessed_form
+
+ from dune.perftool.ufl.preprocess import preprocess_form
+ jacform = preprocess_form(jacform).preprocessed_form
with global_context(form_type="jacobian"):
for measure in set(i.integral_type() for i in jacform.integrals()):
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 5db88540ad6e001d6a7753414924e17ebc543b8d..91a25b3318a047ccfb85a2be19dd5311d80cad1f 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -52,10 +52,18 @@ def name_quadrature_position_in_cell(restriction):
return to_cell_coordinates(name_quadrature_position(), name_quadrature_point(), restriction)
-def name_quadrature_weight():
+def define_quadrature_weight(name):
qp = name_quadrature_point()
- return "{}.weight()".format(qp)
+ return quadrature_preamble(code="auto {} = {}.weight();".format(name, qp),
+ assignees=frozenset({name})
+ )
+
+def name_quadrature_weight():
+ name = 'weight'
+ temporary_variable(name, shape=())
+ define_quadrature_weight(name)
+ return name
def name_order():
diff --git a/python/dune/perftool/ufl/preprocess.py b/python/dune/perftool/ufl/preprocess.py
new file mode 100644
index 0000000000000000000000000000000000000000..540300bb3ee61abaa89b6b958df65e879af69e46
--- /dev/null
+++ b/python/dune/perftool/ufl/preprocess.py
@@ -0,0 +1,36 @@
+""" Preprocessing algorithms for UFL forms """
+
+
+def preprocess_form(form):
+ from ufl.classes import (FacetJacobianDeterminant,
+ FacetNormal,
+ JacobianDeterminant,
+ JacobianInverse,
+ )
+
+ from ufl.algorithms import compute_form_data
+ formdata = compute_form_data(form,
+ do_apply_function_pullbacks=True,
+ do_apply_geometry_lowering=True,
+ do_apply_integral_scaling=True,
+ preserve_geometry_types=(FacetJacobianDeterminant,
+ FacetNormal,
+ JacobianDeterminant,
+ JacobianInverse,
+ )
+ )
+
+ formdata.preprocessed_form = apply_default_transformations(formdata.preprocessed_form)
+
+ return formdata
+
+
+def apply_default_transformations(form):
+ from dune.perftool.ufl.transformations import transform_form
+ from dune.perftool.ufl.transformations.indexpushdown import pushdown_indexed
+ from dune.perftool.ufl.transformations.reindexing import reindexing
+
+ form = transform_form(form, pushdown_indexed)
+ form = transform_form(form, reindexing)
+
+ return form