Major update

Some work on parameter classes, a bugfix on the splitting algorithm
some simplifications in the generation stage

python/dune/perftool/cgen/clazz.py

@@ -39,10 +39,12 @@ class ClassMember(Generable):
         if isinstance(member, str):
             from cgen import Line
             self.member = Line(member)
+        else:
+            from collections import Iterable
+            assert all(isinstance(m, str) for m in self.member)
 
-        # We only consider a Generable or a list thereof as member
-        from collections import Iterable
-        assert isinstance(self.member, Generable) or (isinstance(self.member, Iterable) and all(isinstance(m, Generable) for m in self.member))
+            from cgen import LiteralLines
+            self.member = LiteralLines("\n" + "\n".join(self.member))
 
     def generate(self):
         yield "\n\n"

python/dune/perftool/compile.py

@@ -48,7 +48,7 @@ def read_ufl(uflfile):
     form = transform_form(form, reindexing)
 #    form = transform_form(form, split_arguments)
 
-    return form
+    return form, data.object_names
 
 
 def generate_driver(form, filename):
@@ -78,7 +78,7 @@ def generate_driver(form, filename):
 # This function is the entrypoint of the ufl2pdelab executable
 def compile_form():
     from dune.perftool.options import get_option
-    form = read_ufl(get_option("uflfile"))
+    form, namedata = read_ufl(get_option("uflfile"))
 
     from dune.perftool.generation import cache_context
     if get_option("driver_file"):
@@ -87,7 +87,7 @@ def compile_form():
 
     if get_option("operator_file"):
         from dune.perftool.pdelab.localoperator import generate_localoperator_kernels
-        kernels = generate_localoperator_kernels(form)
+        kernels = generate_localoperator_kernels(form, namedata)
 
         # TODO insert sophisticated analysis/feedback loops here
         if get_option("interactive"):

python/dune/perftool/generation/__init__.py

@@ -29,4 +29,5 @@ from dune.perftool.generation.loopy import (domain,
 from dune.perftool.generation.context import (cache_context,
                                               generic_context,
                                               get_generic_context_value,
+                                              namedata_context,
                                               )

python/dune/perftool/generation/context.py

@@ -53,3 +53,6 @@ def generic_context(key, value):
 
 def get_generic_context_value(key):
     return _generic_context_cache[key]
+
+import functools
+namedata_context=functools.partial(generic_context, "namedata")

python/dune/perftool/loopy/transformer.py

@@ -38,8 +38,22 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper):
     def coefficient(self, o):
         # All trial functions should already be handled
         assert o.count() != 0
-        # TODO implement non-trialfunction coefficients
-        raise NotImplementedError
+
+        # We expect all coefficients to be of type Expression!
+        from dune.perftool.ufl.execution import Expression
+        assert isinstance(o, Expression)
+
+        # Determine the name of the parameter function
+        from dune.perftool.generation import get_generic_context_value
+        name = get_generic_context_value("namedata")[id(o)]
+
+        # Trigger the generation of code for this thing in the parameter class
+        from dune.perftool.pdelab.parameter import parameter_function
+        parameter_function(name, o)
+
+        # And return a symbol
+        from pymbolic.primitives import Variable
+        return Variable(name)
 
     def index_sum(self, o):
         from dune.perftool.ufl.shape import determine_shape

python/dune/perftool/pdelab/basis.py

@@ -181,7 +181,7 @@ def evaluate_basis(element, name):
     temporary_variable(name, shape=(name_lfs_bound(element),), shape_impl=('vec',))
     cache = name_localbasis_cache(element)
     lfs = name_lfs(element)
-    qp = name_quadrature_point()
+    qp = name_quadrature_position()
     instruction(inames=(quadrature_iname(),
                         ),
                 code='auto& {} = {}.evaluateFunction({}, {}.finiteElement().localBasis());'.format(name,

python/dune/perftool/pdelab/geometry.py

@@ -7,9 +7,23 @@ from dune.perftool.pdelab.quadrature import (name_quadrature_position,
                                              )
 
 
+@symbol
+def name_elementgeometry():
+    return 'eg'
+
+
+@symbol
+def name_element():
+    eg = name_elementgeometry()
+    return "{}.element()".format(eg)
+
+
 @preamble
 def define_geometry(name):
-    return "auto {} = eg.geometry();".format(name)
+    eg = name_elementgeometry()
+    return "auto {} = {}.geometry();".format(name,
+                                             eg
+                                             )
 
 
 @symbol
@@ -57,7 +71,8 @@ def define_jacobian_inverse_transposed(name):
     return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
                                                                                geo,
                                                                                pos,
-                                                                               )
+                                                                               ),
+                               assignees=frozenset({name}),
                                )
 
 

python/dune/perftool/pdelab/localoperator.py

@@ -178,7 +178,10 @@ class AssemblyMethod(ClassMember):
     def __init__(self, signature, kernel):
         from loopy import generate_code
         from cgen import LiteralLines, Block
