Merge branch 'feature/code-restructuring' into 'master'

Feature/code restructuring

All the restructuring stuff we talked about.

Tested.

See merge request !38

.gitmodules

-[submodule "python/cgen"]
-	path = python/cgen
-	url = https://github.com/inducer/cgen.git
 [submodule "python/loopy"]
 	path = python/loopy
 	url = https://github.com/dokempf/loopy.git

examples/dg.ufl

-V = FiniteElement('CG',triangle,2)
-u = TrialFunction(V)
-v = TestFunction(V)
-x,y = SpatialCoordinate(triangle)
-n = FacetNormal(triangle)
-
-d_f = FacetArea(triangle)
-
-alpha = 1.3
-
-rho = as_matrix((
-    (1.0,0.2),
-    (0.2,1.0)
-    ))
-
-#r = (avg(u)*inner(n('+'),jump(grad(v))) + inner(n('+'),jump(rho*grad(u)))*avg(v) + alpha / d_f * jump(u)*jump(v))*dS(0)
-r = (avg(u)*jump(v) + avg(v)*jump(u))*dS
-
-forms = [ r ]

examples/mixed.ufl

-V =  VectorElement("CG", triangle, 1, dim=3)
-u = Coefficient(V)
-v = TestFunction(V)
-
-r = inner(u,v)*dx
-
-forms = [r]

examples/stokes2d.ufl

-cell = triangle
-P2 = VectorElement("Lagrange", cell, 2)
-P1 = FiniteElement("Lagrange", cell, 1)
-TH = P2 * P1
-
-(v, q) = TestFunctions(TH)
-(u, p) = TrialFunctions(TH)
-
-f = Coefficient(P2)
-
-r = (inner(grad(v), grad(u)) - div(v)*p + q*div(u) - dot(v, f))*dx
-
-forms = [r]

python/CMakeLists.txt

