Merge branch 'feature/cellwise-constant' into 'master'

Feature/cellwise constant

Make it possible to mark Expressions as cellwise const and put evaluation in pdelab outside the quadrature loop. Closes issue #27.

The implementation is kind of sketchy but I wanted to go to weekend. Some notes/TODO:

- I had to add a statement to the ufl file. For cellwise const expressions a dg0 dummy space is used (-> cell wise const), else we use a dg1 space (in execution.py).
- I now use 'default_declaration' from loopy.py for defining parameter objects.

See merge request !14

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)

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)

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...

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)

test/poisson/poisson_cellwise_constant.mini

+__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

test/poisson/poisson_cellwise_constant.ufl

+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]