[!345] Reimplement selection of correct finite element implementation

Merge branch 'feature/finite-element-implementation-selection' into 'master'

ref:extensions/dune-codegen So far, we were a bit sloppy about which basis we
choose from the UFL counterpart. This commit implements the necessary 1-to-1
mapping.

This fixes [#146]

See merge request [extensions/dune-codegen!345]

  [#146]: gitlab.dune-project.org/NoneNone/issues/146
  [extensions/dune-codegen!345]: gitlab.dune-project.org/extensions/dune-codegen/merge_requests/345


Closes #146

python/dune/codegen/error.py

@@ -23,3 +23,7 @@ class CodegenVectorizationError(CodegenCodegenError):
 
 class CodegenAutotuneError(CodegenVectorizationError):
     pass
+
+
+class CodegenUnsupportedFiniteElementError(CodegenUFLError):
+    pass

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

@@ -22,6 +22,8 @@ from dune.codegen.generation import (generator_factory,
 from dune.codegen.options import (get_form_option,
                                   get_option,
                                   )
+from ufl import TensorProductCell
+
 
 #
 # The following functions are not doing anything useful, but providing easy access
@@ -84,6 +86,9 @@ def isLagrange(fem):
 
 
 def isSimplical(cell):
+    if isinstance(cell, TensorProductCell):
+        return False
+
     # Cells can be identified through strings *or* ufl objects
     if not isinstance(cell, str):
         cell = cell.cellname()
@@ -91,6 +96,9 @@ def isSimplical(cell):
 
 
 def isQuadrilateral(cell):
+    if isinstance(cell, TensorProductCell):
+        return all(tuple(isSimplical(c) for c in cell.sub_cells()))
+
     # Cells can be identified through strings *or* ufl objects
     if not isinstance(cell, str):
         cell = cell.cellname()
@@ -123,20 +131,10 @@ def FEM_name_mangling(fem):
             name = name + FEM_name_mangling(elem)
         return name
     if isinstance(fem, FiniteElement):
-        if isPk(fem):
-            return "P" + str(fem._degree)
-        if isQk(fem):
-            return "Q" + str(fem._degree)
-        if isDG(fem):
-            return "DG" + str(fem._degree)
+        return "{}{}".format(fem._short_name, 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")
+        return "TP_{}".format("_".join(FEM_name_mangling(subel) for subel in fem.sub_elements()))
 
     raise NotImplementedError("FEM NAME MANGLING")
 

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

+from dune.codegen.error import CodegenUnsupportedFiniteElementError
 from dune.codegen.generation import (include_file,
                                      preamble,
                                      )
@@ -9,10 +10,6 @@ from dune.codegen.pdelab.driver import (FEM_name_mangling,
                                         get_dimension,
                                         get_test_element,
                                         get_trial_element,
-                                        isLagrange,
-                                        isDG,
-                                        isPk,
-                                        isQk,
                                         isQuadrilateral,
                                         isSimplical,
                                         name_initree,
@@ -20,7 +17,7 @@ from dune.codegen.pdelab.driver import (FEM_name_mangling,
                                         )
 from dune.codegen.loopy.target import type_floatingpoint
 
-from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement, TensorProductCell
 
 
 @preamble(section="grid")
@@ -123,61 +120,98 @@ def name_leafview():
     return name
 
 
+def get_short_name(element):
+    if isinstance(element, TensorProductElement):
+        assert len(set(subel._short_name for subel in element.sub_elements())) == 1
+        return get_short_name(element.sub_elements()[0])
+
+    return element._short_name
+
+
 @preamble(section="fem")
 def typedef_fem(element, name):
     gv = type_leafview()
     df = type_domainfield()
     r = type_range()
     dim = get_dimension()
+    cell = element.cell()
+    degree = element.degree()
+    short = get_short_name(element)
+
+    # We currently only support TensorProductElement from UFL if it aliases another finite element
+    # available from UFL. Here, we check this condition and recover the aliases element
+    if isinstance(element, TensorProductElement):
+        subels = set(subel._short_name for subel in element.sub_elements())
+        if len(subels) != 1 or len(set(element.degree())) != 1:
+            raise CodegenUnsupportedFiniteElementError(element)
+
+        degree = element.degree()[0]
+        cell = TensorProductCell(*tuple(subel.cell() for subel in element.sub_elements()))
 
+    # The blockstructured code branch has its own handling of finite element selection
     if get_form_option("blockstructured"):
         include_file("dune/codegen/blockstructured/blockstructuredqkfem.hh", filetag="driver")
-        degree = element.degree() * get_form_option("number_of_blocks")
+        degree = degree * get_form_option("number_of_blocks")
         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-codegen")
-    elif isQk(element):
-        include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
-        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())
-    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)
-        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)
-        raise NotImplementedError("Geometry type not known in code generation")
-
-    raise NotImplementedError("FEM not implemented in dune-codegen")
+    # This is a backward-compatibility hack: So far we silently used OPBFem for DG with simplices:
+    if short == "DG" and isSimplical(cell):
+        short = "OPB"
+
+    # Choose the correct finite element implementation
+    if short == "CG":
+        if isSimplical(cell):
+            if dim in (1, 2, 3):
+                include_file("dune/pdelab/finiteelementmap/pkfem.hh", filetag="driver")
+                return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+                    .format(name, gv, df, r, degree)
+            else:
+                raise CodegenUnsupportedFiniteElementError(element)
+        elif isQuadrilateral(cell):
+            if dim in (2, 3) and degree < 3:
+                include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
+                return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+                    .format(name, gv, df, r, degree)
+            else:
+                raise CodegenUnsupportedFiniteElementError(element)
+        else:
+            raise CodegenUnsupportedFiniteElementError(element)
+    elif short == "DG":
+        if isQuadrilateral(cell):
+            if dim < 4:
+                include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
+                return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;" \
+                    .format(name, df, r, degree, dim)
+            else:
+                raise CodegenUnsupportedFiniteElementError(element)
+        else:
+            raise CodegenUnsupportedFiniteElementError(element)
+    elif short == "GL":
+        raise NotImplementedError("Gauss-Legendre polynomials not implemented")
+    elif short == "DGLL":
+        raise NotImplementedError("Discontinuous Gauss-Lobatto-Legendre polynomials not implemented")
+    elif short == "OPB":
+        if isQuadrilateral(cell):
+            gt = "Dune::GeometryType::cube"
+        elif isSimplical(cell):
+            gt = "Dune::GeometryType::simplex"
+        else:
+            raise CodegenUnsupportedFiniteElementError(element)
+
+        include_file("dune/pdelab/finiteelementmap/opbfem.hh", filetag="driver")
+        return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, {}>;" \
+            .format(name, df, r, degree, dim, gt)
+    elif short == "Monom":
+        raise NotImplementedError("Monomials basis DG not implemented")
+    elif short == "RaTu":
+        raise NotImplementedError("Rannacher-Turek elements not implemented")
+    elif short == "RT":
+        raise NotImplementedError("Raviart-Thomas elements not implemented")
+    elif short == "BDM":
+        raise NotImplementedError("Brezzi-Douglas-Marini elements not implemented")
+    else:
+        raise CodegenUnsupportedFiniteElementError(element)
 
 
 def type_fem(element):
@@ -189,23 +223,13 @@ def type_fem(element):
 @preamble(section="fem")
 def define_fem(element, name):
     femtype = type_fem(element)
-    from dune.codegen.pdelab.driver import isDG
-    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))
+
+    # Determine whether the FEM is grid-dependent - currently on the Lagrange elements are
+    if get_short_name(element) == "CG":
         gv = name_leafview()
         return "{} {}({});".format(femtype, name, gv)
+    else:
+        return "{} {};".format(femtype, name)
 
 
 def name_fem(element):