@@ -10,7 +10,6 @@ create_virtualenv_wrapper(COMMANDS python -m IPython notebook
 
 # Install all the external packages that we have as submodules
 dune_install_python_package(PATH pymbolic)
-dune_install_python_package(PATH cgen)
 dune_install_python_package(PATH loopy)
 dune_install_python_package(PATH ufl)
 

cgen @ 716ea0fd

-Subproject commit 716ea0fda102bc5b829122dd578566ce42f37912

python/dune/perftool/loopy/functions.py

-from dune.perftool.generation import function_mangler
-from loopy import CallMangleInfo
-from loopy.symbolic import FunctionIdentifier
-from loopy.types import NumpyType
-
-import numpy
-
-
-class LFSChild(FunctionIdentifier):
-    def __init__(self, lfs):
-        self.lfs = lfs
-
-    def __getinitargs__(self):
-        return (self.lfs,)
-
-    @property
-    def name(self):
-        return '{}.child'.format(self.lfs)
-
-
-@function_mangler
-def lfs_child_mangler(target, func, dtypes):
-    if isinstance(func, LFSChild):
-        return CallMangleInfo(func.name, (NumpyType(str),), (NumpyType(numpy.int32),))
-
-
-class CoefficientAccess(FunctionIdentifier):
-    def __init__(self, container):
-        self.container = container
-
-    def __getinitargs__(self):
-        return (self.container,)
-
-    @property
-    def name(self):
-        return self.container
-
-
-@function_mangler
-def coefficient_mangler(target, func, dtypes):
-    if isinstance(func, CoefficientAccess):
-        return CallMangleInfo(func.name, (NumpyType(numpy.float64),), (NumpyType(str), NumpyType(numpy.int32)))
-
-
-class PDELabAccumulationFunction(FunctionIdentifier):
-    def __init__(self, accumobj, rank):
-        self.accumobj = accumobj
-        self.rank = rank
-
-        assert rank in (1, 2)
-
-    def __getinitargs__(self):
-        return (self.accumobj, self.rank)
-
-    @property
-    def name(self):
-        return '{}.accumulate'.format(self.accumobj)
-
-
-@function_mangler
-def accumulation_mangler(target, func, dtypes):
-    if isinstance(func, PDELabAccumulationFunction):
-        if func.rank == 1:
-            return CallMangleInfo(func.name,
-                                  (),
-                                  (NumpyType(str),
-                                   NumpyType(numpy.int32),
-                                   NumpyType(numpy.float64),
-                                   )
-                                  )
-        if func.rank == 2:
-            return CallMangleInfo(func.name,
-                                  (),
-                                  (NumpyType(str),
-                                   NumpyType(numpy.int32),
-                                   NumpyType(str),
-                                   NumpyType(numpy.int32),
-                                   NumpyType(numpy.float64),
-                                   )
-                                  )

python/dune/perftool/pdelab/__init__.py

@@ -20,15 +20,6 @@ def name_index(index):
     raise NotImplementedError
 
 
-def pymbolic_index(index):
-    from ufl.classes import FixedIndex
-    if isinstance(index, FixedIndex):
-        return index._value
-    else:
-        from pymbolic.primitives import Variable
-        return Variable(name_index(index))
-
-
 def restricted_name(name, restriction):
     from dune.perftool import Restriction
     if restriction == Restriction.NONE:

python/dune/perftool/pdelab/argument.py

-""" Generator functions related to trial and test functions and the accumulation loop"""
+""" Generator functions related to any input and output of accumulation kernels
 
-from dune.perftool.options import get_option
+Namely:
+* Coefficient containers
+* accumulation object (r, jac...)
+"""
 
+from dune.perftool.options import get_option
 from dune.perftool.generation import (domain,
+                                      function_mangler,
                                       iname,
                                       pymbolic_expr,
                                       globalarg,
@@ -15,14 +20,76 @@ from dune.perftool.pdelab import (name_index,
                                   )
 from dune.perftool.pdelab.basis import (evaluate_coefficient,
                                         evaluate_coefficient_gradient,
-                                        lfs_iname,
                                         name_basis,
-                                        name_lfs_bound,
                                         )
+from dune.perftool.pdelab.spaces import (lfs_iname,
+                                         name_lfs_bound,
+                                         )
 from dune.perftool import Restriction
+
 from pymbolic.primitives import Call, Subscript, Variable
 
-import loopy
+from loopy import CallMangleInfo
+from loopy.symbolic import FunctionIdentifier
+from loopy.types import NumpyType
+
+import numpy
+
+
+class CoefficientAccess(FunctionIdentifier):
+    def __init__(self, container):
+        self.container = container
+
+    def __getinitargs__(self):
+        return (self.container,)
+
+    @property
+    def name(self):
+        return self.container
+
+
+@function_mangler
+def coefficient_mangler(target, func, dtypes):
+    if isinstance(func, CoefficientAccess):
+        return CallMangleInfo(func.name, (NumpyType(numpy.float64),), (NumpyType(str), NumpyType(numpy.int32)))
+
+
+class PDELabAccumulationFunction(FunctionIdentifier):
+    def __init__(self, accumobj, rank):
+        self.accumobj = accumobj
+        self.rank = rank
+
+        assert rank in (1, 2)
+
+    def __getinitargs__(self):
+        return (self.accumobj, self.rank)
+
+    @property
+    def name(self):
+        return '{}.accumulate'.format(self.accumobj)
+
+
+@function_mangler
+def accumulation_mangler(target, func, dtypes):
+    if isinstance(func, PDELabAccumulationFunction):
+        if func.rank == 1:
+            return CallMangleInfo(func.name,
+                                  (),
+                                  (NumpyType(str),
+                                   NumpyType(numpy.int32),
+                                   NumpyType(numpy.float64),
+                                   )
+                                  )
+        if func.rank == 2:
+            return CallMangleInfo(func.name,
+                                  (),
+                                  (NumpyType(str),
+                                   NumpyType(numpy.int32),
+                                   NumpyType(str),
+                                   NumpyType(numpy.int32),
+                                   NumpyType(numpy.float64),
+                                   )
+                                  )
 
 
 def name_trialfunction_gradient(element, restriction, component):
@@ -67,22 +134,6 @@ def name_apply_function(element, restriction, component):
     return name
 
 
-def name_testfunctionspace(restriction):
-    return restricted_name("lfsv", restriction)
-
-
-def name_trialfunctionspace(restriction):
-    return restricted_name("lfsu", restriction)
-
-
-def type_testfunctionspace():
-    return "LFSV"
-
-
-def type_trialfunctionspace():
-    return "LFSU"
-
-
 def name_coefficientcontainer(restriction):
     name = restricted_name("x", restriction)
     return name
@@ -97,14 +148,13 @@ def name_applycontainer(restriction):
 def pymbolic_coefficient(container, lfs, index):
     # TODO introduce a proper type for local function spaces!
     if isinstance(lfs, str):
-        valuearg(lfs, dtype=loopy.types.NumpyType("str"))
+        valuearg(lfs, dtype=NumpyType("str"))
 
     # If the LFS is not yet a pymbolic expression, make it one
     from pymbolic.primitives import Expression
     if not isinstance(lfs, Expression):
         lfs = Variable(lfs)
 
-    from dune.perftool.loopy.functions import CoefficientAccess
     return Call(CoefficientAccess(container), (lfs, Variable(index),))
 
 
@@ -112,21 +162,6 @@ def type_coefficientcontainer():
     return "X"
 
 
-def name_argumentspace(ma):
-    from ufl.classes import Argument, Coefficient
-    if isinstance(ma.argexpr, Argument):
-        if ma.argexpr.number() == 0:
-            return name_testfunctionspace(ma.restriction)
-        if ma.argexpr.number() == 1:
-            return name_trialfunctionspace(ma.restriction)
-    if isinstance(ma.argexpr, Coefficient):
-        # Index 0 is reserved for trialfunction, index 1 is reserved for jacobian apply function
-        assert ma.argexpr.count() < 2
-        return name_trialfunctionspace(ma.restriction)
-    # We should never encounter an argument other than 0 or 1
-    assert False
-
-
 def name_jacobian(restriction1, restriction2):
     # Restrictions may only differ if NONE
     if (restriction1 == Restriction.NONE) or (restriction2 == Restriction.NONE):
@@ -147,7 +182,6 @@ def type_residual():
 
 
 def name_accumulation_variable(restrictions=(Restriction.NONE,)):
-    from dune.perftool.generation import get_global_context_value
     ft = get_global_context_value("form_type")
     if ft == 'residual' or ft == 'jacobian_apply':
         return name_residual(*restrictions)
@@ -157,7 +191,6 @@ def name_accumulation_variable(restrictions=(Restriction.NONE,)):
 
 
 def type_accumulation_variable():
-    from dune.perftool.generation import get_global_context_value
     ft = get_global_context_value("form_type")
     if ft == 'residual' or ft == 'jacobian_apply':
         return type_residual()

python/dune/perftool/pdelab/basis.py

@@ -2,14 +2,19 @@
 
 from dune.perftool.generation import (cached,
                                       class_member,
-                                      domain,
                                       generator_factory,
-                                      iname,
                                       include_file,
                                       instruction,
                                       preamble,
                                       temporary_variable,
                                       )
+from dune.perftool.pdelab.spaces import (lfs_child,
+                                         lfs_iname,
+                                         name_leaf_lfs,
+                                         name_lfs,
+                                         name_lfs_bound,
+                                         type_gfs,
+                                         )
 from dune.perftool.pdelab.quadrature import (name_quadrature_position_in_cell,
                                              quadrature_iname,
                                              )
@@ -22,6 +27,7 @@ from dune.perftool.pdelab.localoperator import (lop_template_ansatz_gfs,
                                                 lop_template_test_gfs,
                                                 )
 from dune.perftool.pdelab.driver import FEM_name_mangling
+
 from dune.perftool.pdelab import restricted_name
 from pymbolic.primitives import Product, Subscript, Variable
 from loopy import Reduction
@@ -44,185 +50,6 @@ def name_localbasis_cache(element):
     return name
 
 
-@preamble
-def define_lfs_bound(lfs, bound):
-    return 'auto {} = {}.size();'.format(bound, lfs)
-
-
-def name_lfs_bound(lfs):
-    # LFS might either be given by an UFL element or by a string describig its name
-    from ufl import FiniteElementBase
-    if isinstance(lfs, FiniteElementBase):
-        return name_lfs_bound(name_lfs(lfs))
-
-    bound = '{}_size'.format(lfs)
-    define_lfs_bound(lfs, bound)
-    return bound
-
-
-@preamble
-def using_indices():
-    return "using namespace Dune::Indices;"
-
-
-@preamble
-def define_lfs(name, father, child):
-    using_indices()
-    return "auto {} = child({}, _{});".format(name, father, child)
-
-
-def lfs_child(lfs, children, shape=None, symmetry=False):
-    from pymbolic.primitives import Call, Product, Sum
-    # Old pre-TensorElement implementation kept for comaptibility
-    if shape is None:
-        indices = (Variable(children[0]),)
-    else:
-        if symmetry and len(children) == 2:
-            # I do not want to think about tensors of rank > 2 right now
-            i, j = children
-            if i > j:
-                j, i = i, j
-            i = Variable(i)
-            j = Variable(j)
-            n = len(children)
-
-            indices = (Sum((Product((n - 1, i)), Product((.5, i, 1 - i)), j)),)
-        else:
-            children = tuple(Variable(c) for c in children)
-            indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),)
-
-    from dune.perftool.loopy.functions import LFSChild
-    return Call(LFSChild(lfs), indices)
-
-
-@generator_factory(cache_key_generator=lambda e, r, **kw: (e, r))
-def name_leaf_lfs(leaf_element, restriction, val=None):
-    """ This function just caches leaf lfs names based on the
-    element. The resulting local function spaces are useful only
-    for size information. OTOH, they are available with just the
-    leaf element available (as seen in basis evaluation).
-    """
-    # This generator function should be prepoluted by the lfs tree traversal,
-    # so val should always be None when we actually want the result
-    assert val
-
-    return val
-
-
-@generator_factory(cache_key_generator=lambda e, r, c, **kw: (e, r, c))
-def name_lfs(element, restriction, component, prefix=None):
-    # Omitting the prefix is only valid upon a second call, which will
-    # result in a cache hit.
-    assert prefix
-
-    def _name_lfs(prefix, component):
-        name = prefix
-        if len(component) > 0:
-            name = name + '_' + '_'.join(str(i) for i in component)
-        return name
-
-    name = _name_lfs(prefix, component)
-    if len(component) > 0:
-        father = _name_lfs(prefix, component[:-1])
-        # If this localfunction space is the child of another one, trigger
-        # the extraction preamble. Necessary before going into recursion
-        # for having the correct (top-down) order of preambles
-        define_lfs(name, father, component[-1])
-
-    # Recurse into the given element to define all other local function spaces!
-    from ufl import MixedElement
-    from ufl.functionview import select_subelement
-    from ufl.classes import FixedIndex
-    subel = select_subelement(element, component)
-    if isinstance(subel, MixedElement):
-        for i in range(subel.num_sub_elements()):
-            name_lfs(element, restriction, component + (FixedIndex(i),), prefix=prefix)
-
-    # Cache the name for the subelement
-    name_leaf_lfs(subel, restriction, val=name)
-
-    # Now return the prefix!
-    return name
-
-
-@generator_factory(cache_key_generator=lambda e, **kw: e)
-def type_gfs(element, basetype=None, index_stack=None):
-    # Omitting basetype and index_stack is only valid upon a second call,
-    # which will result in a cache hit.
-    assert basetype
-    assert index_stack is not None
-
-    # Additionally, element is expected to be a ufl finite element
-    from ufl import FiniteElementBase
-    assert isinstance(element, FiniteElementBase)
-
-    # Recurse into the given element to define all other local function spaces!
-    from ufl import MixedElement
-    from ufl.classes import FixedIndex
-    if isinstance(element, MixedElement):
-        for i, subelem in enumerate(element.sub_elements()):
-            type_gfs(subelem, basetype=basetype, index_stack=index_stack + (FixedIndex(i),))
-
-    if len(index_stack) == 0:
-        return basetype
-    else:
-        include_file("dune/typetree/childextraction.hh", filetag="operatorfile")
-        return 'Dune::TypeTree::Child<{},{}>'.format(basetype, ','.join(str(i) for i in index_stack))
-
-
-def traverse_lfs_tree(arg):
-    from dune.perftool.ufl.modified_terminals import ModifiedArgumentDescriptor
-    assert isinstance(arg, ModifiedArgumentDescriptor)
-
-    # First we need to determine the basename as given in the signature of
-    # this kernel method!
-    from dune.perftool.pdelab.argument import name_argumentspace
-    lfs_basename = name_argumentspace(arg)
-    from dune.perftool.pdelab.localoperator import lop_template_gfs
-    gfs_basename = lop_template_gfs(arg)
-
-    # Now start recursively extracting local function spaces and fill the cache with
-    # all those values. That way we can later get a correct local function space with
-    # just the ufl finite element.
-    from ufl.classes import MultiIndex
-    name_lfs(arg.argexpr.ufl_element(), arg.restriction, MultiIndex(()), prefix=lfs_basename)
-    type_gfs(arg.argexpr.ufl_element(), basetype=gfs_basename, index_stack=())
-
-
-@generator_factory(item_tags=("iname",), cache_key_generator=lambda e, r, c: (e, c))
-def _lfs_iname(element, restriction, context):
-    lfs = name_leaf_lfs(element, restriction)
-    bound = name_lfs_bound(lfs)
-
-    name = "{}_{}_index".format(lfs, context)
-    domain(name, bound)
-
-    return name
-
-
-def lfs_iname(element, restriction, count=None, context=''):
-    """ Get the iname to iterate over the local function space given by element
-
-    Arguments:
-    ----------
-    element: ufl.FiniteElementBase
-        The finite element this local function space belongs to
-    argcount: int
-        Use to realize double nesting in case of jacobians
-    context: str
-        Some generation methods will require you to duplicate an iname for
-        a given purpose, see the 'Loops and dependencies' of the loopy docs:
-        https://documen.tician.de/loopy/tutorial.html#loops-and-dependencies
-    """
-    assert not ((context == '') and (count is None))
-    if count is not None:
-        if context != '':
-            context = "{}_{}".format(count, context)
-        else:
-            context = str(count)
-    return _lfs_iname(element, restriction, context)
-
-
 @preamble
 def declare_cache_temporary(element, restriction, name, which):
     t_cache = type_localbasis_cache(element)

