diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 630ed6e871963f0f2483282865586550ed955946..78e872f6b3410c0323f0b91edeca496c7c7132cc 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -152,8 +152,10 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
name = 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(o)
from dune.perftool.pdelab.parameter import intersection_parameter_function
- intersection_parameter_function(name, subdomain_data, t='int')
+ intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int')
predicates = predicates.union(['{} == {}'.format(name, self.subdomain_id)])
@@ -242,14 +244,17 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
from dune.perftool.generation import get_global_context_value
name = get_global_context_value("data").object_names[id(o)]
+ from ufl.checks import is_cellwise_constant
+ cellwise_constant = is_cellwise_constant(o)
+
# Trigger the generation of code for this thing in the parameter class
from dune.perftool.pdelab.parameter import (cell_parameter_function,
intersection_parameter_function,
)
if o.on_intersection:
- intersection_parameter_function(name, o)
+ intersection_parameter_function(name, o, cellwise_constant)
else:
- cell_parameter_function(name, o, self.restriction)
+ cell_parameter_function(name, o, self.restriction, cellwise_constant)
# And return a symbol
return Variable(name)
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index 05c21e823571df0e5973dc34357df2f17fc476d8..79fede9c1eb72e20613b2db92cd6df8b1463050c 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -89,6 +89,47 @@ def define_parameter_function_class_member(name, expr, t, cell):
return result
+@preamble
+def evaluate_cellwise_constant_parameter_function(name, restriction):
+ param = name_paramclass()
+ entity = name_cell(restriction)
+ pos = name_quadrature_position_in_cell(restriction)
+
+ from dune.perftool.generation.loopy import valuearg
+ import numpy
+ valuearg(name, dtype=numpy.float64)
+
+ return '{} = {}.{}({}, {});'.format(name,
+ name_paramclass(),
+ name,
+ entity,
+ pos,
+ )
+
+
+@preamble
+def evaluate_intersectionwise_constant_parameter_function(name):
+ # Check that this is not a volume term, as that would not be well-defined
+ from dune.perftool.generation import get_global_context_value
+ it = get_global_context_value("integral_type")
+ assert it is not 'cell'
+
+ param = name_paramclass()
+ intersection = name_intersection()
+ pos = name_quadrature_position()
+
+ from dune.perftool.generation.loopy import valuearg
+ import numpy
+ valuearg(name, dtype=numpy.float64)
+
+ return '{} = {}.{}({}, {});'.format(name,
+ name_paramclass(),
+ name,
+ intersection,
+ pos,
+ )
+
+
def evaluate_cell_parameter_function(name, restriction):
param = name_paramclass()
entity = name_cell(restriction)
@@ -131,20 +172,30 @@ def construct_nested_fieldvector(t, shape):
@cached
-def cell_parameter_function(name, expr, restriction, t='double'):
+def cell_parameter_function(name, expr, restriction, cellwise_constant, t='double'):
shape = expr.ufl_element().value_shape()
shape_impl = ('fv',) * len(shape)
t = construct_nested_fieldvector(t, shape)
- temporary_variable(name, shape=shape, shape_impl=shape_impl)
define_parameter_function_class_member(name, expr, t, True)
- evaluate_cell_parameter_function(name, restriction)
+ if cellwise_constant:
+ from dune.perftool.generation.loopy import default_declaration
+ default_declaration(name, shape, shape_impl)
+ evaluate_cellwise_constant_parameter_function(name, restriction)
+ else:
+ temporary_variable(name, shape=shape, shape_impl=shape_impl)
+ evaluate_cell_parameter_function(name, restriction)
@cached
-def intersection_parameter_function(name, expr, t='double'):
+def intersection_parameter_function(name, expr, cellwise_constant, t='double'):
shape = expr.ufl_element().value_shape()
shape_impl = ('fv',) * len(shape)
t = construct_nested_fieldvector(t, shape)
- temporary_variable(name, shape=shape, shape_impl=shape_impl)
define_parameter_function_class_member(name, expr, t, False)
- evaluate_intersection_parameter_function(name)
+ if cellwise_constant:
+ from dune.perftool.generation.loopy import default_declaration
+ default_declaration(name, shape, shape_impl)
+ evaluate_intersectionwise_constant_parameter_function(name)
+ else:
+ temporary_variable(name, shape=shape, shape_impl=shape_impl)
+ evaluate_intersection_parameter_function(name)
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 5fe6975b233a73a2629186489742f5e897b6d7f2..fb0aa6b5d008dd028f3d7903c3b57a43ed728842 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -54,17 +54,20 @@ def TrialFunctions(element):
class Expression(Coefficient):
- def __init__(self, expr, is_global=True, on_intersection=False, cell_type="triangle"):
+ def __init__(self, expr, is_global=True, on_intersection=False, cell_type="triangle", cellwise_constant=False):
assert isinstance(expr, str)
self.c_expr = expr
self.is_global = is_global
self.on_intersection = on_intersection
# Avoid ufl complaining about not matching dimension/cells
- _dg0 = FiniteElement("DG", cell_type, 0)
+ if cellwise_constant:
+ _dg = FiniteElement("DG", cell_type, 0)
+ else:
+ _dg = FiniteElement("DG", cell_type, 1)
# Initialize a coefficient with a dummy finite element map.
- Coefficient.__init__(self, _dg0)
+ Coefficient.__init__(self, _dg)
# TODO the subdomain_data code needs a uflid, not idea how to get it here
# The standard way through class decorator fails here...
@@ -73,17 +76,20 @@ class Expression(Coefficient):
class VectorExpression(Coefficient):
- def __init__(self, expr, is_global=True, on_intersection=False, cell_type="triangle"):
+ def __init__(self, expr, is_global=True, on_intersection=False, cell_type="triangle", cellwise_constant=False):
assert isinstance(expr, str)
self.c_expr = expr
self.is_global = is_global
self.on_intersection = on_intersection
# Avoid ufl complaining about not matching dimension/cells
- _dg0vec = VectorElement("DG", cell_type, 0)
+ if cellwise_constant:
+ _dgvec = VectorElement("DG", cell_type, 0)
+ else:
+ _dgvec = VectorElement("DG", cell_type, 1)
# Initialize a coefficient with a dummy finite element map.
- Coefficient.__init__(self, _dg0vec)
+ Coefficient.__init__(self, _dgvec)
# TODO the subdomain_data code needs a uflid, not idea how to get it here
# The standard way through class decorator fails here...
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index 1a7aa445b5fc8cc706650cea3830dbce78b43559..2bc18c5c1e4c0baf87c33f065d79b971b762560b 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -132,6 +132,16 @@ dune_add_system_test(TARGET poisson_symdiff_matrix_free
# SCRIPT dune_vtkcompare.py
# )
+
+# 13. Poisson Test Case: test constant expressions
+add_generated_executable(UFLFILE poisson_cellwise_constant.ufl
+ TARGET poisson_cellwise_constant
+ FORM_COMPILER_ARGS --exact-solution-expression g --compare-l2errorsquared 1e-6
+ )
+dune_add_system_test(TARGET poisson_cellwise_constant
+ INIFILE poisson_cellwise_constant.mini
+ )
+
# Add the reference code with the PDELab localoperator that produced
# the reference vtk file
add_executable(poisson_dg_ref reference_main.cc)
diff --git a/test/poisson/poisson_cellwise_constant.mini b/test/poisson/poisson_cellwise_constant.mini
new file mode 100644
index 0000000000000000000000000000000000000000..5b1ceb7dbd66eab1366e0cf89dda846e9d269fac
--- /dev/null
+++ b/test/poisson/poisson_cellwise_constant.mini
@@ -0,0 +1,11 @@
+__name = poisson_cellwise_constant
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_ref
+extension = vtu
diff --git a/test/poisson/poisson_cellwise_constant.ufl b/test/poisson/poisson_cellwise_constant.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..b1f3c08ebb2aab9a0b040e24bbb589c2ea767bcc
--- /dev/null
+++ b/test/poisson/poisson_cellwise_constant.ufl
@@ -0,0 +1,9 @@
+a = Expression("return 1.0;", cellwise_constant=True)
+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());")
+
+V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+forms = [(a*inner(grad(u), grad(v)) - f*v)*dx]