Merge branch 'feature/tensor-product-element-2' into 'master'

TensorProductElement

See merge request !174

patches/apply_patches.sh

@@ -12,3 +12,7 @@ popd
 pushd python/ufl
 git apply ../../patches/ufl/conditional-uflid.patch
 popd
+
+pushd python/ufl
+git apply ../../patches/ufl/tensor-product-element.patch
+popd

patches/ufl/tensor-product-element.patch

+commit f87dcd18d765b0200808b79b2e7374f82a0c6199
+Author: René Heß <rene.hess@iwr.uni-heidelberg.de>
+Date:   Tue Aug 29 14:56:17 2017 +0200
+
+    Patch for TensorProductElements
+
+diff --git a/ufl/algorithms/compute_form_data.py b/ufl/algorithms/compute_form_data.py
+index 3388bbfc..1cef3924 100644
+--- a/ufl/algorithms/compute_form_data.py
++++ b/ufl/algorithms/compute_form_data.py
+@@ -56,7 +56,7 @@ def _auto_select_degree(elements):
+     """
+     # Use max degree of all elements, at least 1 (to work with
+     # Lagrange elements)
+-    return max({e.degree() for e in elements} - {None} | {1})
++    return max({e.degree() if not isinstance(e.degree(), tuple) else max(e.degree()) for e in elements} - {None} | {1})
+ 
+ 
+ def _compute_element_mapping(form):

python/dune/perftool/generation/loopy.py

@@ -50,10 +50,21 @@ def valuearg(name, **kw):
 
 @generator_factory(item_tags=("domain",), context_tags="kernel")
 def domain(iname, shape):
-    if isinstance(iname, tuple):
+    if isinstance(iname, tuple) and isinstance(shape, tuple):
+        assert(len(iname) == len(shape))
+        condition = ""
+        for index, (i, s) in enumerate(zip(iname, shape)):
+            if index > 0:
+                condition += " and "
+            else:
+                condition += " "
+            condition += "0<={}<{}".format(i, s)
         iname = ",".join(iname)
-    if isinstance(shape, str):
-        valuearg(shape)
+        return "{{ [{}] : {} }}".format(iname, condition)
+    elif isinstance(iname, tuple):
+        iname = ",".join(iname)
+        if isinstance(shape, str):
+            valuearg(shape)
     return "{{ [{0}] : 0<={0}<{1} }}".format(iname, shape)
 
 

python/dune/perftool/loopy/transformations/vectorize_quad.py

@@ -136,6 +136,7 @@ def _vectorize_quadrature_loop(knl, inames, suffix):
             # Do precomputation of the quantity
             prec_quantity = "{}_precomputed".format(quantity)
             prec_quantities.append(prec_quantity)
+
             knl = lp.precompute(knl, subst_name, inames,
                                 temporary_name=prec_quantity,
                                 )

python/dune/perftool/options.py

@@ -63,8 +63,8 @@ class PerftoolOptionsArray(ImmutableRecord):
 
     # Arguments that are mainly to be set by logic depending on other options
     max_vector_width = PerftoolOption(default=256, helpstr=None)
-    unroll_dimension_loops = PerftoolOption(default=False, helpstr="whether loops over the gemetric dimension should be unrolled.")
-    precompute_quadrature_info = PerftoolOption(default=True, helpstr="whether loops over the gemetric dimension should be unrolled.")
+    unroll_dimension_loops = PerftoolOption(default=False, helpstr="whether loops over the geometric dimension should be unrolled")
+    precompute_quadrature_info = PerftoolOption(default=True, helpstr="compute quadrature points and weights in the constructor of the local operator")
     blockstructured = PerftoolOption(default=False, helpstr="Use block structure")
     number_of_blocks = PerftoolOption(default=1, helpstr="Number of sub blocks in one direction")
 

python/dune/perftool/pdelab/driver/__init__.py

@@ -136,7 +136,7 @@ def isDG(fem):
 
 
 def FEM_name_mangling(fem):
-    from ufl import MixedElement, VectorElement, FiniteElement, TensorElement
+    from ufl import MixedElement, VectorElement, FiniteElement, TensorElement, TensorProductElement
     if isinstance(fem, VectorElement):
         return FEM_name_mangling(fem.sub_elements()[0]) + "_" + str(fem.num_sub_elements())
     if isinstance(fem, TensorElement):
@@ -155,6 +155,14 @@ def FEM_name_mangling(fem):
             return "Q" + str(fem._degree)
         if isDG(fem):
             return "DG" + str(fem._degree)
+    if isinstance(fem, TensorProductElement):
+        assert(len(set(subel._short_name for subel in fem.sub_elements())) == 1)
+
+        if isLagrange(fem.sub_elements()[0]):
+            return "TensorQ" + '_'.join(map(str, fem._degree))
+        if isDG(fem.sub_elements()[0]):
+            return "TensorDG" + '_'.join(map(str, fem._degree))
+        raise NotImplementedError("fem name mangling")
 
     raise NotImplementedError("FEM NAME MANGLING")
 

python/dune/perftool/pdelab/driver/constraints.py

@@ -14,7 +14,7 @@ from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
                                                            preprocess_leaf_data,
                                                            )
 
-from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
 
 
 def name_assembled_constraints():
@@ -53,7 +53,7 @@ def name_bctype_function(element, is_dirichlet):
         define_composite_bctype_function(element, is_dirichlet, name, tuple(childs))
         return name
     else:
-        assert isinstance(element, FiniteElement)
+        assert isinstance(element, (FiniteElement, TensorProductElement))
         name = "{}_bctype".format(FEM_name_mangling(element).lower())
         define_bctype_function(element, is_dirichlet[0], name)
         return name

python/dune/perftool/pdelab/driver/gridfunctionspace.py

@@ -7,6 +7,7 @@ from dune.perftool.pdelab.driver import (FEM_name_mangling,
                                          get_dimension,
                                          get_test_element,
                                          get_trial_element,
+                                         isLagrange,
                                          isDG,
                                          isPk,
                                          isQk,
@@ -16,7 +17,7 @@ from dune.perftool.pdelab.driver import (FEM_name_mangling,
                                          preprocess_leaf_data,
                                          )
 
-from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
 
 
 @preamble
@@ -124,24 +125,52 @@ def typedef_fem(element, name):
     df = type_domainfield()
     r = type_range()
     dim = get_dimension()
+
     if get_option("blockstructured"):
         include_file("dune/perftool/blockstructured/blockstructuredqkfem.hh", filetag="driver")
         degree = element.degree() * get_option("number_of_blocks")
-        return "using {} = Dune::PDELab::BlockstructuredQkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, degree)
-    if isQk(element):
+        return "using {} = Dune::PDELab::BlockstructuredQkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+            .format(name, gv, df, r, degree)
+
+    if isinstance(element, TensorProductElement):
+        # Only allow TensorProductElements where all subelements are
+        # of the same type ('CG' or 'DG')
+        assert(len(set(subel._short_name for subel in element.sub_elements())) == 1)
+
+        # Anisotropic degree is not yet supported in Dune
+        degrees = element.degree()
+        for deg in degrees:
+            assert (deg == degrees[0])
+
+        # TensorProductElements have Qk structure -> no Pk
+        if isLagrange(element.sub_elements()[0]):
+            include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
+            return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+                .format(name, gv, df, r, degrees[0])
+        elif isDG(element.sub_elements()[0]):
+            include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
+            # TODO allow switching the basis here!
+            return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;" \
+                .format(name, df, r, degrees[0], dim)
+        raise NotImplementedError("FEM not implemented in dune-perftool")
+    elif isQk(element):
         include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
-        return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
-    if isPk(element):
+        return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+            .format(name, gv, df, r, element.degree())
+    elif isPk(element):
         include_file("dune/pdelab/finiteelementmap/pkfem.hh", filetag="driver")
-        return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
-    if isDG(element):
+        return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+            .format(name, gv, df, r, element.degree())
+    elif isDG(element):
         if isQuadrilateral(element.cell()):
             include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
             # TODO allow switching the basis here!
-            return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, df, r, element.degree(), dim)
+            return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;" \
+                .format(name, df, r, element.degree(), dim)
         if isSimplical(element.cell()):
             include_file("dune/pdelab/finiteelementmap/opbfem.hh", filetag="driver")
-            return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;".format(name, df, r, element.degree(), dim)
+            return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;" \
+                .format(name, df, r, element.degree(), dim)
         raise NotImplementedError("Geometry type not known in code generation")
 
     raise NotImplementedError("FEM not implemented in dune-perftool")
@@ -157,15 +186,26 @@ def type_fem(element):
 def define_fem(element, name):
     femtype = type_fem(element)
     from dune.perftool.pdelab.driver import isDG
-    if isDG(element):
+    if isinstance(element, TensorProductElement):
+        # Only allow TensorProductElements where all subelements are
+        # of the same type ('CG' or 'DG')
+        assert(len(set(subel._short_name for subel in element.sub_elements())) == 1)
+        if isDG(element.sub_elements()[0]):
+            return "{} {};".format(femtype, name)
+        else:
+            assert(isLagrange(element.sub_elements()[0]))
+            gv = name_leafview()
+            return "{} {}({});".format(femtype, name, gv)
+    elif isDG(element):
         return "{} {};".format(femtype, name)
     else:
+        assert(isLagrange(element))
         gv = name_leafview()
         return "{} {}({});".format(femtype, name, gv)
 
 
 def name_fem(element):
-    assert isinstance(element, FiniteElement)
+    assert isinstance(element, (FiniteElement, TensorProductElement))
     name = "{}_fem".format(FEM_name_mangling(element).lower())
     define_fem(element, name)
     return name
@@ -192,7 +232,7 @@ def name_gfs(element, is_dirichlet, treepath=()):
                                        "_".join(str(t) for t in treepath))
         define_power_gfs(element, is_dirichlet, name, subgfs)
         return name
-    if isinstance(element, MixedElement):
+    elif isinstance(element, MixedElement):
         k = 0
         subgfs = []
         for i, subel in enumerate(element.sub_elements()):
@@ -203,7 +243,7 @@ def name_gfs(element, is_dirichlet, treepath=()):
         define_composite_gfs(element, is_dirichlet, name, tuple(subgfs))
         return name
     else:
-        assert isinstance(element, FiniteElement)
+        assert isinstance(element, (FiniteElement, TensorProductElement))
         name = "{}{}_gfs_{}".format(FEM_name_mangling(element).lower(),
                                     "_dirichlet" if is_dirichlet[0] else "",
                                     "_".join(str(t) for t in treepath))
@@ -230,7 +270,7 @@ def type_gfs(element, is_dirichlet):
         name = "{}_POW{}GFS".format(subgfs, element.num_sub_elements())
         typedef_power_gfs(element, is_dirichlet, name, subgfs)
         return name
-    if isinstance(element, MixedElement):
+    elif isinstance(element, MixedElement):
         k = 0
         subgfs = []
         for subel in element.sub_elements():
@@ -240,7 +280,7 @@ def type_gfs(element, is_dirichlet):
         typedef_composite_gfs(element, name, tuple(subgfs))
         return name
     else:
-        assert isinstance(element, FiniteElement)
+        assert isinstance(element, (FiniteElement, TensorProductElement))
         name = "{}{}_GFS".format(FEM_name_mangling(element).upper(),
                                  "_dirichlet" if is_dirichlet[0] else "",
                                  )

python/dune/perftool/pdelab/driver/interpolate.py

@@ -15,7 +15,7 @@ from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
                                                            )
 from dune.perftool.pdelab.driver.gridoperator import (name_parameters,)
 
-from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
 
 
 def _do_interpolate(dirichlet):
@@ -55,7 +55,7 @@ def name_boundary_function(element, func):
         define_compositegfs_parameterfunction(name, tuple(childs))
         return name
     else:
-        assert isinstance(element, FiniteElement)
+        assert isinstance(element, (FiniteElement, TensorProductElement))
         name = get_counted_variable("func")
         define_boundary_function(name, func[0])
         return name

python/dune/perftool/pdelab/driver/vtk.py

@@ -41,7 +41,10 @@ def type_vtkwriter():
 @preamble
 def define_subsamplinglevel(name):
     ini = name_initree()
-    return "int {} = {}.get<int>(\"vtk.subsamplinglevel\", {});".format(name, ini, max(get_trial_element().degree() - 1, 0))
+    degree = get_trial_element().degree()
+    if isinstance(degree, tuple):
+        degree = max(degree)
+    return "int {} = {}.get<int>(\"vtk.subsamplinglevel\", {});".format(name, ini, max(degree - 1, 0))
 
 
 def name_subsamplinglevel():

python/dune/perftool/pdelab/localoperator.py

@@ -775,9 +775,9 @@ def generate_localoperator_kernels(formdata, data):
 
                 operator_kernels[(measure, 'residual')] = kernel
 
-    logger.info("generate_localoperator_kernels: create jacobian methods")
     # Generate the necessary jacobian methods
     if not get_option("numerical_jacobian"):
+        logger.info("generate_localoperator_kernels: create jacobian methods")
         from ufl import derivative
         jacform = derivative(formdata.original_form, formdata.original_form.coefficients()[0])
 

python/dune/perftool/pdelab/quadrature.py

@@ -189,9 +189,25 @@ def _estimate_quadrature_order():
     # Estimate polynomial degree of integrals of current type (eg 'Cell')
     integral_type = get_global_context_value("integral_type")
     integrals = form.integrals_by_type(integral_type)
-    polynomial_degree = 0
+
+    # Degree could be a tuple (for TensorProductElements)
+    degree = integrals[0].metadata()['estimated_polynomial_degree']
+    if isinstance(degree, int):
+        degree = (degree,)
+    polynomial_degree = [0, ] * len(degree)
+
     for i in integrals:
-        polynomial_degree = max(polynomial_degree, i.metadata()['estimated_polynomial_degree'])
+        degree = i.metadata()['estimated_polynomial_degree']
+        if isinstance(degree, int):
+            degree = [degree, ]
+        assert(len(polynomial_degree) == len(degree))
+        for i in range(len(polynomial_degree)):
+            polynomial_degree[i] = max(polynomial_degree[i], degree[i])
+
+    # Return either a tuple or an int
+    polynomial_degree = tuple(polynomial_degree)
+    if len(polynomial_degree) == 1:
+        polynomial_degree = polynomial_degree[0]
 
     return polynomial_degree
 
@@ -199,14 +215,31 @@ def _estimate_quadrature_order():
 def quadrature_order():
     """Get quadrature order
 