python/dune/perftool/pdelab/geometry.py

@@ -15,56 +15,6 @@ from ufl.algorithms import MultiFunction
 from pymbolic.primitives import Variable
 
 
-class GeometryMapper(MultiFunction):
-    """
-    A collection of visitors for geometry related UFL nodes
-    NB: This is kind of 'abstract' as it needs to be combined
-        with a ModifiedTerminalTracker through multi inheritance.
-    """
-    def __init__(self):
-        super(GeometryMapper, self).__init__()
-
-    def facet_normal(self, o):
-        # The normal must be restricted to be well-defined
-        assert self.restriction is not Restriction.NONE
-
-        from pymbolic.primitives import Variable
-        if self.restriction == Restriction.POSITIVE:
-            return Variable(name_unit_outer_normal())
-        if self.restriction == Restriction.NEGATIVE:
-            # It is highly unnatural to have this generator function,
-            # but I do run into subtle trouble with return -1*outer
-            # as the indexing into the normal happens only later.
-            # Not investing more time into this cornercase right now.
-            return Variable(name_unit_inner_normal())
-
-    # TODO This one was just copied over so, it might need some tweaking
-    def facet_area(self, o):
-        from dune.perftool.pdelab.geometry import name_facetarea
-        return Variable(name_facetarea())
-
-    def quadrature_weight(self, o):
-        from dune.perftool.pdelab.quadrature import name_quadrature_weight
-        return Variable(name_quadrature_weight())
-
-    def jacobian_determinant(self, o):
-        return Variable(name_jacobian_determinant())
-
-    def jacobian_inverse(self, o):
-        restriction = self.restriction
-        if self.measure == 'exterior_facet':
-            restriction = Restriction.NEGATIVE
-
-        self.transpose_necessary = True
-        return Variable(name_jacobian_inverse_transposed(restriction))
-
-    def jacobian(self, o):
-        raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
-
-    def facet_jacobian_determinant(self, o):
-        return Variable(name_facet_jacobian_determinant())
-
-
 @iname
 def _dimension_iname(context, count):
     if context:

