diff --git a/CMakeLists.txt b/CMakeLists.txt
index 54d56ae9580085d059f0e1ca5be81ce6065419ee..e7ab965f0024f104cfa5a156df6d5eec1c1cfc67 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -28,6 +28,7 @@ add_subdirectory(python)
add_subdirectory(test)
add_subdirectory(bin)
+add_subdirectory(applications)
# finalize the dune project, e.g. generating config.h etc.
finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..22200ad351968fbeef9ddb19da07db127477f505
--- /dev/null
+++ b/applications/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(poisson_dg)
diff --git a/applications/poisson_dg/CMakeLists.txt b/applications/poisson_dg/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9bb13626e9b948e987c98188a1c9e61c1e0385e0
--- /dev/null
+++ b/applications/poisson_dg/CMakeLists.txt
@@ -0,0 +1,5 @@
+dune_add_formcompiler_system_test(UFLFILE poisson_dg.ufl
+ BASENAME app_poisson_dg
+ INIFILE poisson_dg.mini
+ NO_TESTS
+ )
diff --git a/applications/poisson_dg/poisson_dg.mini b/applications/poisson_dg/poisson_dg.mini
new file mode 100644
index 0000000000000000000000000000000000000000..e7a0eedf624a78760143f777449033c3fcafefb2
--- /dev/null
+++ b/applications/poisson_dg/poisson_dg.mini
@@ -0,0 +1,15 @@
+__name = app_poisson_dg_{__exec_suffix}
+__exec_suffix = deg{formcompiler.ufl_variants.degree}
+
+extension = 1.0 1.0 1.0
+cells = 16 16 16
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+sumfact = 1
+
+[formcompiler.ufl_variants]
+degree = 1, 2 | expand
diff --git a/applications/poisson_dg/poisson_dg.ufl b/applications/poisson_dg/poisson_dg.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..6a7e50f7e9f28f333304e713214f77516d00fe45
--- /dev/null
+++ b/applications/poisson_dg/poisson_dg.ufl
@@ -0,0 +1,28 @@
+cell = hexahedron
+
+x = SpatialCoordinate(cell)
+f = -6.
+g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
+
+V = FiniteElement("DG", cell, degree)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+gamma = 1.0
+theta = 1.0
+
+r = inner(grad(u), grad(v))*dx \
+ + inner(n, avg(grad(u)))*jump(v)*dS \
+ + gamma*jump(u)*jump(v)*dS \
+ + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+ - inner(n, grad(u))*v*ds \
+ + gamma*u*v*ds \
+ - theta*u*inner(grad(v), n)*ds \
+ - f*v*dx \
+ + theta*g*inner(grad(v), n)*ds \
+ - gamma*g*v*ds
+
+forms = [r]
diff --git a/cmake/modules/GeneratedSystemtests.cmake b/cmake/modules/GeneratedSystemtests.cmake
index 7f3e11367dd91900a33ae9df57e782d4cf46c036..2e3de48d33d9939aed33ffbec4d90035f50eb42b 100644
--- a/cmake/modules/GeneratedSystemtests.cmake
+++ b/cmake/modules/GeneratedSystemtests.cmake
@@ -1,9 +1,9 @@
-# System testing for generated executables. All ideas taken from dune-testtools.
+# System testing for generated executables. All ideas taken from dune-testtools.
#
function(dune_add_formcompiler_system_test)
# parse arguments
- set(OPTION DEBUG)
+ set(OPTION DEBUG NO_TESTS)
set(SINGLE INIFILE BASENAME SCRIPT UFLFILE)
set(MULTI CREATED_TARGETS)
cmake_parse_arguments(SYSTEMTEST "${OPTION}" "${SINGLE}" "${MULTI}" ${ARGN})
@@ -60,17 +60,19 @@ function(dune_add_formcompiler_system_test)
add_dependencies(${SYSTEMTEST_BASENAME} ${tname})
endif()
- # and have it depend on the metatarget build_tests
- add_dependencies(build_tests ${tname})
+ if(NOT ${SYSTEMTEST_NO_TESTS})
+ # and have it depend on the metatarget build_tests
+ add_dependencies(build_tests ${tname})
- _add_test(NAME ${tname}
- COMMAND ${CMAKE_BINARY_DIR}/dune-env ${SYSTEMTEST_SCRIPT}
- --exec ${tname}
- --ini "${CMAKE_CURRENT_BINARY_DIR}/${inifile}"
- --source ${CMAKE_CURRENT_SOURCE_DIR}
- )
+ _add_test(NAME ${tname}
+ COMMAND ${CMAKE_BINARY_DIR}/dune-env ${SYSTEMTEST_SCRIPT}
+ --exec ${tname}
+ --ini "${CMAKE_CURRENT_BINARY_DIR}/${inifile}"
+ --source ${CMAKE_CURRENT_SOURCE_DIR}
+ )
- set_tests_properties(${tname} PROPERTIES SKIP_RETURN_CODE 77)
- set_tests_properties(${tname} PROPERTIES TIMEOUT 20)
+ set_tests_properties(${tname} PROPERTIES SKIP_RETURN_CODE 77)
+ set_tests_properties(${tname} PROPERTIES TIMEOUT 20)
+ endif()
endforeach()
-endfunction()
\ No newline at end of file
+endfunction()
diff --git a/patches/ufl/0001-Give-conditional-an-ufl_id.patch b/patches/ufl/0001-Give-conditional-an-ufl_id.patch
new file mode 100644
index 0000000000000000000000000000000000000000..7080209276521bcb5ff181a5ee7f94afbf4d1bb1
--- /dev/null
+++ b/patches/ufl/0001-Give-conditional-an-ufl_id.patch
@@ -0,0 +1,26 @@
+From 8caa33075fba04cfc5bdcb47452a3c5ca400c0c5 Mon Sep 17 00:00:00 2001
+From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
+Date: Thu, 15 Dec 2016 23:30:03 +0100
+Subject: [PATCH] Give conditional an ufl_id
+
+---
+ ufl/conditional.py | 3 +++
+ 1 file changed, 3 insertions(+)
+
+diff --git a/ufl/conditional.py b/ufl/conditional.py
+index faffbb5..e48ee19 100644
+--- a/ufl/conditional.py
++++ b/ufl/conditional.py
+@@ -240,6 +240,9 @@ class Conditional(Operator):
+ def __repr__(self):
+ return "Conditional(%r, %r, %r)" % self.ufl_operands
+
++ def ufl_id(self):
++ return 0
++
+
+ #--- Specific functions higher level than a conditional ---
+
+--
+2.1.4
+
diff --git a/patches/vectorclass/0001-Correct-preprocessor-gurads-in-vectormath_exp.h.patch b/patches/vectorclass/0001-Correct-preprocessor-gurads-in-vectormath_exp.h.patch
new file mode 100644
index 0000000000000000000000000000000000000000..1b5d17a19ad4abab2c3680933be1a95f1c351e8a
--- /dev/null
+++ b/patches/vectorclass/0001-Correct-preprocessor-gurads-in-vectormath_exp.h.patch
@@ -0,0 +1,52 @@
+From d4c2e050005c06ca3d423ce45279a99550eff497 Mon Sep 17 00:00:00 2001
+From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
+Date: Mon, 19 Dec 2016 13:45:41 +0100
+Subject: [PATCH] Correct preprocessor gurads in vectormath_exp.h
+
+---
+ vectormath_exp.h | 8 ++++----
+ 1 file changed, 4 insertions(+), 4 deletions(-)
+
+diff --git a/vectormath_exp.h b/vectormath_exp.h
+index 49f0e09..729983a 100644
+--- a/vectormath_exp.h
++++ b/vectormath_exp.h
+@@ -1261,8 +1261,6 @@ static inline Vec16f square_cbrt(Vec16f const & x) {
+ return cbrt_f<Vec16f, Vec16ui, Vec16fb, 2> (x);
+ }
+
+-#endif // MAX_VECTOR_SIZE >= 512
+-
+ // Helper function for power function: insert special values of pow(x,y) when x=0:
+ // y<0 -> inf, y=0 -> 1, y>0 -> 0, y=nan -> nan
+ static inline Vec8d wm_pow_case_x0(Vec8db const & xiszero, Vec8d const & y, Vec8d const & z) {
+@@ -1274,6 +1272,8 @@ static inline Vec8d wm_pow_case_x0(Vec8db const & xiszero, Vec8d const & y, Vec8
+ #endif
+ }
+
++#endif // MAX_VECTOR_SIZE >= 512
++
+ static inline Vec4d wm_pow_case_x0(Vec4db const & xiszero, Vec4d const & y, Vec4d const & z) {
+ //#if defined __AVX512VL__ && defined ?
+ // const __m256i table = Vec4q(0x85858A00);
+@@ -1527,8 +1527,6 @@ inline Vec8d pow<float>(Vec8d const & x, float const & y) {
+ return pow_template_d<Vec8d, Vec8q, Vec8db>(x, (double)y);
+ }
+
+-#endif // MAX_VECTOR_SIZE >= 512
+-
+ // Helper function for power function: insert special values of pow(x,y) when x=0:
+ // y<0 -> inf, y=0 -> 1, y>0 -> 0, y=nan -> nan
+ static inline Vec16f wm_pow_case_x0(Vec16fb const & xiszero, Vec16f const & y, Vec16f const & z) {
+@@ -1540,6 +1538,8 @@ static inline Vec16f wm_pow_case_x0(Vec16fb const & xiszero, Vec16f const & y, V
+ #endif
+ }
+
++#endif // MAX_VECTOR_SIZE >= 512
++
+ static inline Vec8f wm_pow_case_x0(Vec8fb const & xiszero, Vec8f const & y, Vec8f const & z) {
+ return select(xiszero, select(y < 0.f, infinite_vec<Vec8f>(), select(y == 0.f, Vec8f(1.f), Vec8f(0.f))), z);
+ }
+--
+2.1.4
+
diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index 7e09f2fe96794c283ccfa61f06ce33af907e8006..1a3f725b4d21cd8d4f8e043082a6ae613d5dbf10 100644
--- a/python/dune/perftool/compile.py
+++ b/python/dune/perftool/compile.py
@@ -20,13 +20,22 @@ from dune.perftool.pdelab.localoperator import (generate_localoperator_basefile,
name_localoperator_file)
from dune.perftool.ufl.preprocess import preprocess_form
-import os.path
+from os.path import splitext, basename
# Disable loopy caching before we do anything else!
loopy.CACHING_ENABLED = False
+def type_guessing(val):
+ for t in [int, float]:
+ try:
+ return t(val)
+ except TypeError:
+ pass
+ return val
+
+
def read_ufl(uflfile):
"""Read uflfile file, extract and preprocess forms
@@ -41,13 +50,23 @@ def read_ufl(uflfile):
"""
# Read the given ufl file and execute it
uflcode = read_ufl_file(uflfile)
+
+ # Prepopulate a namespace with variation information
namespace = globals()
+ ini = get_option("ini_file")
+ if ini:
+ from dune.common.parametertree.parser import parse_ini_file
+ ini = parse_ini_file(ini)
+
+ for k, v in ini.get("formcompiler.ufl_variants", {}).items():
+ namespace[k] = type_guessing(v)
+
try:
exec("from dune.perftool.ufl.execution import *\n" + uflcode, namespace)
except:
- basename = os.path.splitext(os.path.basename(uflfile))[0]
- basename = "{}_debug".format(basename)
- pyname = "{}.py".format(basename)
+ name = splitext(basename(uflfile))[0]
+ name = "{}_debug".format(name)
+ pyname = "{}.py".format(name)
print(pyname)
pycode = "#!/usr/bin/env python\nfrom dune.perftool.ufl.execution import *\nset_level(DEBUG)\n" + uflcode
with file(pyname, "w") as f:
@@ -56,6 +75,11 @@ def read_ufl(uflfile):
# Extract and preprocess the forms
data = interpret_ufl_namespace(namespace)
+
+ # Enrich data by some additional objects we deem worth keeping
+ if get_option("exact_solution_expression"):
+ data.object_by_name[get_option("exact_solution_expression")] = namespace[get_option("exact_solution_expression")]
+
formdatas = []
forms = data.forms
for index, form in enumerate(forms):
diff --git a/python/dune/perftool/loopy/mangler.py b/python/dune/perftool/loopy/mangler.py
new file mode 100644
index 0000000000000000000000000000000000000000..c0172170f8ea3dbeeda9871a209795cd38eed7d4
--- /dev/null
+++ b/python/dune/perftool/loopy/mangler.py
@@ -0,0 +1,23 @@
+""" Function manglers for math functions in C++ """
+
+from dune.perftool.generation import (function_mangler,
+ include_file,
+ )
+
+from loopy import CallMangleInfo
+
+
+@function_mangler
+def dune_math_manglers(kernel, name, arg_dtypes):
+ if name == "exp":
+ include_file("dune/perftool/vectorclass/vectormath_exp.h", filetag="operatorfile")
+ return CallMangleInfo("exp",
+ arg_dtypes,
+ arg_dtypes,
+ )
+
+ if name == "sqrt":
+ return CallMangleInfo("sqrt",
+ arg_dtypes,
+ arg_dtypes,
+ )
diff --git a/python/dune/perftool/loopy/target.py b/python/dune/perftool/loopy/target.py
index 0d46b7e2d9c955fcb345b4f448158c732fc67c77..952d6e6429d230fed00c295e133d4ca4a083cda1 100644
--- a/python/dune/perftool/loopy/target.py
+++ b/python/dune/perftool/loopy/target.py
@@ -86,7 +86,9 @@ class DuneExpressionToCExpressionMapper(ExpressionToCExpressionMapper):
if self.codegen_state.vectorization_info:
# If this is vectorized we call the VCL function mul_add
return prim.Call(prim.Variable("mul_add"),
- (self.rec(expr.mul_op1), self.rec(expr.mul_op2), self.rec(expr.add_op)))
+ (self.rec(expr.mul_op1, type_context),
+ self.rec(expr.mul_op2, type_context),
+ self.rec(expr.add_op, type_context)))
else:
# Default implementation that discards the node in favor of the resp.
# additions and multiplications.
diff --git a/python/dune/perftool/loopy/transformations/collect_rotate.py b/python/dune/perftool/loopy/transformations/collect_rotate.py
index 424056de75c539eeef83ccba3c1a660932ec403d..61282fd4d422536b12cd5775dd0b561251a998ef 100644
--- a/python/dune/perftool/loopy/transformations/collect_rotate.py
+++ b/python/dune/perftool/loopy/transformations/collect_rotate.py
@@ -124,17 +124,40 @@ def collect_vector_data_rotate(knl):
)
)
elif isinstance(insn, lp.CInstruction):
- # Rip apart the code and change the assignee
- assignee, expression = insn.code.split("=")
-
- # TODO This is a bit whacky: It only works for scalar assignees
- # OTOH this code is on its way out anyway, because of CInstruction
- assignee = prim.Subscript(prim.Variable(arrname), (prim.Variable('rotate_index'),))
-
- code = "{} ={}".format(str(assignee), expression)
- new_insns.append(insn.copy(code=code,
- depends_on_is_final=True,
- ))
+ # This entire code path should go away as we either
+ # * switch CInstructions to implicit iname assignments (see https://github.com/inducer/loopy/issues/55)
+ # * switch to doing geometry stuff for sum factorization ourselves
+ if len(shape) == 0:
+ # Rip apart the code and change the assignee
+ assignee, expression = insn.code.split("=")
+ assignee = prim.Subscript(prim.Variable(arrname), (prim.Variable('rotate_index'),))
+ code = "{} ={}".format(str(assignee), expression)
+ new_insns.append(insn.copy(code=code,
+ depends_on_is_final=True,
+ ))
+ else:
+ # This is a *very* unfortunate code path
+
+ # Get inames to assign to the vector buffer
+ cinsn_inames = tuple("{}_assign_{}".format(quantity, i) for i in range(len(shape)))
+ domains = frozenset("{{ [{0}] : 0<={0}<{1} }}".format(iname, shape[i]) for i, iname in enumerate(cinsn_inames))
+ for dom in domains:
+ domain = parse_domains(dom, {})
+ knl = knl.copy(domains=knl.domains + domain)
+
+ # We keep the old writing instructions
+ new_insns.append(insn)
+
+ # and write a new one
+ cinsn_id = "{}_assign_id".format(quantity)
+ new_insns.append(lp.Assignment(prim.Subscript(prim.Variable(arrname), tuple(prim.Variable(i) for i in cinsn_inames) + (prim.Variable('rotate_index'),)),
+ prim.Subscript(prim.Variable(quantity), tuple(prim.Variable(i) for i in cinsn_inames)),
+ within_inames=common_inames.union(inames).union(frozenset(cinsn_inames)),
+ within_inames_is_final=True,
+ depends_on=frozenset({lp.match.Writes(quantity)}),
+ id=cinsn_id,
+ ))
+ all_writers.append(cinsn_id)
else:
raise NotImplementedError
elif quantity in knl.temporary_variables:
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 6bbde90cae5fffbb4cb34871340402f5b3ab1474..76f894bfe0a5cef91d547557173ea15b18a25c6a 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -16,6 +16,7 @@ from dune.perftool.pdelab.geometry import (dimension_iname,
name_jacobian_inverse_transposed,
name_unit_inner_normal,
name_unit_outer_normal,
+ to_global,
)
from dune.perftool.pdelab.index import (name_index,
)
@@ -23,10 +24,12 @@ from dune.perftool.pdelab.parameter import (cell_parameter_function,
intersection_parameter_function,
)
from dune.perftool.pdelab.quadrature import (pymbolic_quadrature_weight,
+ pymbolic_quadrature_position_in_cell,
quadrature_inames,
)
from dune.perftool.pdelab.spaces import (lfs_inames,
)
+from dune.perftool.pdelab.tensors import pymbolic_list_tensor
class PDELabInterface(object):
@@ -80,10 +83,20 @@ class PDELabInterface(object):
def cell_parameter_function(self, name, expr, restriction, cellwise_constant):
return cell_parameter_function(name, expr, restriction, cellwise_constant)
+ #
+ # Tensor expression related generator functions
+ #
+
+ def pymbolic_list_tensor(self, o, visitor):
+ return pymbolic_list_tensor(o, visitor)
+
#
# Geometry related generator functions
#
+ def pymbolic_spatial_coordinate(self, restriction):
+ return to_global(pymbolic_quadrature_position_in_cell(restriction))
+
def name_facet_jacobian_determinant(self):
return name_facet_jacobian_determinant()
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 39ef2e247a8206c78339d3e6c6dad0ec242d8c50..548d6c08ce88c198c86a135760b794023f8a9704 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -8,6 +8,7 @@ gained there.
"""
from dune.perftool.error import PerftoolCodegenError
from dune.perftool.generation import (generator_factory,
+ global_context,
include_file,
cached,
preamble,
@@ -429,12 +430,27 @@ def name_constraintscontainer(expr):
@preamble
def define_intersection_lambda(expression, name):
+ from dune.perftool.ufl.execution import Expression
+ from ufl.classes import Expr
if expression is None:
- return "auto {} = [&](const auto& x){{ return 0; }};".format(name)
- if expression.is_global:
- return "auto {} = [&](const auto& x){{ {} }};".format(name, expression.c_expr[0])
- else:
- return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr[0])
+ return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
+ elif isinstance(expression, Expression):
+ if expression.is_global:
+ return "auto {} = [&](const auto& x){{ {} }};".format(name, expression.c_expr[0])
+ else:
+ return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr[0])
+ elif isinstance(expression, Expr):
+ # Set up a visitor
+ with global_context(integral_type="exterior_facet", formdata=_driver_data["formdata"], driver=True):
+ from dune.perftool.ufl.visitor import UFL2LoopyVisitor
+ from dune.perftool.pdelab import PDELabInterface
+ visitor = UFL2LoopyVisitor(PDELabInterface(), "exterior_facet", {})
+ from pymbolic.mapper.c_code import CCodeMapper
+ ccm = CCodeMapper()
+ expr = visitor(expression)
+ return "auto {} = [&](const auto& x){{ return {}; }};".format(name, ccm(expr))
+
+ raise ValueError("Expression not understood")
def name_bctype_lambda(name, dirichlet):
@@ -748,13 +764,27 @@ def define_vector(name, formdata):
@preamble
def define_boundary_lambda(boundary, name):
+ from dune.perftool.ufl.execution import Expression
+ from ufl.classes import Expr
if boundary is None:
return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
- if boundary.is_global:
- return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr[0])
- else:
- return "auto {} = [&](const auto& e, const auto& x){{ {} }};".format(name, boundary.c_expr[0])
-
+ elif isinstance(boundary, Expression):
+ if boundary.is_global:
+ return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr[0])
+ else:
+ return "auto {} = [&](const auto& e, const auto& x){{ {} }};".format(name, boundary.c_expr[0])
+ elif isinstance(boundary, Expr):
+ # Set up a visitor
+ with global_context(integral_type="exterior_facet", formdata=_driver_data["formdata"], driver=True):
+ from dune.perftool.ufl.visitor import UFL2LoopyVisitor
+ from dune.perftool.pdelab import PDELabInterface
+ visitor = UFL2LoopyVisitor(PDELabInterface(), "exterior_facet", {})
+ from dune.perftool.loopy.target import numpy_to_cpp_dtype
+ return "auto {} = [&](const auto& x){{ return ({}){}; }};".format(name,
+ numpy_to_cpp_dtype("float64"),
+ visitor(boundary))
+
+ raise ValueError("Expression not understood")
def name_boundary_lambda(boundary, name):
define_boundary_lambda(boundary, name + "lambda")
@@ -790,8 +820,14 @@ def expression_splitter(expr, length):
if expr is None:
return (None,) * length
else:
- from dune.perftool.ufl.execution import split_expression
- return split_expression(expr)
+ from dune.perftool.ufl.execution import Expression, split_expression
+ from ufl.classes import ListTensor
+ if isinstance(expr, Expression):
+ return split_expression(expr)
+ elif isinstance(expr, ListTensor):
+ return expr.ufl_operands
+
+ raise TypeError("No idea how to split {} in components".format(type(expr)))
@cached
@@ -829,9 +865,17 @@ def name_solution_function(tree_path=()):
from dune.perftool.generation import get_global_context_value
expr = get_global_context_value("data").object_by_name[get_option("exact_solution_expression")]
- from dune.perftool.ufl.execution import split_expression
+ from dune.perftool.ufl.execution import Expression, split_expression
+ from ufl.classes import ListTensor
for i in tree_path:
- expr = split_expression(expr)[int(i)]
+ if isinstance(expr, Expression):
+ expr = split_expression(expr)[int(i)]
+ elif isinstance(expr, tuple):
+ expr = expr[int(i)]
+ elif isinstance(expr, ListTensor):
+ expr = expr.ufl_operands[int(i)]
+ else:
+ raise TypeError("No idea how to split {}".format(type(expr)))
from dune.perftool.generation import get_global_context_value
try:
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 1930ed923104649e6410b573fb878c8b9680c9d1..1abc8bf3627f539fed2b44439fa9201454dc0767 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -24,6 +24,7 @@ from pymbolic.primitives import Expression as PymbolicExpression
from loopy.match import Writes
import numpy as np
+import pymbolic.primitives as prim
@preamble
@@ -383,3 +384,23 @@ def name_facet_jacobian_determinant():
name = 'fdetjac'
get_backend(interface="detjac", selector=option_switch("constant_transformation_matrix"))(name)
return name
+
+
+def apply_to_global_transformation(name, local):
+ temporary_variable(name, shape=(world_dimension(),), shape_impl=("fv",))
+ geo = name_geometry()
+ code = "{} = {}.global({});".format(name,
+ geo,
+ local,
+ )
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(local)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(local))})
+ )
+
+def to_global(local):
+ assert isinstance(local, prim.Expression)
+ name = get_pymbolic_basename(local) + "_global"
+ apply_to_global_transformation(name, local)
+ return prim.Variable(name)
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index e214b5239eaa628832eaf8a3c30585d2db59bff2..0686f44c49001424ef3cb819c9df134243ab5b18 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -26,7 +26,11 @@ from dune.perftool.cgen.clazz import (AccessModifier,
ClassMember,
)
from dune.perftool.ufl.modified_terminals import Restriction
+
+import dune.perftool.loopy.mangler
+
from pymbolic.primitives import Variable
+import pymbolic.primitives as prim
from pytools import Record
import loopy as lp
@@ -247,7 +251,7 @@ def determine_accumulation_space(expr, number, measure):
)
-def boundary_predicates(expr, measure, subdomain_id):
+def boundary_predicates(expr, measure, subdomain_id, visitor):
predicates = frozenset([])
if subdomain_id not in ['everywhere', 'otherwise']:
@@ -273,16 +277,22 @@ def boundary_predicates(expr, measure, subdomain_id):
assert measure in subdomains
subdomain_data = subdomains[measure]
- # Determine the name of the parameter function
- name = get_global_context_value("data").object_names[id(subdomain_data)]
+ from ufl.classes import Expr
+ if isinstance(subdomain_data, Expr):
+ cond = visitor(subdomain_data)
+ else:
+ # Determine the name of the parameter function
+ cond = get_global_context_value("data").object_names[id(subdomain_data)]
+
+ # Trigger the generation of code for this thing in the parameter class
+ from ufl.checks import is_cellwise_constant
+ cellwise_constant = is_cellwise_constant(expr)
+ from dune.perftool.pdelab.parameter import intersection_parameter_function
+ intersection_parameter_function(cond, subdomain_data, cellwise_constant, t='int32')
- # Trigger the generation of code for this thing in the parameter class
- from ufl.checks import is_cellwise_constant
- cellwise_constant = is_cellwise_constant(expr)
- from dune.perftool.pdelab.parameter import intersection_parameter_function
- intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int32')
+ cond = prim.Variable(cond)
- predicates = predicates.union(['{} == {}'.format(name, subdomain_id)])
+ predicates = predicates.union([prim.Comparison(cond, '==', subdomain_id)])
return predicates
@@ -328,7 +338,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
from dune.perftool.pdelab.argument import name_accumulation_variable
accumvar = name_accumulation_variable((ansatz_lfs.get_restriction() + test_lfs.get_restriction()))
- predicates = boundary_predicates(accterm.term, measure, subdomain_id)
+ predicates = boundary_predicates(accterm.term, measure, subdomain_id, visitor)
rank = 1 if ansatz_lfs.lfs is None else 2
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 36ec2a44d54ec2a1558229c0196c8d5fbc1d8978..81a918a2ddf29789bfec40d440c30896f291d98f 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -105,7 +105,8 @@ def name_quadrature_points():
"""Name of vector storing quadrature points as class member"""
from dune.perftool.pdelab.geometry import local_dimension
dim = local_dimension()
- name = "qp_dim" + str(dim)
+ order = quadrature_order()
+ name = "qp_dim{}_order{}".format(dim, order)
shape = (name_quadrature_bound(), dim)
globalarg(name, shape=shape, dtype=numpy.float64, managed=False)
define_quadrature_points(name)
@@ -166,7 +167,8 @@ def name_quadrature_weights():
""""Name of vector storing quadrature weights as class member"""
from dune.perftool.pdelab.geometry import local_dimension
dim = local_dimension()
- name = "qw_dim" + str(dim)
+ order = quadrature_order()
+ name = "qw_dim{}_order{}".format(dim, order)
define_quadrature_weights(name)
fill_quadrature_weights_cache(name)
diff --git a/python/dune/perftool/pdelab/tensors.py b/python/dune/perftool/pdelab/tensors.py
new file mode 100644
index 0000000000000000000000000000000000000000..16cf6bb6a4d7bb3c1815ca372f3cdfdd53617167
--- /dev/null
+++ b/python/dune/perftool/pdelab/tensors.py
@@ -0,0 +1,27 @@
+""" Code generation for explicitly specified tensors """
+
+from dune.perftool.generation import (get_counted_variable,
+ kernel_cached,
+ instruction,
+ temporary_variable,
+ )
+
+import pymbolic.primitives as prim
+import numpy as np
+
+
+def define_list_tensor(name, expr, visitor):
+ for i, child in enumerate(expr.ufl_operands):
+ instruction(assignee=prim.Subscript(prim.Variable(name), (i,)),
+ expression=visitor.call(child))
+
+
+@kernel_cached
+def pymbolic_list_tensor(expr, visitor):
+ name = get_counted_variable("listtensor")
+ temporary_variable(name,
+ shape=expr.ufl_shape,
+ dtype=np.float64,
+ )
+ define_list_tensor(name, expr, visitor)
+ return prim.Variable(name)
diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py
index cba13702efbf8f8a7333dc919e48c8de87048232..1900e7b51987870a2a5baeadac09f4f01dfd2a22 100644
--- a/python/dune/perftool/sumfact/__init__.py
+++ b/python/dune/perftool/sumfact/__init__.py
@@ -1,5 +1,6 @@
from dune.perftool.sumfact.quadrature import (quadrature_inames,
quadrature_weight,
+ pymbolic_quadrature_position_in_cell,
)
from dune.perftool.sumfact.basis import (lfs_inames,
@@ -38,3 +39,7 @@ class SumFactInterface(PDELabInterface):
def pymbolic_quadrature_weight(self):
return quadrature_weight()
+
+ def pymbolic_spatial_coordinate(self, restriction):
+ from dune.perftool.pdelab.geometry import to_global
+ return to_global(pymbolic_quadrature_position_in_cell(restriction))
diff --git a/python/dune/perftool/sumfact/amatrix.py b/python/dune/perftool/sumfact/amatrix.py
index 3822ba43818301ea8a4240163cdefb0239f5ae7e..6cadb6193d8c22b954e0db80825b5ad1897198e4 100644
--- a/python/dune/perftool/sumfact/amatrix.py
+++ b/python/dune/perftool/sumfact/amatrix.py
@@ -70,10 +70,11 @@ class LargeAMatrix(ImmutableRecord):
@property
def name(self):
- name = "ThetaLarge{}{}_{}".format("face{}_".format(self.face) if self.face is not None else "",
- "T" if self.transpose else "",
- "_".join(tuple("d" if d else "" for d in self.derivative))
- )
+ name = "ThetaLarge{}{}_{}_qp{}".format("face{}_".format(self.face) if self.face is not None else "",
+ "T" if self.transpose else "",
+ "_".join(tuple("d" if d else "" for d in self.derivative)),
+ quadrature_points_per_direction(),
+ )
for i, d in enumerate(self.derivative):
define_theta(name, (self.rows, self.cols), self.transpose, d, face=self.face, additional_indices=(i,))
@@ -125,7 +126,8 @@ def define_oned_quadrature_weights(name):
def name_oned_quadrature_weights():
- name = "qw"
+ num = quadrature_points_per_direction()
+ name = "qw_num{}".format(num)
define_oned_quadrature_weights(name)
return name
@@ -139,7 +141,8 @@ def define_oned_quadrature_points(name):
def name_oned_quadrature_points():
- name = "qp"
+ num = quadrature_points_per_direction()
+ name = "qp_num{}".format(num)
define_oned_quadrature_points(name)
return name
@@ -182,12 +185,10 @@ def name_polynomials():
@preamble(kernel="operator")
-def sort_quadrature_points_weights():
+def sort_quadrature_points_weights(qp, qw):
range_field = lop_template_range_field()
domain_field = name_domain_field()
number_qp = quadrature_points_per_direction()
- qp = name_oned_quadrature_points()
- qw = name_oned_quadrature_weights()
include_file("dune/perftool/sumfact/onedquadrature.hh", filetag="operatorfile")
return "onedQuadraturePointsWeights<{}, {}, {}>({}, {});".format(range_field, domain_field, number_qp, qp, qw)
@@ -219,9 +220,10 @@ def polynomial_lookup_mangler(target, func, dtypes):
def define_theta(name, shape, transpose, derivative, face=None, additional_indices=()):
- sort_quadrature_points_weights()
- polynomials = name_polynomials()
qp = name_oned_quadrature_points()
+ qw = name_oned_quadrature_weights()
+ sort_quadrature_points_weights(qp, qw)
+ polynomials = name_polynomials()
dim_tags = "f,f"
if additional_indices:
@@ -255,10 +257,11 @@ def define_theta(name, shape, transpose, derivative, face=None, additional_indic
def name_theta(transpose=False, derivative=False, face=None):
- name = "{}{}Theta{}".format("face{}_".format(face) if face is not None else "",
- "d" if derivative else "",
- "T" if transpose else "",
- )
+ name = "{}{}Theta{}_qp_{}".format("face{}_".format(face) if face is not None else "",
+ "d" if derivative else "",
+ "T" if transpose else "",
+ quadrature_points_per_direction(),
+ )
shape = [quadrature_points_per_direction(), basis_functions_per_direction()]
if face is not None:
diff --git a/python/dune/perftool/sumfact/basis.py b/python/dune/perftool/sumfact/basis.py
index d3452f31a0391493ce21a571a9c06c98fe4a6625..86a9b354efa1d1b6bc07d90edc45b8c531ff577c 100644
--- a/python/dune/perftool/sumfact/basis.py
+++ b/python/dune/perftool/sumfact/basis.py
@@ -68,7 +68,6 @@ def pymbolic_trialfunction_gradient(element, restriction, component, visitor):
dim = world_dimension()
buffers = []
- insn_dep = None
for i in range(dim):
# Construct the matrix sequence for this sum factorization
a_matrices = construct_amatrix_sequence(derivative=i,
@@ -90,8 +89,7 @@ def pymbolic_trialfunction_gradient(element, restriction, component, visitor):
).get_temporary(shape=shape,
name=inp,
)
- if insn_dep is None:
- insn_dep = frozenset({Writes(inp)})
+ insn_dep = frozenset({Writes(inp)})
if get_option('fastdg'):
# Name of direct input, shape and globalarg is set in sum_factorization_kernel
@@ -110,6 +108,7 @@ def pymbolic_trialfunction_gradient(element, restriction, component, visitor):
preferred_position=i,
insn_dep=insn_dep,
restriction=restriction,
+ outshape=tuple(mat.rows for mat in a_matrices if mat.face is None),
direct_input=direct_input,
)
buffers.append(var)
@@ -161,7 +160,6 @@ def pymbolic_trialfunction(element, restriction, component, visitor):
name=inp,
)
- # TODO: fastdg and vectorization
if get_option('fastdg'):
# Name of direct input, shape and globalarg is set in sum_factorization_kernel
direct_input = name_coefficientcontainer(restriction)
@@ -177,7 +175,7 @@ def pymbolic_trialfunction(element, restriction, component, visitor):
1,
preferred_position=None,
insn_dep=frozenset({Writes(inp)}),
- outshape=tuple(mat.rows for mat in a_matrices if mat.rows != 1),
+ outshape=tuple(mat.rows for mat in a_matrices if mat.face is None),
restriction=restriction,
direct_input=direct_input,
)
diff --git a/python/dune/perftool/sumfact/sumfact.py b/python/dune/perftool/sumfact/sumfact.py
index e135fb75a1a8650a088852f95d67d1da03a4b77e..2429da4ef3d688ebe651f7d4de140ca76a8a779a 100644
--- a/python/dune/perftool/sumfact/sumfact.py
+++ b/python/dune/perftool/sumfact/sumfact.py
@@ -144,7 +144,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
a_matrices, buf, inp, index, padding = get_vectorization_info(a_matrices, (accterm.argument.restriction, restriction))
# Initialize a base storage for this buffer and get a temporay pointing to it
- shape = tuple(mat.cols for mat in a_matrices if mat.cols != 1)
+ shape = tuple(mat.cols for mat in a_matrices if mat.face is None)
dim_tags = ",".join(['f'] * local_dimension())
if index is not None:
shape = shape + (4,)
@@ -434,7 +434,9 @@ def sum_factorization_kernel(a_matrices,
})
if outshape is None:
- outshape = tuple(mat.rows for mat in a_matrices if mat.rows != 1)
+ assert stage == 3
+ outshape = tuple(mat.rows for mat in a_matrices)
+
dim_tags = ",".join(['f'] * len(outshape))
if next(iter(a_matrices)).vectorized:
diff --git a/python/dune/perftool/sumfact/switch.py b/python/dune/perftool/sumfact/switch.py
index d2361c3b052b0ed51fcc584eccfa74ba0f97ce42..6b81cfaa8a0d451b96b0f4cd96965077c5cdfabb 100644
--- a/python/dune/perftool/sumfact/switch.py
+++ b/python/dune/perftool/sumfact/switch.py
@@ -56,6 +56,11 @@ def decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
if not get_option("diagonal_transformation_matrix"):
return True
+ # The PDELab machineries visit-once policy combined with Yasp avoids any visits
+ # with facemod_s being True
+ if facemod_s:
+ return False
+
# A codim1 entity can never be on the upper resp. lower side of the ref element
# in both inside and outside cell in a YaspGrid
if facemod_n == facemod_s:
@@ -84,7 +89,7 @@ def generate_exterior_facet_switch():
for facedir_s in range(dim):
for facemod_s in range(2):
- block.append(" case {}: {}({}); break;".format(dim * facedir_s + facemod_s,
+ block.append(" case {}: {}({}); break;".format(2 * facedir_s + facemod_s,
get_kernel_name(facedir_s=facedir_s,
facemod_s=facemod_s,
),
@@ -114,7 +119,7 @@ def generate_interior_facet_switch():
for facedir_n in range(dim):
for facemod_n in range(2):
if decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
- block.append(" case {}: {}({}); break;".format((dim * facedir_s + facemod_s) * (2 * dim) + dim * facedir_n + facemod_n,
+ block.append(" case {}: {}({}); break;".format((2 * facedir_s + facemod_s) * (2 * dim) + 2 * facedir_n + facemod_n,
get_kernel_name(facedir_s=facedir_s,
facemod_s=facemod_s,
facedir_n=facedir_n,
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 57d337fbe541beb5de6afe429c5a3930dccf962e..93769fe2c0ed2bcb8f61db1ff9180c2006bd6346 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -145,13 +145,14 @@ class FiniteElement(ufl.FiniteElement):
self.dirichlet_constraints = kwargs.pop('dirichlet_constraints', 'return true;')
if isinstance(self.dirichlet_constraints, str):
self.dirichlet_constraints = Expression(self.dirichlet_constraints)
- assert isinstance(self.dirichlet_constraints, Expression)
+ from ufl.classes import Expr
+ assert isinstance(self.dirichlet_constraints, (Expression, Expr))
# Get dirichlet_constraints and convert it to Expression if necessary!
self.dirichlet_expression = kwargs.pop('dirichlet_expression', 'return 0.0;')
if isinstance(self.dirichlet_expression, str):
self.dirichlet_expression = Expression(self.dirichlet_expression)
- assert isinstance(self.dirichlet_expression, Expression)
+ assert isinstance(self.dirichlet_expression, (Expression, Expr))
# Initialize the original finite element from ufl
ufl.FiniteElement.__init__(self, *args, **kwargs)
@@ -164,13 +165,14 @@ class VectorElement(ufl.VectorElement):
self.dirichlet_constraints = kwargs.pop('dirichlet_constraints', 'return true;')
if isinstance(self.dirichlet_constraints, str):
self.dirichlet_constraints = Expression(self.dirichlet_constraints)
- assert isinstance(self.dirichlet_constraints, Expression)
+ from ufl.classes import Expr
+ assert isinstance(self.dirichlet_constraints, (Expression, Expr))
# Get dirichlet_constraints and convert it to Expression if necessary!
self.dirichlet_expression = kwargs.pop('dirichlet_expression', 'return 0.0;')
if isinstance(self.dirichlet_expression, str):
self.dirichlet_expression = Expression(self.dirichlet_expression)
- assert isinstance(self.dirichlet_expression, Expression)
+ assert isinstance(self.dirichlet_expression, (Expression, Expr))
# Initialize the original finite element from ufl
ufl.VectorElement.__init__(self, *args, **kwargs)
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index c1c9c8659af11007ed2cd904a6a601920cf31aab..14e83916eaf4223bf13ff1dfb5574b8478581de5 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -33,6 +33,8 @@ from ufl.classes import (FixedIndex,
JacobianDeterminant,
)
+import pymbolic.primitives as prim
+
class UFL2LoopyVisitor(ModifiedTerminalTracker):
def __init__(self, interface, measure, dimension_indices):
@@ -49,6 +51,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
def __call__(self, o):
# Reset some state variables that are reinitialized for each accumulation term
self.indices = None
+ self._indices_backup = []
self.inames = ()
return self.call(o)
@@ -149,6 +152,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
# This is necessary in some places where the locality of the tree transformation
# is not given (easiest example: jacobian_inverse handler, complex example:
# sum factorziation handler for trial function gradient)
+ self._indices_backup.append(self.indices)
self.indices = self.call(o.ufl_operands[1])
# Handle the aggregate!
@@ -156,10 +160,11 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
# self.indices being None means they are already handled during recursion!
if self.indices is None:
+ self.indices = self._indices_backup.pop()
return aggr
else:
indices = self.indices
- self.indices = None
+ self.indices = self._indices_backup.pop()
return maybe_wrap_subscript(aggr, indices)
def index_sum(self, o, additional_inames=()):
@@ -208,6 +213,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
def index(self, o):
return self._index_or_fixed_index(o)
+ def list_tensor(self, o):
+ return self.interface.pymbolic_list_tensor(o, self)
+
#
# Handlers for arithmetic operators and functions
# Those handlers would be valid in any code going from UFL to pymbolic
@@ -235,12 +243,80 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
if isinstance(o.ufl_operands[0], JacobianDeterminant):
return self.call(o.ufl_operands[0])
else:
- return Call('abs', self.call(o.ufl_operands[0]))
+ return Call(Variable('abs'), (self.call(o.ufl_operands[0]),))
+
+ def exp(self, o):
+ return Call(Variable('exp'), (self.call(o.ufl_operands[0]),))
+
+ def sqrt(self, o):
+ return Call(Variable('sqrt'), (self.call(o.ufl_operands[0]),))
+
+ def power(self, o):
+ from ufl.constantvalue import IntValue
+ if isinstance(o.ufl_operands[1], IntValue):
+ if o.ufl_operands[1].value() == 2:
+ p = self.call(o.ufl_operands[0])
+ return Product((p, p))
+
+ raise NotImplementedError("Power function not really implemented")
+
+ #
+ # Handler for conditionals, use pymbolic base implementation
+ #
+
+ def conditional(self, o):
+ return prim.If(*tuple(self.call(op) for op in o.ufl_operands))
+
+ def eq(self, o):
+ return prim.Comparison(self.call(o.ufl_operands[0]),
+ "==",
+ right = self.call(o.ufl_operands[1]))
+
+ def ge(self, o):
+ return prim.Comparison(self.call(o.ufl_operands[0]),
+ ">=",
+ right = self.call(o.ufl_operands[1]))
+
+ def gt(self, o):
+ return prim.Comparison(self.call(o.ufl_operands[0]),
+ ">",
+ right = self.call(o.ufl_operands[1]))
+
+ def le(self, o):
+ return prim.Comparison(self.call(o.ufl_operands[0]),
+ "<=",
+ right = self.call(o.ufl_operands[1]))
+
+ def lt(self, o):
+ return prim.Comparison(self.call(o.ufl_operands[0]),
+ "<",
+ right = self.call(o.ufl_operands[1]))
+
+ def ne(self, o):
+ return prim.Comparison(self.call(o.ufl_operands[0]),
+ "!=",
+ right = self.call(o.ufl_operands[1]))
+
+ def and_condition(self, o):
+ return prim.LogicalAnd((self.call(o.ufl_operands[0]), self.call(o.ufl_operands[1])))
+
+ def or_condition(self, o):
+ return prim.LogicalOr((self.call(o.ufl_operands[0]), self.call(o.ufl_operands[1])))
+
+ def not_condition(self, o):
+ return prim.LogicalNot(self.call(o.ufl_operands[0]))
#
# Handlers for geometric quantities
#
+ def spatial_coordinate(self, o):
+ # If this is called from the driver, we just want to return x
+ if get_global_context_value("driver", False):
+ return prim.Variable("x")
+ else:
+ return self.interface.pymbolic_spatial_coordinate(self.restriction)
+
def facet_normal(self, o):
# The normal must be restricted to be well-defined
assert self.restriction is not Restriction.NONE
diff --git a/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl b/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl
index 6bda534f16236d66b351c0ac1c0671a41f03f9de..d2b9134176d9533068e52d995f17bf1ad78a3523 100644
--- a/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_1d_cg_interval.ufl
@@ -1,10 +1,12 @@
cell = "interval"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", cell=cell)
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2
+g = exp(-1.*c)
+f = 2*(1.-2*c)*g
V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
u = TrialFunction(V)
v = TestFunction(V)
-forms = [(inner(grad(u), grad(v)) - f*v)*dx]
\ No newline at end of file
+forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl b/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
index 907f19841d51695f3944617ff162b33fce7b5053..1b3189637a4481f9147274286fe974255e02e000 100644
--- a/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
@@ -1,7 +1,9 @@
cell = "interval"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2
+g = exp(-1.*c)
+f = 2*(1.-2*c)*g
V = FiniteElement("DG", cell, 1)
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini
index 8f5e193a613a7caeb1a9fc2796e9d88b0f5c9b22..ccbd4d0af1011db9fbd5b9107e28a743ef9057ab 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini
+++ b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini
@@ -9,5 +9,4 @@ extension = vtu
[formcompiler]
exact_solution_expression = g
-compare_dofs = 1e-14
compare_l2errorsquared = 1e-6
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl
index 5a2c570583a2fe200e91a1dbf20b0dde699ec156..74d3f2b1485e4d971e13a0a0601213f9b9b0dd60 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.ufl
@@ -1,7 +1,9 @@
cell = "quadrilateral"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", cell=cell)
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 2*(2.-2*c)*g
V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
u = TrialFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_cg_triangle.ufl b/test/poisson/dimension-grid-variations/poisson_2d_cg_triangle.ufl
index 902a7f2a36f739b0088ded8b08d12f3d906832da..2bc0536e848da384b713db8661d11ef503bb4b57 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_cg_triangle.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_cg_triangle.ufl
@@ -1,5 +1,8 @@
-f = Expression("return -2.0*x.size();")
-g = Expression("return x.two_norm2();")
+cell = triangle
+
+x = SpatialCoordinate(cell)
+f = -4.
+g = x[0]*x[0] + x[1]*x[1]
V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g)
u = TrialFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
index 1ca65e1eaf16d4ff71444bd71626c277b82328c9..3551ab6e77b2643c24c65bd4d0ea036d877a3be2 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
@@ -1,7 +1,9 @@
cell = "quadrilateral"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 2*(2.-2*c)*g
V = FiniteElement("DG", cell, 1)
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl b/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl
index bf606f9cdb7b3726b457688db4921c4aacfa99ab..54b9a2946ffc173b1173ffe911587fe6fed63b6e 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl
@@ -1,7 +1,9 @@
-f = Expression("return -2.0*x.size();")
-g = Expression("return x.two_norm2();", on_intersection=True)
-
cell = triangle
+
+x = SpatialCoordinate(cell)
+f = -4.
+g = x[0]*x[0] + x[1]*x[1]
+
V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl
index 4cad67e31f354a7eb0c4c6d25b95713f31067612..e9ff908ebe67d357db20509d42f03b46afd8d90a 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.ufl
@@ -1,7 +1,8 @@
cell = "hexahedron"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", cell=cell)
+x = SpatialCoordinate(cell)
+f = -6.
+g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
u = TrialFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl
index d8f43c0941267480a07eba857b93d635571af43d..27ee0fd5d63f4398cf1a6d92f750f07634b73546 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_cg_tetrahedron.ufl
@@ -1,7 +1,8 @@
cell = "tetrahedron"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", cell=cell)
+x = SpatialCoordinate(cell)
+f = -6.
+g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
V = FiniteElement("CG", "tetrahedron", 1, dirichlet_expression=g)
u = TrialFunction(V)
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
index bf5cc370a5651c1a408d2f1b6687d4c3255c0172..0ac2968b5fa7155dd4920774fe8583489ba3349a 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
@@ -1,7 +1,9 @@
cell = "hexahedron"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
+x = SpatialCoordinate(cell)
+c = (0.5 - x[0])**2 + (0.5 - x[1])**2 + (0.5 - x[2])**2
+g = exp(-1.*c)
+f = 2*(3.-2*c)*g
V = FiniteElement("DG", cell, 1)
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
index 4c73c7debf2ac3a5f81941d2e6d12cb5ee4ff4fa..9ab31b81e0deda97b098158614ac95998a162386 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
@@ -1,7 +1,8 @@
cell = "tetrahedron"
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
+x = SpatialCoordinate(cell)
+f = -6.
+g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
V = FiniteElement("DG", cell, 1)
diff --git a/test/poisson/poisson.ufl b/test/poisson/poisson.ufl
index 8c62dd3a24d2f64a3aec350c2096fd65967498d8..73515675141be1bd95dfd4d3d65d7c7195bb5bb3 100644
--- a/test/poisson/poisson.ufl
+++ b/test/poisson/poisson.ufl
@@ -1,8 +1,14 @@
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
+cell = triangle
-V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g)
+x = SpatialCoordinate(cell)
+
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 4*(1.-c)*g
+
+V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
u = TrialFunction(V)
v = TestFunction(V)
+
forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/poisson/poisson_dg.ufl b/test/poisson/poisson_dg.ufl
index 0ee5a857931d2ceea62b8cf5e52c5aed2b696851..06f3eb0f953695ccb9bf948fa947080901956e97 100644
--- a/test/poisson/poisson_dg.ufl
+++ b/test/poisson/poisson_dg.ufl
@@ -1,7 +1,10 @@
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", on_intersection=True)
-
cell = triangle
+
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 4*(1.-c)*g
+
V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
diff --git a/test/poisson/poisson_dg_neumann.ufl b/test/poisson/poisson_dg_neumann.ufl
index eeda053a2f2df24a59f6d88a88e075c69e95c754..6d1bd594e43ee40fbd40ce50f16052c5376c724f 100644
--- a/test/poisson/poisson_dg_neumann.ufl
+++ b/test/poisson/poisson_dg_neumann.ufl
@@ -1,9 +1,13 @@
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", on_intersection=True)
-j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return -2.*s*(x[1]-0.5)*std::exp(-1.*c.two_norm2());", on_intersection=True)
-bctype = Expression("if ((x[1]<1e-8) || (x[1]>1.-1e-8)) return 0; else return 1;", on_intersection=True)
-
cell = triangle
+
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 4*(1.-c)*g
+sgn = conditional(x[1] > 0.5, 1., -1.)
+j = -2.*sgn*(x[1]-0.5)*g
+bctype = conditional(Or(x[1]<1e-8, x[1]>1.-1e-8), 0, 1)
+
V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
diff --git a/test/poisson/poisson_dg_quadrilateral.ufl b/test/poisson/poisson_dg_quadrilateral.ufl
index e5372870c524abdc02bb9eb061dc7a8b13b6d197..f92ed970768ac6d9585e72d333dacb1884c24bf6 100644
--- a/test/poisson/poisson_dg_quadrilateral.ufl
+++ b/test/poisson/poisson_dg_quadrilateral.ufl
@@ -1,7 +1,9 @@
cell = "quadrilateral"
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", on_intersection=True, cell=cell)
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 4*(1.-c)*g
V = FiniteElement("DG", cell, 1)
diff --git a/test/poisson/poisson_matrix_free.mini b/test/poisson/poisson_matrix_free.mini
index 491f9c2a22b5c5d8360d8644a0dd6b5a41f99f6b..364a2723edf7ae4901370a3f68089685cafbb322 100644
--- a/test/poisson/poisson_matrix_free.mini
+++ b/test/poisson/poisson_matrix_free.mini
@@ -1,5 +1,4 @@
-__name = poisson_matrix_free_{__exec_suffix}
-__exec_suffix = numdiff, symdiff | expand num
+__name = poisson_matrix_free
lowerleft = 0.0 0.0
upperright = 1.0 1.0
@@ -12,7 +11,6 @@ reference = poisson_ref
extension = vtu
[formcompiler]
-numerical_jacobian = 1, 0 | expand num
matrix_free = 1
exact_solution_expression = g
compare_l2errorsquared = 1e-6
diff --git a/test/poisson/poisson_neumann.ufl b/test/poisson/poisson_neumann.ufl
index e86bb53b86f06c190736be2ab7b5fc0ea92e2792..f1a46662b51945393cba6f4576770bf8e372a36c 100644
--- a/test/poisson/poisson_neumann.ufl
+++ b/test/poisson/poisson_neumann.ufl
@@ -1,7 +1,13 @@
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
-j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return -2.*s*(x[1]-0.5)*std::exp(-1.*c.two_norm2());", on_intersection=True)
-bctype = Expression("if ((x[1]<1e-8) || (x[1]>1.-1e-8)) return 0; else return 1;", on_intersection=True)
+cell = triangle
+x = SpatialCoordinate(cell)
+
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 4*(1.-c)*g
+sgn = conditional(x[1] > 0.5, 1., -1.)
+j = -2.*sgn*(x[1]-0.5)*g
+
+bctype = conditional(Or(x[1]<1e-8, x[1]>1.-1e-8), 0, 1)
V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g, dirichlet_constraints=bctype)
u = TrialFunction(V)
diff --git a/test/stokes/stokes.ufl b/test/stokes/stokes.ufl
index da632caf9ecbf0f93eb5aa1c70cd5df17aba01b6..595d3ba80a40d99ae8a6d390b3c2a0cd2ee3a394 100644
--- a/test/stokes/stokes.ufl
+++ b/test/stokes/stokes.ufl
@@ -1,10 +1,11 @@
-v_bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
-g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"))
-g_p = Expression("8*(1.-x[0])")
+cell = triangle
-g = g_v * g_p
+x = SpatialCoordinate(cell)
+v_bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
+g_v = as_vector((4.*x[1]*(1.-x[1]), 0.0))
+g_p = 8.*(1.-x[0])
+g = (g_v, g_p)
-cell = triangle
P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
P1 = FiniteElement("Lagrange", cell, 1)
TH = P2 * P1
diff --git a/test/stokes/stokes_dg.ufl b/test/stokes/stokes_dg.ufl
index b63893f0fe98622e1b24019e930026c615074010..f1995c8b2f928d34a1965b62f3fee391dec4caa3 100644
--- a/test/stokes/stokes_dg.ufl
+++ b/test/stokes/stokes_dg.ufl
@@ -1,9 +1,11 @@
-g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"), on_intersection=True)
-g_p = Expression("8*(1.-x[0])")
-g = g_v * g_p
-bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
-
cell = triangle
+
+x = SpatialCoordinate(cell)
+g_v = as_vector((4*x[1]*(1.-x[1]), 0.0))
+g_p = 8*(1.-x[0])
+g = (g_v, g_p)
+bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
+
P2 = VectorElement("DG", cell, 2)
P1 = FiniteElement("DG", cell, 1)
TH = P2 * P1
diff --git a/test/sumfact/poisson/CMakeLists.txt b/test/sumfact/poisson/CMakeLists.txt
index 1dc2483f40d0e3a4af5542c2dc2a12dabbce53f4..889b29853ba92c628f91cfbc6a19523189d69fad 100644
--- a/test/sumfact/poisson/CMakeLists.txt
+++ b/test/sumfact/poisson/CMakeLists.txt
@@ -1,42 +1,35 @@
-# Poisson CG 2D
-dune_add_formcompiler_system_test(UFLFILE poisson_2d_order1.ufl
- BASENAME sumfact_poisson_2d_order1
- INIFILE poisson_2d_order1.mini
+# 1. Poisson Test Case in 2D
+dune_add_formcompiler_system_test(UFLFILE poisson_2d.ufl
+ BASENAME sumfact_poisson_2d
+ INIFILE poisson_2d.mini
)
-dune_add_formcompiler_system_test(UFLFILE poisson_2d_order2.ufl
- BASENAME sumfact_poisson_2d_order2
- INIFILE poisson_2d_order2.mini
+# 2. Poisson Test case in 3D
+dune_add_formcompiler_system_test(UFLFILE poisson_3d.ufl
+ BASENAME sumfact_poisson_3d
+ INIFILE poisson_3d.mini
)
-
-# Poisson CG 3D
-dune_add_formcompiler_system_test(UFLFILE poisson_3d_order1.ufl
- BASENAME sumfact_poisson_3d_order1
- INIFILE poisson_3d_order1.mini
- )
-dune_add_formcompiler_system_test(UFLFILE poisson_3d_order2.ufl
- BASENAME sumfact_poisson_3d_order2
- INIFILE poisson_3d_order2.mini
- )
-
-
# Operator counting test
dune_add_formcompiler_system_test(UFLFILE opcount_poisson_2d_order2.ufl
BASENAME opcount_sumfact_poisson_2d_order2
INIFILE opcount_poisson_2d_order2.mini
)
-
-# Poisson DG
-dune_add_formcompiler_system_test(UFLFILE poisson_dg.ufl
- BASENAME sumfact_poisson_dg
- INIFILE poisson_dg.mini
+# 3. Poisson Test Case: DG in 2D
+dune_add_formcompiler_system_test(UFLFILE poisson_dg_2d.ufl
+ BASENAME sumfact_poisson_dg_2d
+ INIFILE poisson_dg_2d.mini
)
+# 4. Poisson Test Case: DG in 3D
+dune_add_formcompiler_system_test(UFLFILE poisson_dg_3d.ufl
+ BASENAME sumfact_poisson_dg_3d
+ INIFILE poisson_dg_3d.mini
+ )
# Poisson DG using FastDGGridOperator
-dune_add_formcompiler_system_test(UFLFILE poisson_dg.ufl
+dune_add_formcompiler_system_test(UFLFILE poisson_dg_2d.ufl
BASENAME sumfact_poisson_fastdg
INIFILE poisson_fastdg.mini
)
diff --git a/test/sumfact/poisson/poisson_2d_order1.mini b/test/sumfact/poisson/poisson_2d.mini
similarity index 63%
rename from test/sumfact/poisson/poisson_2d_order1.mini
rename to test/sumfact/poisson/poisson_2d.mini
index 85e91b2135c564ac9ab743fec5502040b2c34a65..6f70b71dae74402374423f4fc7518f92031300a6 100644
--- a/test/sumfact/poisson/poisson_2d_order1.mini
+++ b/test/sumfact/poisson/poisson_2d.mini
@@ -1,8 +1,9 @@
-__name = sumfact_poisson_2d_order1_{__exec_suffix}
-__exec_suffix = {diff_suffix}_{vec_suffix}
+__name = sumfact_poisson_2d_{__exec_suffix}
+__exec_suffix = {deg_suffix}_{diff_suffix}_{vec_suffix}
diff_suffix = numdiff, symdiff | expand num
vec_suffix = vec, nonvec | expand vec
+deg_suffix = deg{formcompiler.ufl_variants.degree}
cells = 8 8
extension = 1. 1.
@@ -18,3 +19,6 @@ exact_solution_expression = g
compare_l2errorsquared = 1e-4
sumfact = 1
vectorize_quad = 1, 0 | expand vec
+
+[formcompiler.ufl_variants]
+degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_2d.ufl b/test/sumfact/poisson/poisson_2d.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..a2c74a593f694bfedfbe2cc947a581c28ce0c07e
--- /dev/null
+++ b/test/sumfact/poisson/poisson_2d.ufl
@@ -0,0 +1,12 @@
+cell = "quadrilateral"
+
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 2*(2.-2*c)*g
+
+V = FiniteElement("CG", cell, degree, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/sumfact/poisson/poisson_2d_order1.ufl b/test/sumfact/poisson/poisson_2d_order1.ufl
deleted file mode 100644
index 576c603e5c5f8b8d00bad8ea4cdcb6ec19ebd332..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_2d_order1.ufl
+++ /dev/null
@@ -1,10 +0,0 @@
-cell = "quadrilateral"
-
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", cell=cell)
-
-V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
-u = TrialFunction(V)
-v = TestFunction(V)
-
-forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/sumfact/poisson/poisson_2d_order2.mini b/test/sumfact/poisson/poisson_2d_order2.mini
deleted file mode 100644
index 5cdd15ff4bf39a8839ec1b1aa74f623ff4ad5585..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_2d_order2.mini
+++ /dev/null
@@ -1,20 +0,0 @@
-__name = sumfact_poisson_2d_order2_{__exec_suffix}
-__exec_suffix = {diff_suffix}_{vec_suffix}
-
-diff_suffix = numdiff, symdiff | expand num
-vec_suffix = vec, nonvec | expand vec
-
-cells = 8 8
-extension = 1. 1.
-
-[wrapper.vtkcompare]
-name = {__name}
-reference = poisson_ref
-extension = vtu
-
-[formcompiler]
-numerical_jacobian = 1, 0 | expand num
-exact_solution_expression = g
-compare_l2errorsquared = 1e-8
-sumfact = 1
-vectorize_quad = 1, 0 | expand vec
diff --git a/test/sumfact/poisson/poisson_2d_order2.ufl b/test/sumfact/poisson/poisson_2d_order2.ufl
deleted file mode 100644
index 702ba911668471600d003213c66745fb865922f0..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_2d_order2.ufl
+++ /dev/null
@@ -1,10 +0,0 @@
-cell = "quadrilateral"
-
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", cell=cell)
-
-V = FiniteElement("CG", cell, 2, dirichlet_expression=g)
-u = TrialFunction(V)
-v = TestFunction(V)
-
-forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/sumfact/poisson/poisson_3d_order1.mini b/test/sumfact/poisson/poisson_3d.mini
similarity index 62%
rename from test/sumfact/poisson/poisson_3d_order1.mini
rename to test/sumfact/poisson/poisson_3d.mini
index d794e0dca0f7f6b2147b22893a98339e33e3ee52..a5e2d647fe524bdb2ecd8855869ca1bb2c8e7bfa 100644
--- a/test/sumfact/poisson/poisson_3d_order1.mini
+++ b/test/sumfact/poisson/poisson_3d.mini
@@ -1,9 +1,10 @@
-__name = sumfact_poisson_3d_order1_{__exec_suffix}
-__exec_suffix = {diff_suffix}_{vecq_suffix}_{vecg_suffix}
+__name = sumfact_poisson_3d_{__exec_suffix}
+__exec_suffix = {deg_suffix}_{diff_suffix}_{vecq_suffix}_{vecg_suffix}
diff_suffix = numdiff, symdiff | expand num
vecq_suffix = quadvec, nonquadvec | expand vec_q
vecg_suffix = gradvec, nongradvec | expand vec_g
+deg_suffix = deg{formcompiler.ufl_variants.degree}
cells = 8 8 8
extension = 1. 1. 1.
@@ -16,7 +17,10 @@ extension = vtu
[formcompiler]
numerical_jacobian = 1, 0 | expand num
exact_solution_expression = g
-compare_l2errorsquared = 1e-4
+compare_l2errorsquared = 2e-4
sumfact = 1
vectorize_quad = 1, 0 | expand vec_q
vectorize_grads = 1, 0 | expand vec_g
+
+[formcompiler.ufl_variants]
+degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_3d.ufl b/test/sumfact/poisson/poisson_3d.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..7cee4c5b15c4c555edeec547de99c2dd95016956
--- /dev/null
+++ b/test/sumfact/poisson/poisson_3d.ufl
@@ -0,0 +1,12 @@
+cell = "hexahedron"
+
+x = SpatialCoordinate(cell)
+c = (0.5 - x[0])**2 + (0.5 - x[1])**2 + (0.5 - x[2])**2
+g = exp(-1.*c)
+f = 2*(3.-2*c)*g
+
+V = FiniteElement("CG", cell, degree, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/sumfact/poisson/poisson_3d_order1.ufl b/test/sumfact/poisson/poisson_3d_order1.ufl
deleted file mode 100644
index 4de9c2392cb46ee2dd31c110e88312acbbbfc978..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_3d_order1.ufl
+++ /dev/null
@@ -1,10 +0,0 @@
-cell = "hexahedron"
-
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", cell=cell)
-
-V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
-u = TrialFunction(V)
-v = TestFunction(V)
-
-forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/sumfact/poisson/poisson_3d_order2.mini b/test/sumfact/poisson/poisson_3d_order2.mini
deleted file mode 100644
index c81eae557caa00dbd1268b54292fee936cd29213..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_3d_order2.mini
+++ /dev/null
@@ -1,22 +0,0 @@
-__name = sumfact_poisson_3d_order2_{__exec_suffix}
-__exec_suffix = {diff_suffix}_{vecq_suffix}_{vecg_suffix}
-
-diff_suffix = numdiff, symdiff | expand num
-vecq_suffix = quadvec, nonquadvec | expand vec_q
-vecg_suffix = gradvec, nongradvec | expand vec_g
-
-cells = 8 8 8
-extension = 1. 1. 1.
-
-[wrapper.vtkcompare]
-name = {__name}
-reference = poisson_ref
-extension = vtu
-
-[formcompiler]
-numerical_jacobian = 1, 0 | expand num
-exact_solution_expression = g
-compare_l2errorsquared = 1e-8
-sumfact = 1
-vectorize_quad = 1, 0 | expand vec_q
-vectorize_grads = 1, 0 | expand vec_g
diff --git a/test/sumfact/poisson/poisson_3d_order2.ufl b/test/sumfact/poisson/poisson_3d_order2.ufl
deleted file mode 100644
index 66870b07d109d4fbe3fe342df81691651232c16e..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_3d_order2.ufl
+++ /dev/null
@@ -1,10 +0,0 @@
-cell = "hexahedron"
-
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", cell=cell)
-
-V = FiniteElement("CG", cell, 2, dirichlet_expression=g)
-u = TrialFunction(V)
-v = TestFunction(V)
-
-forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/sumfact/poisson/poisson_dg.mini b/test/sumfact/poisson/poisson_dg_2d.mini
similarity index 65%
rename from test/sumfact/poisson/poisson_dg.mini
rename to test/sumfact/poisson/poisson_dg_2d.mini
index 747a5b7dac42a836c8a4eb0acdfc59b418f58581..ed6b6795cb4d720560b6761adcee293b70ab98e3 100644
--- a/test/sumfact/poisson/poisson_dg.mini
+++ b/test/sumfact/poisson/poisson_dg_2d.mini
@@ -1,9 +1,10 @@
-__name = sumfact_poisson_dg_{__exec_suffix}
-__exec_suffix = {diff_suffix}_{quadvec_suffix}_{gradvec_suffix}
+__name = sumfact_poisson_dg_2d_{__exec_suffix}
+__exec_suffix = {deg_suffix}_{diff_suffix}_{quadvec_suffix}_{gradvec_suffix}
diff_suffix = numdiff, symdiff | expand num
quadvec_suffix = quadvec, nonquadvec | expand quad
gradvec_suffix = gradvec, nongradvec | expand grad
+deg_suffix = deg{formcompiler.ufl_variants.degree}
cells = 16 16
extension = 1. 1.
@@ -19,3 +20,6 @@ exact_solution_expression = g
compare_l2errorsquared = 1e-4
vectorize_quad = 1, 0 | expand quad
vectorize_grads = 1, 0 | expand grad
+
+[formcompiler.ufl_variants]
+degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg.ufl b/test/sumfact/poisson/poisson_dg_2d.ufl
similarity index 58%
rename from test/sumfact/poisson/poisson_dg.ufl
rename to test/sumfact/poisson/poisson_dg_2d.ufl
index 9019467539150010b6aff3eaebc74216984bbe1e..5ef70869cb556cb6ed3d31ae736a7abd31919edf 100644
--- a/test/sumfact/poisson/poisson_dg.ufl
+++ b/test/sumfact/poisson/poisson_dg_2d.ufl
@@ -1,9 +1,11 @@
cell = "quadrilateral"
-f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", on_intersection=True, cell=cell)
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-1.*c)
+f = 2*(2.-2*c)*g
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
@@ -24,4 +26,4 @@ r = inner(grad(u), grad(v))*dx \
+ theta*g*inner(grad(v), n)*ds \
- gamma*g*v*ds
-forms = [r]
\ No newline at end of file
+forms = [r]
diff --git a/test/sumfact/poisson/poisson_dg_3d.mini b/test/sumfact/poisson/poisson_dg_3d.mini
new file mode 100644
index 0000000000000000000000000000000000000000..7bab7f81c1e5a946738540a9ee0b903e38f22140
--- /dev/null
+++ b/test/sumfact/poisson/poisson_dg_3d.mini
@@ -0,0 +1,25 @@
+__name = sumfact_poisson_dg_3d_{__exec_suffix}
+__exec_suffix = {deg_suffix}_{diff_suffix}_{quadvec_suffix}_{gradvec_suffix}
+
+diff_suffix = numdiff, symdiff | expand num
+quadvec_suffix = quadvec, nonquadvec | expand quad
+gradvec_suffix = gradvec, nongradvec | expand grad
+deg_suffix = deg{formcompiler.ufl_variants.degree}
+
+cells = 8 8 8
+extension = 1. 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand num
+sumfact = 1
+exact_solution_expression = g
+compare_l2errorsquared = 1e-4
+vectorize_quad = 1, 0 | expand quad
+vectorize_grads = 1, 0 | expand grad
+
+[formcompiler.ufl_variants]
+degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_3d.ufl b/test/sumfact/poisson/poisson_dg_3d.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..d9172368b8c8a4c0d8ab81afe3778dbd2af5075f
--- /dev/null
+++ b/test/sumfact/poisson/poisson_dg_3d.ufl
@@ -0,0 +1,29 @@
+cell = "hexahedron"
+
+x = SpatialCoordinate(cell)
+c = (0.5-x[0])**2 + (0.5-x[1])**2 + (0.5-x[2])**2
+g = exp(-1.*c)
+f = 2*(3.-2*c)*g
+
+V = FiniteElement("DG", cell, degree)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+gamma = 1.0
+theta = 1.0
+
+r = inner(grad(u), grad(v))*dx \
+ + inner(n, avg(grad(u)))*jump(v)*dS \
+ + gamma*jump(u)*jump(v)*dS \
+ + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+ - inner(n, grad(u))*v*ds \
+ + gamma*u*v*ds \
+ - theta*u*inner(grad(v), n)*ds \
+ - f*v*dx \
+ + theta*g*inner(grad(v), n)*ds \
+ - gamma*g*v*ds
+
+forms = [r]
diff --git a/test/sumfact/poisson/poisson_dg_only_volume.mini b/test/sumfact/poisson/poisson_dg_only_volume.mini
deleted file mode 100644
index c0153667d89981fe2310cffa23beef6d34e4079a..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_dg_only_volume.mini
+++ /dev/null
@@ -1,20 +0,0 @@
-__name = sumfact_poisson_dg_only_volume_{__exec_suffix}
-
-__exec_suffix = {__num_suffix}_{__sumfact_suffix}
-__num_suffix = numdiff, symdiff |expand num
-__sumfact_suffix = normal, sumfact | expand sumf
-
-cells = 1 1
-extension = 1. 1.
-printresidual = 1
-printmatrix = 1
-
-[wrapper.vtkcompare]
-name = {__name}
-extension = vtu
-
-[formcompiler]
-numerical_jacobian = 1, 0 | expand num
-sumfact = 0, 1 | expand sumf
-print_transformations = 1
-print_transformations_dir = .
diff --git a/test/sumfact/poisson/poisson_dg_only_volume.ufl b/test/sumfact/poisson/poisson_dg_only_volume.ufl
deleted file mode 100644
index 8c8caa9ccdacfb53d980822cbf3e956463fd1e37..0000000000000000000000000000000000000000
--- a/test/sumfact/poisson/poisson_dg_only_volume.ufl
+++ /dev/null
@@ -1,19 +0,0 @@
-cell = "quadrilateral"
-
-f = Expression("return -2.0*x.size();", cell=cell)
-g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
-
-V = FiniteElement("DG", cell, 1)
-
-u = TrialFunction(V)
-v = TestFunction(V)
-
-n = FacetNormal(cell)('+')
-
-gamma = 1.0
-theta = 1.0
-
-r = inner(grad(u), grad(v))*dx \
- - f*v*dx
-
-forms = [r]
\ No newline at end of file
diff --git a/test/sumfact/poisson/poisson_fastdg.mini b/test/sumfact/poisson/poisson_fastdg.mini
index 87e642e977d462e8f5b1577c12e7512187b1e351..e0cf2f268ddb66f043503e0ca00c1112b369884f 100644
--- a/test/sumfact/poisson/poisson_fastdg.mini
+++ b/test/sumfact/poisson/poisson_fastdg.mini
@@ -19,3 +19,6 @@ compare_l2errorsquared = 1e-4
vectorize_quad = 1, 0 | expand quadvec
vectorize_grads = 1, 0 | expand gradvec
fastdg = 1
+
+[formcompiler.ufl_variants]
+degree = 1