-    Note: In PDELab quadrature order m means that integration of
-    polynomials of degree m is excat.
+    Notes:
+
+    - In PDELab quadrature order m means that integration of
+      polynomials of degree m is exact.
+
+    - If you use sum factorization and TensorProductElement it is
+      possible to use a different quadrature_order per direction.
     """
     if get_option("quadrature_order"):
-        quadrature_order = int(get_option("quadrature_order"))
+        quadrature_order = tuple(map(int, get_option("quadrature_order").split(',')))
     else:
         quadrature_order = _estimate_quadrature_order()
 
+    # TensorProductElements can have different quadrature order for
+    # every direction so quadrature_order may be a tuple. If we do not
+    # use TensorProductElements we want to return an int.
+    if isinstance(quadrature_order, tuple):
+        if len(quadrature_order) == 1:
+            quadrature_order = quadrature_order[0]
+    if isinstance(quadrature_order, tuple):
+        if not get_option('sumfact'):
+            raise NotImplementedError("Different quadrature order per direction is only implemented for kernels using sum factorization.")
+        from dune.perftool.pdelab.geometry import world_dimension
+        assert(len(quadrature_order) == world_dimension())
+
     return quadrature_order
 
 

python/dune/perftool/pdelab/spaces.py

@@ -6,6 +6,7 @@ from dune.perftool.generation import (class_member,
                                       generator_factory,
                                       include_file,
                                       preamble,
+                                      valuearg,
                                       )
 from dune.perftool.pdelab.restriction import restricted_name
 from dune.perftool.ufl.modified_terminals import Restriction
@@ -73,6 +74,7 @@ def define_lfs(name, father, child):
 @preamble
 def define_lfs_size(lfs, element, restriction):
     name = name_lfs_bound(lfs)
+    valuearg(name, dtype=numpy.int32)
     return "auto {} = {}.size();".format(name, lfs)
 
 

python/dune/perftool/sumfact/accumulation.py

@@ -41,6 +41,7 @@ import loopy as lp
 import numpy as np
 import pymbolic.primitives as prim
 import ufl.classes as uc
+from ufl import FiniteElement, MixedElement, TensorProductElement
 
 
 @iname
@@ -138,8 +139,7 @@ def get_accumulation_info(expr, visitor):
 
 
 def _get_childs(element):
-    from ufl import FiniteElement
-    if isinstance(element, FiniteElement):
+    if isinstance(element, (FiniteElement, TensorProductElement)):
         yield (0, element)
     else:
         for i in range(element.value_size()):
@@ -199,10 +199,18 @@ def generate_accumulation_instruction(expr, visitor):
     test_info = visitor.test_info
     trial_info = visitor.trial_info
 
+    # Number of basis functions per direction
     leaf_element = test_info.element
-    if leaf_element.num_sub_elements() > 0:
+    from ufl import MixedElement
+    if isinstance(leaf_element, MixedElement):
         leaf_element = leaf_element.extract_component(test_info.element_index)[1]
-    basis_size = leaf_element.degree() + 1
+    degree = leaf_element._degree
+    if isinstance(degree, int):
+        degree = (degree,) * dim
+    basis_size = tuple(deg + 1 for deg in degree)
+
+    # Anisotropic finite elements are not (yet) supported by Dune
+    assert(size == basis_size[0] for size in basis_size)
 
     from dune.perftool.pdelab.localoperator import boundary_predicates
     predicates = boundary_predicates(expr,
@@ -316,7 +324,7 @@ def generate_accumulation_instruction(expr, visitor):
 
     # Collect the lfs and lfs indices for the accumulate call
     test_lfs.index = flatten_index(tuple(prim.Variable(i) for i in inames),
-                                   (basis_size,) * dim,
+                                   basis_size,
                                    order="f"
                                    )
 
@@ -327,7 +335,7 @@ def generate_accumulation_instruction(expr, visitor):
         from dune.perftool.sumfact.basis import _basis_functions_per_direction
         ansatz_lfs.index = flatten_index(tuple(prim.Variable(jacobian_inames[i])
                                                for i in range(world_dimension())),
-                                         (_basis_functions_per_direction(trial_leaf_element),) * dim,
+                                         _basis_functions_per_direction(trial_leaf_element),
                                          order="f"
                                          )
 

python/dune/perftool/sumfact/basis.py

@@ -42,7 +42,7 @@ from dune.perftool.tools import maybe_wrap_subscript
 from dune.perftool.pdelab.basis import shape_as_pymbolic
 from dune.perftool.sumfact.accumulation import sumfact_iname
 
-from ufl import VectorElement, TensorElement
+from ufl import MixedElement, VectorElement, TensorElement, TensorProductElement
 
 from pytools import product, ImmutableRecord
 
@@ -96,20 +96,30 @@ class LFSSumfactKernelInput(SumfactKernelInputBase, ImmutableRecord):
 
 def _basis_functions_per_direction(element):
     """Number of basis functions per direction """
-    from ufl import FiniteElement
-    assert isinstance(element, FiniteElement)
-    return element.degree() + 1
+    from ufl import FiniteElement, TensorProductElement
+    assert isinstance(element, (FiniteElement, TensorProductElement))
+    degree = element.degree()
+    if isinstance(degree, int):
+        degree = (degree,) * world_dimension()
+
+    basis_size = tuple(deg + 1 for deg in degree)
+
+    # Anisotropic finite elements are not (yet) supported by Dune
+    for size in basis_size:
+        assert(size == basis_size[0])
+
+    return basis_size
 
 
 @kernel_cached
 def pymbolic_coefficient_gradient(element, restriction, index, coeff_func, visitor_indices):
     sub_element = element
     grad_index = visitor_indices[0]
-    if element.num_sub_elements() > 0:
+    if isinstance(element, MixedElement):
         sub_element = element.extract_component(index)[1]
 
     from ufl import FiniteElement
-    assert isinstance(sub_element, FiniteElement)
+    assert isinstance(sub_element, (FiniteElement, TensorProductElement))
 
     # Number of basis functions per direction
     basis_size = _basis_functions_per_direction(sub_element)
@@ -182,6 +192,7 @@ def pymbolic_coefficient(element, restriction, index, coeff_func, visitor_indice
 
 @iname
 def sumfact_lfs_iname(element, bound, dim):
+    assert(isinstance(bound, int))
     from dune.perftool.pdelab.driver import FEM_name_mangling
     name = "sumfac_lfs_{}_{}".format(FEM_name_mangling(element), dim)
     domain(name, bound)
@@ -190,13 +201,13 @@ def sumfact_lfs_iname(element, bound, dim):
 
 @backend(interface="lfs_inames", name="sumfact")
 def lfs_inames(element, restriction, number=1, context=''):
-    from ufl import FiniteElement
-    assert isinstance(element, FiniteElement)
+    from ufl import FiniteElement, TensorProductElement
+    assert isinstance(element, (FiniteElement, TensorProductElement))
     if number == 0:
         return ()
     else:
         dim = world_dimension()
-        return tuple(sumfact_lfs_iname(element, _basis_functions_per_direction(element), d) for d in range(dim))
+        return tuple(sumfact_lfs_iname(element, _basis_functions_per_direction(element)[d], d) for d in range(dim))
 
 
 @backend(interface="evaluate_basis", name="sumfact")
@@ -209,10 +220,18 @@ def evaluate_basis(element, name, restriction):
 
     # Collect the pairs of lfs/quad inames that are in use
     # On facets, the normal direction of the facet is excluded
-    tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element))
-    prod = tuple(tab.pymbolic((prim.Variable(i), prim.Variable(j)))
-                 for (i, j) in zip(quad_inames, tuple(iname for i, iname in enumerate(inames) if i != facedir))
-                 )
+    #
+    # If facedir is not none, the length of inames and quad_inames is
+    # different. For inames we want to skip the facedir direction, for
+    # quad_inames we need all entries. Thats the reason for the
+    # help_index.
+    basis_per_dir = _basis_functions_per_direction(element)
+    prod = ()
+    help_index = 0
+    for direction in range(len(inames)):
+        if direction != facedir:
+            prod = prod + (BasisTabulationMatrix(basis_size=basis_per_dir[direction]).pymbolic((prim.Variable(quad_inames[help_index]), prim.Variable(inames[direction]))),)
+            help_index += 1
 
     # Add the missing direction on facedirs by evaluating at either 0 or 1
     if facedir is not None:
@@ -263,7 +282,7 @@ def evaluate_reference_gradient(element, name, restriction, index):
     prod = []
     for i in range(dim):
         if i != facedir:
-            tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element),
+            tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element)[i],
                                         derivative=index == i)
             prod.append(tab.pymbolic((prim.Variable(quadinamemapping[i]), prim.Variable(inames[i]))))
 

python/dune/perftool/sumfact/quadrature.py

@@ -10,13 +10,16 @@ from dune.perftool.generation import (backend,
                                       )
 
 from dune.perftool.sumfact.tabulation import (quadrature_points_per_direction,
+                                              local_quadrature_points_per_direction,
                                               name_oned_quadrature_points,
                                               name_oned_quadrature_weights,
                                               )
 from dune.perftool.pdelab.argument import name_accumulation_variable
 from dune.perftool.pdelab.geometry import (local_dimension,
+                                           world_dimension,
                                            )
 from dune.perftool.options import get_option
+from dune.perftool.sumfact.switch import get_facedir
 
 from loopy import CallMangleInfo
 from loopy.symbolic import FunctionIdentifier
@@ -80,52 +83,59 @@ def quadrature_inames(element):
     if element is None:
         names = tuple("quad_{}".format(d) for d in range(local_dimension()))
     else:
-        from ufl import FiniteElement
-        assert isinstance(element, FiniteElement)
+        from ufl import FiniteElement, TensorProductElement
+        assert isinstance(element, (FiniteElement, TensorProductElement))
         from dune.perftool.pdelab.driver import FEM_name_mangling
         names = tuple("quad_{}_{}".format(FEM_name_mangling(element), d) for d in range(local_dimension()))
-    domain(names, quadrature_points_per_direction())
+
+    local_qps_per_dir = local_quadrature_points_per_direction()
+    domain(names, local_qps_per_dir)
     return names
 
 
 @iname(kernel="operator")
 def constructor_quad_iname(name, d, bound):
-    name = "{}_{}".format(name, d)
-    domain(name, quadrature_points_per_direction(), kernel="operator")
+    name = "{}_localdim{}".format(name, d)
+    domain(name, bound, kernel="operator")
     return name
 
 
-def constructor_quadrature_inames(dim, num1d):
-    name = "quadiname_dim{}_num{}".format(dim, num1d)
-    return tuple(constructor_quad_iname(name, d, quadrature_points_per_direction()) for d in range(local_dimension()))
+def constructor_quadrature_inames(dim):
+    local_qps_per_dir = local_quadrature_points_per_direction()
+    local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
+    name = "quadiname_dim{}_num{}".format(dim, local_qps_per_dir_str)
+    return tuple(constructor_quad_iname(name, d, local_qps_per_dir[d]) for d in range(local_dimension()))
 
 
-def define_recursive_quadrature_weight(visitor, name, dir):
+def define_recursive_quadrature_weight(visitor, name, direction):
     info = visitor.current_info[1]
     if info is None:
         element = None
     else:
         element = info.element
-    iname = quadrature_inames(element)[dir]
+    iname = quadrature_inames(element)[direction]
     temporary_variable(name, shape=(), shape_impl=())
-    instruction(expression=Product((recursive_quadrature_weight(dir + 1),
-                                    Subscript(Variable(name_oned_quadrature_weights()),
+
+    qps_per_dir = quadrature_points_per_direction()
+    qp_bound = qps_per_dir[direction]
+    instruction(expression=Product((recursive_quadrature_weight(direction + 1),
+                                    Subscript(Variable(name_oned_quadrature_weights(qp_bound)),
                                               (Variable(iname),)
                                               ))
                                    ),
                 assignee=Variable(name),
-                forced_iname_deps=frozenset(quadrature_inames()[dir:]),
+                forced_iname_deps=frozenset(quadrature_inames()[direction:]),
                 forced_iname_deps_is_final=True,
                 tags=frozenset({"quad"}),
                 )
 
 
-def recursive_quadrature_weight(visitor, dir=0):
-    if dir == local_dimension():
+def recursive_quadrature_weight(visitor, direction=0):
+    if direction == local_dimension():
         return pymbolic_base_weight()
     else:
-        name = 'weight_{}'.format(dir)
-        define_recursive_quadrature_weight(visitor, name, dir)
+        name = 'weight_{}'.format(direction)
+        define_recursive_quadrature_weight(visitor, name, direction)
         return Variable(name)
 
 
@@ -134,14 +144,16 @@ def quadrature_weight(visitor):
     if not get_option("precompute_quadrature_info"):
         return recursive_quadrature_weight(visitor)
 
+    # Quadrature points per (local) direction
     dim = local_dimension()
-    num1d = quadrature_points_per_direction()
-    name = "quad_weights_dim{}_num{}".format(dim, num1d)
+    local_qps_per_dir = local_quadrature_points_per_direction()
+    local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
+    name = "quad_weights_dim{}_num{}".format(dim, local_qps_per_dir_str)
 
     # Add a class member
     loopy_class_member(name,
                        dtype=np.float64,
-                       shape=(num1d,) * dim,
+                       shape=local_qps_per_dir,
                        classtag="operator",
                        dim_tags=",".join(["f"] * dim),
                        managed=True,
@@ -149,9 +161,12 @@ def quadrature_weight(visitor):
                        )
 
     # Precompute it in the constructor
-    inames = constructor_quadrature_inames(dim, num1d)
-    instruction(assignee=prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames)),
-                expression=prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights()), (prim.Variable(i),)) for i in inames)),
+    inames = constructor_quadrature_inames(dim)
+    assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
+    expression = prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights(local_qps_per_dir[index])),
+                                              (prim.Variable(iname),)) for index, iname in enumerate(inames)))
+    instruction(assignee=assignee,
+                expression=expression,
                 within_inames=frozenset(inames),
                 kernel="operator",
                 )
@@ -165,7 +180,11 @@ def quadrature_weight(visitor):
 
 
 def define_quadrature_position(name, index):
-    instruction(expression=Subscript(Variable(name_oned_quadrature_points()), (Variable(quadrature_inames()[index]),)),
+    qps_per_dir = quadrature_points_per_direction()
+    qp_bound = qps_per_dir[index]
+    expression = Subscript(Variable(name_oned_quadrature_points(qp_bound)),
+                           (Variable(quadrature_inames()[index]),))
+    instruction(expression=expression,
                 assignee=Subscript(Variable(name), (index,)),
                 forced_iname_deps=frozenset(quadrature_inames(element)),
                 forced_iname_deps_is_final=True,
@@ -184,12 +203,14 @@ def pymbolic_quadrature_position(index, visitor):
 
     assert isinstance(index, int)
     dim = local_dimension()
-    num1d = quadrature_points_per_direction()
-    name = "quad_points_dim{}_num{}_dir{}".format(dim, num1d, index)
+    qps_per_dir = quadrature_points_per_direction()
+    local_qps_per_dir = local_quadrature_points_per_direction()
+    local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
+    name = "quad_points_dim{}_num{}_dir{}".format(dim, local_qps_per_dir_str, index)
 
     loopy_class_member(name,
                        dtype=np.float64,
-                       shape=(num1d,) * dim,
+                       shape=local_qps_per_dir,
                        classtag="operator",
                        dim_tags=",".join(["f"] * dim),
                        managed=True,
@@ -197,9 +218,12 @@ def pymbolic_quadrature_position(index, visitor):
                        )
 
     # Precompute it in the constructor
-    inames = constructor_quadrature_inames(dim, num1d)
-    instruction(assignee=prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames)),
-                expression=Subscript(Variable(name_oned_quadrature_points()), (prim.Variable(inames[index]))),
+    inames = constructor_quadrature_inames(dim)
+    assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
+    expression = Subscript(Variable(name_oned_quadrature_points(qps_per_dir[index])),
+                           (prim.Variable(inames[index])))
+    instruction(assignee=assignee,
+                expression=expression,
                 within_inames=frozenset(inames),
                 kernel="operator",
                 )

python/dune/perftool/sumfact/symbolic.py

@@ -4,7 +4,7 @@ from dune.perftool.options import get_option
 from dune.perftool.generation import get_counted_variable
 from dune.perftool.pdelab.geometry import local_dimension, world_dimension
 from dune.perftool.sumfact.quadrature import quadrature_inames
-from dune.perftool.sumfact.tabulation import BasisTabulationMatrixArray
+from dune.perftool.sumfact.tabulation import BasisTabulationMatrixBase, BasisTabulationMatrixArray
 
 from pytools import ImmutableRecord, product
 
@@ -107,7 +107,6 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
         """
         # Assert the inputs!
         assert isinstance(matrix_sequence, tuple)