python/dune/perftool/pdelab/localoperator.py

@@ -257,7 +257,7 @@ def generate_kernel(integrals):
         accterms = split_into_accumulation_terms(integrand)
 
         # Get a transformer instance for this kernel
-        from dune.perftool.loopy.transformer import UFL2LoopyVisitor
+        from dune.perftool.ufl.visitor import UFL2LoopyVisitor
         visitor = UFL2LoopyVisitor(measure, subdomain_id, dimension_indices)
 
         # Iterate over the terms and generate a kernel

python/dune/perftool/pdelab/signatures.py

@@ -9,11 +9,12 @@ from dune.perftool.pdelab.argument import (name_accumulation_variable,
                                            name_coefficientcontainer,
                                            type_coefficientcontainer,
                                            name_applycontainer,
-                                           name_testfunctionspace,
-                                           type_testfunctionspace,
-                                           name_trialfunctionspace,
-                                           type_trialfunctionspace,
                                            )
+from dune.perftool.pdelab.spaces import (name_testfunctionspace,
+                                         type_testfunctionspace,
+                                         name_trialfunctionspace,
+                                         type_trialfunctionspace,
+                                         )
 
 
 def alpha_volume_signature():

python/dune/perftool/pdelab/spaces.py

+""" Generator functions for PDELab local/grid function spaces etc. """
+
+from dune.perftool.generation import (domain,
+                                      function_mangler,
+                                      generator_factory,
+                                      include_file,
+                                      preamble,
+                                      )
+from dune.perftool.pdelab import restricted_name
+
+from loopy import CallMangleInfo
+from loopy.symbolic import FunctionIdentifier
+from loopy.types import NumpyType
+
+from pymbolic.primitives import Variable
+
+import numpy
+
+
+class LFSChild(FunctionIdentifier):
+    def __init__(self, lfs):
+        self.lfs = lfs
+
+    def __getinitargs__(self):
+        return (self.lfs,)
+
+    @property
+    def name(self):
+        return '{}.child'.format(self.lfs)
+
+
+@function_mangler
+def lfs_child_mangler(target, func, dtypes):
+    if isinstance(func, LFSChild):
+        return CallMangleInfo(func.name, (NumpyType(str),), (NumpyType(numpy.int32),))
+
+
+@preamble
+def define_lfs_bound(lfs, bound):
+    return 'auto {} = {}.size();'.format(bound, lfs)
+
+
+def name_lfs_bound(lfs):
+    # LFS might either be given by an UFL element or by a string describig its name
+    from ufl import FiniteElementBase
+    if isinstance(lfs, FiniteElementBase):
+        return name_lfs_bound(name_lfs(lfs))
+
+    bound = '{}_size'.format(lfs)
+    define_lfs_bound(lfs, bound)
+    return bound
+
+
+@preamble
+def using_indices():
+    return "using namespace Dune::Indices;"
+
+
+@preamble
+def define_lfs(name, father, child):
+    using_indices()
+    return "auto {} = child({}, _{});".format(name, father, child)
+
+
+def lfs_child(lfs, children, shape=None, symmetry=False):
+    from pymbolic.primitives import Call, Product, Sum
+    # Old pre-TensorElement implementation kept for comaptibility
+    if shape is None:
+        indices = (Variable(children[0]),)
+    else:
+        if symmetry and len(children) == 2:
+            # I do not want to think about tensors of rank > 2 right now
+            i, j = children
+            if i > j:
+                j, i = i, j
+            i = Variable(i)
+            j = Variable(j)
+            n = len(children)
+
+            indices = (Sum((Product((n - 1, i)), Product((.5, i, 1 - i)), j)),)
+        else:
+            children = tuple(Variable(c) for c in children)
+            indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),)
+
+    return Call(LFSChild(lfs), indices)
+
+
+@generator_factory(cache_key_generator=lambda e, r, **kw: (e, r))
+def name_leaf_lfs(leaf_element, restriction, val=None):
+    """ This function just caches leaf lfs names based on the
+    element. The resulting local function spaces are useful only
+    for size information. OTOH, they are available with just the
+    leaf element available (as seen in basis evaluation).
+    """
+    # This generator function should be prepoluted by the lfs tree traversal,
+    # so val should always be None when we actually want the result
+    assert val
+
+    return val
+
+
+@generator_factory(cache_key_generator=lambda e, r, c, **kw: (e, r, c))
+def name_lfs(element, restriction, component, prefix=None):
+    # Omitting the prefix is only valid upon a second call, which will
+    # result in a cache hit.
+    assert prefix
+
+    def _name_lfs(prefix, component):
+        name = prefix
+        if len(component) > 0:
+            name = name + '_' + '_'.join(str(i) for i in component)
+        return name
+
+    name = _name_lfs(prefix, component)
+    if len(component) > 0:
+        father = _name_lfs(prefix, component[:-1])
+        # If this localfunction space is the child of another one, trigger
+        # the extraction preamble. Necessary before going into recursion
+        # for having the correct (top-down) order of preambles
+        define_lfs(name, father, component[-1])
+
+    # Recurse into the given element to define all other local function spaces!
+    from ufl import MixedElement
+    from ufl.functionview import select_subelement
+    from ufl.classes import FixedIndex
+    subel = select_subelement(element, component)
+    if isinstance(subel, MixedElement):
+        for i in range(subel.num_sub_elements()):
+            name_lfs(element, restriction, component + (FixedIndex(i),), prefix=prefix)
+
+    # Cache the name for the subelement
+    name_leaf_lfs(subel, restriction, val=name)
+
+    # Now return the prefix!
+    return name
+
+
+@generator_factory(cache_key_generator=lambda e, **kw: e)
+def type_gfs(element, basetype=None, index_stack=None):
+    # Omitting basetype and index_stack is only valid upon a second call,
+    # which will result in a cache hit.
+    assert basetype
+    assert index_stack is not None
+
+    # Additionally, element is expected to be a ufl finite element
+    from ufl import FiniteElementBase
+    assert isinstance(element, FiniteElementBase)
+
+    # Recurse into the given element to define all other local function spaces!
+    from ufl import MixedElement
+    from ufl.classes import FixedIndex
+    if isinstance(element, MixedElement):
+        for i, subelem in enumerate(element.sub_elements()):
+            type_gfs(subelem, basetype=basetype, index_stack=index_stack + (FixedIndex(i),))
+
+    if len(index_stack) == 0:
+        return basetype
+    else:
+        include_file("dune/typetree/childextraction.hh", filetag="operatorfile")
+        return 'Dune::TypeTree::Child<{},{}>'.format(basetype, ','.join(str(i) for i in index_stack))
+
+
+def traverse_lfs_tree(arg):
+    from dune.perftool.ufl.modified_terminals import ModifiedArgumentDescriptor
+    assert isinstance(arg, ModifiedArgumentDescriptor)
+
+    # First we need to determine the basename as given in the signature of
+    # this kernel method!
+    lfs_basename = name_argumentspace(arg)
+    from dune.perftool.pdelab.localoperator import lop_template_gfs
+    gfs_basename = lop_template_gfs(arg)
+
+    # Now start recursively extracting local function spaces and fill the cache with
+    # all those values. That way we can later get a correct local function space with
+    # just the ufl finite element.
+    from ufl.classes import MultiIndex
+    name_lfs(arg.argexpr.ufl_element(), arg.restriction, MultiIndex(()), prefix=lfs_basename)
+    type_gfs(arg.argexpr.ufl_element(), basetype=gfs_basename, index_stack=())
+
+
+@generator_factory(item_tags=("iname",), cache_key_generator=lambda e, r, c: (e, c))
+def _lfs_iname(element, restriction, context):
+    lfs = name_leaf_lfs(element, restriction)
+    bound = name_lfs_bound(lfs)
+
+    name = "{}_{}_index".format(lfs, context)
+    domain(name, bound)
+
+    return name
+
+
+def lfs_iname(element, restriction, count=None, context=''):
+    """ Get the iname to iterate over the local function space given by element
+
+    Arguments:
+    ----------
+    element: ufl.FiniteElementBase
+        The finite element this local function space belongs to
+    argcount: int
+        Use to realize double nesting in case of jacobians
+    context: str
+        Some generation methods will require you to duplicate an iname for
+        a given purpose, see the 'Loops and dependencies' of the loopy docs:
+        https://documen.tician.de/loopy/tutorial.html#loops-and-dependencies
+    """
+    assert not ((context == '') and (count is None))
+    if count is not None:
+        if context != '':
+            context = "{}_{}".format(count, context)
+        else:
+            context = str(count)
+    return _lfs_iname(element, restriction, context)
+
+
+def name_testfunctionspace(restriction):
+    return restricted_name("lfsv", restriction)
+
+
+def name_trialfunctionspace(restriction):
+    return restricted_name("lfsu", restriction)
+
+
+def name_argumentspace(ma):
+    from ufl.classes import Argument, Coefficient
+    if isinstance(ma.argexpr, Argument):
+        if ma.argexpr.number() == 0:
+            return name_testfunctionspace(ma.restriction)
+        if ma.argexpr.number() == 1:
+            return name_trialfunctionspace(ma.restriction)
+    if isinstance(ma.argexpr, Coefficient):
+        # Index 0 is reserved for trialfunction, index 1 is reserved for jacobian apply function
+        assert ma.argexpr.count() < 2
+        return name_trialfunctionspace(ma.restriction)
+    # We should never encounter an argument other than 0 or 1
+    assert False
+
+
+def type_testfunctionspace():
+    return "LFSV"
+
+
+def type_trialfunctionspace():
+    return "LFSU"