-        content = [LiteralLines('\n' + '\n'.join(signature)), Block([LiteralLines('\n' + generate_code(kernel)[0])])]
+        content = signature
+        content.append('{')
+        content.extend('  ' + l for l in generate_code(kernel)[0].split('\n'))
+        content.append('}')
         ClassMember.__init__(self, content)
 
 
@@ -201,7 +204,7 @@ def cgen_class_from_cache(tag, members=[]):
     return Class(basename, base_classes=base_classes, members=members + pm, constructors=[constructor], tparam_decls=tparams)
 
 
-def generate_localoperator_kernels(form):
+def generate_localoperator_kernels(form, namedata):
     # For the moment, I do assume that there is but one integral of each type. This might differ
     # if you use different quadrature orders for different terms.
     assert len(form.integrals()) == len(set(i.integral_type() for i in form.integrals()))
@@ -220,15 +223,20 @@ def generate_localoperator_kernels(form):
     localoperator_basename()
     lop_template_ansatz_gfs()
     lop_template_test_gfs()
+    from dune.perftool.pdelab.parameter import parameterclass_basename
+    parameterclass_basename()
 
     base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")
 
     # Have a data structure collect the generated kernels
     operator_kernels = {}
 
+    from dune.perftool.generation import namedata_context
+
     # Generate the necessary residual methods
     for integral in form.integrals():
-        kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
+        with namedata_context(namedata):
+            kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
         operator_kernels[(integral.integral_type(), 'residual')] = kernel
 
     # Generate the necessary jacobian methods
@@ -241,7 +249,8 @@ def generate_localoperator_kernels(form):
         jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
 
         for integral in jacform.integrals():
-            kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
+            with namedata_context(namedata):
+                kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
             operator_kernels[(integral.integral_type(), 'jacobian')] = kernel
 
     # TODO: JacobianApply for matrix-free computations.

python/dune/perftool/pymbolic/inameset.py

@@ -9,6 +9,12 @@ class INameMapper(CombineMapper):
         else:
             return frozenset([str(i) for i in e.index])
 
+    def map_constant(self, e):
+        return frozenset({})
+
+    def map_algebraic_leaf(self, e):
+        return frozenset({})
+
     def combine(self, values):
         return frozenset().union(*values)
 

python/dune/perftool/ufl/execution.py

@@ -32,3 +32,13 @@ def Coefficients(element):
 
 def TrialFunctions(element):
     return split(TrialFunction(element))
+
+
+class Expression(ufl.Coefficient):
+    def __init__(self, expr, is_global=True):
+        assert isinstance(expr, str)
+        self.c_expr = expr
+        self.is_global = is_global
+
+        # Initialize a coefficient with a dummy finite element map.
+        ufl.Coefficient.__init__(self, FiniteElement("DG", "triangle", 0))

python/dune/perftool/ufl/transformations/extract_accumulation_terms.py

@@ -16,14 +16,13 @@ import itertools
 
 
 class _ReplacementDict(dict):
-    def __init__(self, *args):
+    def __init__(self, good=[], bad=[]):
         dict.__init__(self)
-        for a in args:
+        for a in bad:
+            self[a] = Zero()
+        for a in good:
             self[a] = a
 
-    def __getitem__(self, key):
-        return self.get(key, Zero())
-
 
 @ufl_transformation(name="accterms2", extraction_lambda=lambda l: l)
 def split_into_accumulation_terms(expr):
@@ -35,13 +34,23 @@ def split_into_accumulation_terms(expr):
     if len(filter(lambda ma: ma.argexpr.count() == 1, mod_args)) == 0:
         for arg in mod_args:
             # Do the replacement on the expression
-            accumulation_terms.append(replace_expression(expr, replacemap=_ReplacementDict(arg.expr)))
+            accumulation_terms.append(replace_expression(expr,
+                                                         replacemap=_ReplacementDict(good=(arg.expr,),
+                                                                                     bad=[ma.expr for ma in mod_args],
+                                                                                     )
+                                                         )
+                                      )
     # and now the case of a rank 2 form:
     else:
         for arg1, arg2 in itertools.product(filter(lambda ma: ma.argexpr.count() == 0, mod_args),
                                             filter(lambda ma: ma.argexpr.count() == 1, mod_args)
                                             ):
-            accumulation_terms.append(replace_expression(expr, replacemap=_ReplacementDict(arg1.expr, arg2.expr)))
+            accumulation_terms.append(replace_expression(expr,
+                                                         replacemap=_ReplacementDict(good=(arg1.expr, arg2.expr),
+                                                                                     bad=[ma.expr for ma in mod_args],
+                                                                                     )
+                                                         )
+                                      )
 
     # and return the result
     return accumulation_terms

test/CMakeLists.txt

@@ -5,3 +5,12 @@ add_generated_executable(UFLFILE ../examples/laplace.ufl
 
 dune_add_system_test(TARGET laplace
                      INIFILE laplace.mini)
+
+
+add_generated_executable(UFLFILE ../examples/poisson.ufl
+                         TARGET poisson
+                         FORM_COMPILER_ARGS --numerical-jacobian
+                         )
+
+dune_add_system_test(TARGET poisson
+                     INIFILE laplace.mini)