-        from dune.perftool.sumfact.tabulation import BasisTabulationMatrixBase
         assert all(isinstance(m, BasisTabulationMatrixBase) for m in matrix_sequence)
 
         assert stage in (1, 3)

python/dune/perftool/sumfact/tabulation.py

@@ -3,7 +3,7 @@ from dune.perftool.ufl.modified_terminals import Restriction
 from dune.perftool.options import get_option
 
 from dune.perftool.pdelab.argument import name_coefficientcontainer
-from dune.perftool.pdelab.geometry import world_dimension
+from dune.perftool.pdelab.geometry import world_dimension, local_dimension
 from dune.perftool.generation import (class_member,
                                       domain,
                                       function_mangler,
@@ -51,11 +51,13 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
                  slice_size=None,
                  slice_index=None,
                  ):
+        assert(isinstance(basis_size, int))
         if quadrature_size is None:
             quadrature_size = quadrature_points_per_direction()
+            assert(qp == quadrature_size[0] for qp in quadrature_size)
+            quadrature_size = quadrature_size[0]
         if slice_size is not None:
             quadrature_size = ceildiv(quadrature_size, slice_size)
-
         ImmutableRecord.__init__(self,
                                  quadrature_size=quadrature_size,
                                  basis_size=basis_size,
@@ -81,13 +83,14 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
             return self.basis_size
 
     def pymbolic(self, indices):
-        name = "{}{}Theta{}{}_qp{}_dof{}".format("face{}_".format(self.face) if self.face is not None else "",
-                                                 "d" if self.derivative else "",
-                                                 "T" if self.transpose else "",
-                                                 "_slice{}".format(self.slice_index) if self.slice_size is not None else "",
-                                                 self.quadrature_size,
-                                                 self.basis_size,
-                                                 )
+        name = "{}{}Theta{}{}_qp{}_dof{}" \
+               .format("face{}_".format(self.face) if self.face is not None else "",
+                       "d" if self.derivative else "",
+                       "T" if self.transpose else "",
+                       "_slice{}".format(self.slice_index) if self.slice_size is not None else "",
+                       self.quadrature_size,
+                       self.basis_size,
+                       )
         define_theta(name, self)
         return prim.Subscript(prim.Variable(name), indices)
 
@@ -199,51 +202,68 @@ class BasisTabulationMatrixArray(BasisTabulationMatrixBase):
 
 
 def quadrature_points_per_direction():
-    # Quadrature order
+    # Quadrature order per direction
     q = quadrature_order()
+    if isinstance(q, int):
+        q = (q,) * world_dimension()
 
     # Quadrature points in per direction
-    nb_qp = q // 2 + 1
+    nb_qp = tuple(order // 2 + 1 for order in q)
 
     return nb_qp
 
 
+def local_quadrature_points_per_direction():
+    """On a volume this simply gives the number of quadrature points per
+    direction. On a face it only returns the worlddim - 1 number of
+    quadrature points belonging to this face. (Delete normal direction).
+    """
+    qps_per_dir = quadrature_points_per_direction()
+    if local_dimension() != world_dimension():
+        facedir = get_global_context_value('facedir_s')
+        assert(get_global_context_value('facedir_s') == get_global_context_value('facedir_n') or
+               get_global_context_value('integral_type') == 'exterior_facet')
+        qps_per_dir = qps_per_dir[:facedir] + qps_per_dir[facedir + 1:]
+    return qps_per_dir
+
+
 def polynomial_degree():
     form = get_global_context_value("formdata").preprocessed_form
-    return form.coefficients()[0].ufl_element()._degree
+    degree = form.coefficients()[0].ufl_element().degree()
+    if isinstance(degree, int):
+        degree = (degree,) * world_dimension()
+    return degree
 
 
 def basis_functions_per_direction():
-    return polynomial_degree() + 1
+    return tuple(degree + 1 for degree in polynomial_degree())
 
 
-def define_oned_quadrature_weights(name):
+def define_oned_quadrature_weights(name, bound):
     loopy_class_member(name,
                        dtype=np.float64,
                        classtag="operator",
-                       shape=(quadrature_points_per_direction(),),
+                       shape=(bound,),
                        )
 
 
-def name_oned_quadrature_weights():
-    num = quadrature_points_per_direction()
-    name = "qw_num{}".format(num)
-    define_oned_quadrature_weights(name)
+def name_oned_quadrature_weights(bound):
+    name = "qw_num{}".format(bound)
+    define_oned_quadrature_weights(name, bound)
     return name
 
 
-def define_oned_quadrature_points(name):
+def define_oned_quadrature_points(name, bound):
     loopy_class_member(name,
                        dtype=np.float64,
                        classtag="operator",
-                       shape=(quadrature_points_per_direction(),),
+                       shape=(bound,),
                        )
 
 
-def name_oned_quadrature_points():
-    num = quadrature_points_per_direction()
-    name = "qp_num{}".format(num)
-    define_oned_quadrature_points(name)
+def name_oned_quadrature_points(bound):
+    name = "qp_num{}".format(bound)
+    define_oned_quadrature_points(name, bound)
     return name
 
 
@@ -284,12 +304,12 @@ def name_polynomials(degree):
 
 
 @preamble(kernel="operator")
-def sort_quadrature_points_weights(qp, qw):
+def sort_quadrature_points_weights(qp, qw, bound):
     range_field = lop_template_range_field()
     domain_field = name_domain_field()
-    number_qp = quadrature_points_per_direction()
     include_file("dune/perftool/sumfact/onedquadrature.hh", filetag="operatorfile")
-    return "onedQuadraturePointsWeights<{}, {}, {}>({}, {});".format(range_field, domain_field, number_qp, qp, qw)
+    return "onedQuadraturePointsWeights<{}, {}, {}>({}, {});" \
+        .format(range_field, domain_field, bound, qp, qw)
 
 
 @iname(kernel="operator")
@@ -325,9 +345,12 @@ def polynomial_lookup_mangler(target, func, dtypes):
 
 def define_theta(name, tabmat, additional_indices=(), width=None):
     assert isinstance(tabmat, BasisTabulationMatrix)
-    qp = name_oned_quadrature_points()
-    qw = name_oned_quadrature_weights()
-    sort_quadrature_points_weights(qp, qw)
+    bound = tabmat.quadrature_size
+    if tabmat.slice_size is not None:
+        bound *= tabmat.slice_size
+    qp = name_oned_quadrature_points(bound)
+    qw = name_oned_quadrature_weights(bound)
+    sort_quadrature_points_weights(qp, qw, bound)
     degree = tabmat.basis_size - 1
     polynomials = name_polynomials(degree)
 
@@ -370,17 +393,18 @@ def define_theta(name, tabmat, additional_indices=(), width=None):
 def construct_basis_matrix_sequence(transpose=False, derivative=None, facedir=None, facemod=None, basis_size=None):
     dim = world_dimension()
     result = [None] * dim
+    quadrature_size = quadrature_points_per_direction()
+    assert (basis_size is not None)
+    if facedir is not None:
+        quadrature_size = quadrature_size[:facedir] + (1,) + quadrature_size[facedir:]
 
     for i in range(dim):
-        quadrature_size = quadrature_points_per_direction()
-
         onface = None
         if facedir == i:
-            quadrature_size = 1
             onface = facemod
 
-        result[i] = BasisTabulationMatrix(quadrature_size=quadrature_size,
-                                          basis_size=basis_size,
+        result[i] = BasisTabulationMatrix(quadrature_size=quadrature_size[i],
+                                          basis_size=basis_size[i],
                                           transpose=transpose,
                                           derivative=derivative == i,
                                           face=onface)

python/dune/perftool/sumfact/vectorization.py

@@ -9,7 +9,6 @@ from dune.perftool.generation import (generator_factory,
 from dune.perftool.pdelab.restriction import (Restriction,
                                               restricted_name,
                                               )
-from dune.perftool.sumfact.tabulation import BasisTabulationMatrixArray
 from dune.perftool.error import PerftoolError
 from dune.perftool.options import get_option
 

python/dune/perftool/ufl/visitor.py

@@ -23,7 +23,9 @@ from pymbolic.primitives import (Call,
 from ufl.algorithms import MultiFunction
 from ufl.checks import is_cellwise_constant
 from ufl import (VectorElement,
+                 MixedElement,
                  TensorElement,
+                 TensorProductElement,
                  )
 from ufl.classes import (FixedIndex,
                          IndexSum,
@@ -98,7 +100,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         leaf_element = o.ufl_element()
 
         # Select the correct leaf element in the case of this being a mixed finite element
-        if o.ufl_element().num_sub_elements() > 0:
+        if isinstance(o.ufl_element(), MixedElement):
             index = self.indices[0]
             assert isinstance(index, int)
             self.indices = self.indices[1:]
@@ -127,7 +129,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             self.interface.initialize_function_spaces(o, self)
 
             index = None
-            if o.ufl_element().num_sub_elements() > 0:
+            if isinstance(o.ufl_element(), MixedElement):
                 index = self.indices[0]
                 assert isinstance(index, int)
                 self.indices = self.indices[1:]