python/dune/perftool/pymbolic/uflmapper.py

-""" A multi function mapping UFL expressions to pymbolic expressions
-
-This mapper only defines handler for those UFL expression types that have an
-equivalent in pymbolic (mainly algebraic operations).
-
-This is mainly intended as a base class for anything mapping ufl to pymbolic.
-"""
-
-from __future__ import absolute_import
-from ufl.algorithms import MultiFunction
-from pymbolic.primitives import Product, Quotient, Subscript, Sum, Variable, Call
-
-# The constructed pymbolic expressions use n-ary operators instead of ufls binary operators
-from dune.perftool.ufl.flatoperators import get_operands
-
-
-class UFL2PymbolicMapper(MultiFunction):
-    def __init__(self):
-        super(UFL2PymbolicMapper, self).__init__()
-
-    call = MultiFunction.__call__
-
-    def product(self, o):
-        return Product(tuple(self.call(op) for op in get_operands(o)))
-
-    def float_value(self, o):
-        return o.value()
-
-    def int_value(self, o):
-        return o.value()
-
-    def division(self, o):
-        assert len(o.ufl_operands) == 2
-
-        return Quotient(self.call(o.ufl_operands[0]), self.call(o.ufl_operands[1]))
-
-    def sum(self, o):
-        return Sum(tuple(self.call(op) for op in get_operands(o)))
-
-    def zero(self, o):
-        return 0
-
-    def abs(self, o):
-        from ufl.classes import JacobianDeterminant
-        if isinstance(o.ufl_operands[0], JacobianDeterminant):
-            return self.call(o.ufl_operands[0])
-        else:
-            return Call('abs', self.call(o.ufl_operands[0]))

