TensorProductElements (with same degree in all directions)

- Use TensorProductElement in one example.

- We still use the sumfact option to swicth to the sum factorization
  code branch. This way it is still possible to easily swicth between
  sumfact and non sumfact code.

- Only TensorProductElements with the same degree in all directions
  will work. Anisotropie and adaption of quadrature rule will happen
  in the next commits.

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,47 @@ 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)
+
+        # 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, max(element.degree()))
+        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, max(element.degree()), 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 +181,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 +227,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 +238,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 +265,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 +275,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/quadrature.py

@@ -191,7 +191,10 @@ def _estimate_quadrature_order():
     integrals = form.integrals_by_type(integral_type)
     polynomial_degree = 0
     for i in integrals:
-        polynomial_degree = max(polynomial_degree, i.metadata()['estimated_polynomial_degree'])
+        degree = i.metadata()['estimated_polynomial_degree']
+        if isinstance(degree, tuple):
+            degree = max(degree)
+        polynomial_degree = max(polynomial_degree, degree)
 
     return polynomial_degree
 

python/dune/perftool/sumfact/accumulation.py

@@ -41,7 +41,7 @@ import loopy as lp
 import numpy as np
 import pymbolic.primitives as prim
 import ufl.classes as uc
-
+from ufl import FiniteElement, TensorProductElement
 
 @iname
 def _sumfact_iname(bound, _type, count):
@@ -138,8 +138,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()):
@@ -200,9 +199,12 @@ def generate_accumulation_instruction(expr, visitor):
     trial_info = visitor.trial_info
 
     leaf_element = test_info.element
-    if leaf_element.num_sub_elements() > 0:
+    if leaf_element.num_sub_elements() > 0 and not isinstance(leaf_element, TensorProductElement):
         leaf_element = leaf_element.extract_component(test_info.element_index)[1]
-    basis_size = leaf_element.degree() + 1
+    degree = leaf_element._degree
+    if isinstance(degree, tuple):
+        degree = max(degree)
+    basis_size = degree + 1
 
     from dune.perftool.pdelab.localoperator import boundary_predicates
     predicates = boundary_predicates(expr,

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 VectorElement, TensorElement, TensorProductElement
 
 from pytools import product, ImmutableRecord
 
@@ -96,20 +96,23 @@ 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, tuple):
+        degree = max(degree)
+    return degree + 1
 
 
 @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 element.num_sub_elements() > 0 and not isinstance(element, TensorProductElement):
         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)
@@ -190,8 +193,8 @@ 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:

python/dune/perftool/ufl/visitor.py

@@ -24,6 +24,7 @@ from ufl.algorithms import MultiFunction
 from ufl.checks import is_cellwise_constant
 from ufl import (VectorElement,
                  TensorElement,
+                 TensorProductElement,
                  )
 from ufl.classes import (FixedIndex,
                          IndexSum,
@@ -98,7 +99,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 o.ufl_element().num_sub_elements() > 0 and not isinstance(o.ufl_element(),TensorProductElement):
             index = self.indices[0]
             assert isinstance(index, int)
             self.indices = self.indices[1:]
@@ -127,7 +128,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             self.interface.initialize_function_spaces(o, self)
 
             index = None
-            if o.ufl_element().num_sub_elements() > 0:
+            if o.ufl_element().num_sub_elements() > 0 and not isinstance(o.ufl_element(), TensorProductElement):
                 index = self.indices[0]
                 assert isinstance(index, int)
                 self.indices = self.indices[1:]

test/sumfact/poisson/poisson_2d.ufl

@@ -5,7 +5,10 @@ c = (0.5-x[0])**2 + (0.5-x[1])**2
 g = exp(-1.*c)
 f = 2*(2.-2*c)*g
 
-V = FiniteElement("CG", cell, degree)
+V_0 = FiniteElement("CG", interval, degree)
+V_1 = FiniteElement("CG", interval, degree)
+V = TensorProductElement(V_0, V_1, cell=cell)
+
 u = TrialFunction(V)
 v = TestFunction(V)