python/dune/perftool/pymbolic/simplify.py → python/dune/perftool/sympy.py

-from __future__ import absolute_import
-
 from pymbolic.interop.sympy import SympyToPymbolicMapper, PymbolicToSympyMapper
 
 import pymbolic.primitives as prim

python/dune/perftool/loopy/transformer.py → python/dune/perftool/ufl/visitor.py

 """
-This is the module that contains the main transformation from ufl to loopy
-(with pdelab as the hardcoded generation target)
+This module defines the main visitor algorithm transforming ufl expressions
+to pymbolic and loopy.
 """
 
-from __future__ import absolute_import
-
 from dune.perftool import Restriction
 from dune.perftool.ufl.modified_terminals import ModifiedTerminalTracker
-from dune.perftool.pymbolic.uflmapper import UFL2PymbolicMapper
-from dune.perftool.pdelab.geometry import GeometryMapper
 from dune.perftool.generation import (domain,
                                       get_temporary_name,
                                       global_context,
@@ -19,17 +15,18 @@ from dune.perftool.generation import (domain,
                                       valuearg,
                                       )
 
-from dune.perftool.pdelab.basis import (lfs_iname,
-                                        name_leaf_lfs,
-                                        name_lfs,
-                                        name_lfs_bound,
-                                        traverse_lfs_tree,
-                                        )
+from dune.perftool.pdelab.spaces import (lfs_iname,
+                                         name_leaf_lfs,
+                                         name_lfs,
+                                         name_lfs_bound,
+                                         traverse_lfs_tree,
+                                         )
 from dune.perftool.pdelab.quadrature import quadrature_iname
 from pymbolic.primitives import Subscript, Variable
+from ufl.algorithms import MultiFunction
 
 
-class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapper):
+class UFL2LoopyVisitor(ModifiedTerminalTracker):
     def __init__(self, measure, subdomain_id, dimension_indices):
         # Some variables describing the integral measure of this integral
         self.measure = measure
@@ -39,6 +36,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
         # Call base class constructors
         super(UFL2LoopyVisitor, self).__init__()
 
+    # Allow recursion through self.call(..)
+    call = MultiFunction.__call__
+
     def __call__(self, o):
         # Reset some state variables that are reinitialized for each accumulation term
         self.argshape = 0
@@ -59,6 +59,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
             if self.measure == 'exterior_facet':
                 ma.restriction = Restriction.NEGATIVE
 
+            from dune.perftool.pdelab.spaces import traverse_lfs_tree
             traverse_lfs_tree(ma)
 
         # Determine the rank of the term
@@ -166,6 +167,11 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
                     predicates=predicates
                     )
 
+    #
+    # Form argument/coefficients handlers:
+    # This is where the actual domain specific work happens
+    #
+
     def argument(self, o):
         # Correct the restriction on boundary integrals
         restriction = self.restriction
@@ -199,6 +205,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
             raise ValueError("Gradients should have been transformed to reference gradients!!!")
 
         # Have the issued instruction depend on the iname for this localfunction space
+        from dune.perftool.pdelab.spaces import lfs_iname
         iname = lfs_iname(leaf_element, restriction, o.number())
         self.inames.append(iname)
 
@@ -260,6 +267,10 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
             # And return a symbol
             return Variable(name)
 
+    #
+    # Handlers for all indexing related stuff
+    #
+
     def indexed(self, o):
         """ Although indexed is something that we want to handle in the pymbolic mapper we
         need to handle it here, as we have to give special treatment to indices into vector
@@ -334,3 +345,89 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
 
     def index(self, o):
         return self._index_or_fixed_index(o)
+
+    #
+    # Handlers for arithmetic operators and functions
+    # Those handlers would be valid in any code going from UFL to pymbolic
+    #
+
+    def product(self, o):
+        from dune.perftool.ufl.flatoperators import get_operands
+        from pymbolic.primitives import Product
+        return Product(tuple(self.call(op) for op in get_operands(o)))
+
+    def float_value(self, o):
+        return o.value()
+
+    def int_value(self, o):
+        return o.value()
+
+    def division(self, o):
+        assert len(o.ufl_operands) == 2
+
+        from pymbolic.primitives import Quotient
+        return Quotient(self.call(o.ufl_operands[0]), self.call(o.ufl_operands[1]))
+
+    def sum(self, o):
+        from dune.perftool.ufl.flatoperators import get_operands
+        from pymbolic.primitives import Sum
+        return Sum(tuple(self.call(op) for op in get_operands(o)))
+
+    def zero(self, o):
+        return 0
+
+    def abs(self, o):
+        from ufl.classes import JacobianDeterminant
+        if isinstance(o.ufl_operands[0], JacobianDeterminant):
+            return self.call(o.ufl_operands[0])
+        else:
+            from pymbolic.primitives import Call
+            return Call('abs', self.call(o.ufl_operands[0]))
+
+    #
+    # Handlers for geometric quantities
+    #
+
+    def facet_normal(self, o):
+        # The normal must be restricted to be well-defined
+        assert self.restriction is not Restriction.NONE
+
+        from pymbolic.primitives import Variable
+        if self.restriction == Restriction.POSITIVE:
+            from dune.perftool.pdelab.geometry import name_unit_outer_normal
+            return Variable(name_unit_outer_normal())
+        if self.restriction == Restriction.NEGATIVE:
+            # It is highly unnatural to have this generator function,
+            # but I do run into subtle trouble with return -1*outer
+            # as the indexing into the normal happens only later.
+            # Not investing more time into this cornercase right now.
+            from dune.perftool.pdelab.geometry import name_unit_inner_normal
+            return Variable(name_unit_inner_normal())
+
+    def facet_area(self, o):
+        from dune.perftool.pdelab.geometry import name_facetarea
+        return Variable(name_facetarea())
+
+    def quadrature_weight(self, o):
+        from dune.perftool.pdelab.quadrature import name_quadrature_weight
+        return Variable(name_quadrature_weight())
+
+    def jacobian_determinant(self, o):
+        from dune.perftool.pdelab.geometry import name_jacobian_determinant
+        return Variable(name_jacobian_determinant())
+
+    def jacobian_inverse(self, o):
+        restriction = self.restriction
+        if self.measure == 'exterior_facet':
+            restriction = Restriction.NEGATIVE
+
+        self.transpose_necessary = True
+        from dune.perftool.pdelab.geometry import name_jacobian_inverse_transposed
+        return Variable(name_jacobian_inverse_transposed(restriction))
+
+    def jacobian(self, o):
+        raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
+
+    def facet_jacobian_determinant(self, o):
+        from dune.perftool.pdelab.geometry import name_facet_jacobian_determinant
+        return Variable(name_facet_jacobian_determinant())
\ No newline at end of file