My first attempt to apply mixins to blockstructured stuff

python/dune/codegen/blockstructured/__init__.py

+import dune.codegen.blockstructured.accumulation
 import dune.codegen.blockstructured.quadrature
 import dune.codegen.blockstructured.argument
 import dune.codegen.blockstructured.geometry

python/dune/codegen/blockstructured/basis.py

 from dune.codegen.generation import (backend,
+                                     basis_mixin,
                                      kernel_cached,
                                      get_backend,
                                      instruction,
@@ -9,7 +10,8 @@ from dune.codegen.generation import (backend,
                                      include_file,)
 from dune.codegen.tools import get_pymbolic_basename
 from dune.codegen.loopy.target import type_floatingpoint
-from dune.codegen.pdelab.basis import (declare_cache_temporary,
+from dune.codegen.pdelab.basis import (GenericBasisMixin,
+                                       declare_cache_temporary,
                                        name_localbasis_cache,
                                        type_localbasis,
                                        FEM_name_mangling)
@@ -27,6 +29,56 @@ from ufl import MixedElement
 import pymbolic.primitives as prim
 
 
+@basis_mixin("blockstructured")
+class BlockStructuredBasisMixin(GenericBasisMixin):
+    def lfs_inames(self, element, restriction, number, context=""):
+        return lfs_inames(element, restriction, number, context=context)
+
+    def implement_basis(self, element, restriction, number, context=''):
+        assert not isinstance(element, MixedElement)
+        name = "phi_{}".format(FEM_name_mangling(element))
+        name = restricted_name(name, restriction)
+        self.evaluate_basis(element, name, restriction)
+        inames = self.lfs_inames(element, restriction, number, context=context)
+
+        return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, element.degree() + 1), 0))
+
+    @kernel_cached
+    def evaluate_basis(self, element, name, restriction):
+        temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1),
+                           decl_method=declare_cache_temporary(element, restriction, 'Function'))
+        cache = name_localbasis_cache(element)
+        qp = self.to_cell(self.quadrature_position())
+        localbasis = name_localbasis(element)
+        instruction(inames=self.quadrature_inames(),
+                    code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
+                    assignees=frozenset({name}),
+                    read_variables=frozenset({get_pymbolic_basename(qp)}),
+                    )
+
+    def implement_reference_gradient(self, element, restriction, number, context=''):
+        assert not isinstance(element, MixedElement)
+        name = "js_{}".format(FEM_name_mangling(element))
+        name = restricted_name(name, restriction)
+        self.evaluate_reference_gradient(element, name, restriction)
+        inames = self.lfs_inames(element, restriction, number, context=context)
+
+        return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, element.degree() + 1), 0))
+
+    @kernel_cached
+    def evaluate_reference_gradient(self, element, name, restriction):
+        temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1, world_dimension()),
+                           decl_method=declare_cache_temporary(element, restriction, 'Jacobian'))
+        cache = name_localbasis_cache(element)
+        qp = self.to_cell(self.quadrature_position())
+        localbasis = name_localbasis(element)
+        instruction(inames=self.quadrature_inames(),
+                    code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
+                    assignees=frozenset({name}),
+                    read_variables=frozenset({get_pymbolic_basename(qp)}),
+                    )
+
+
 # define FE basis explicitly in localoperator
 @backend(interface="typedef_localbasis", name="blockstructured")
 @class_member(classtag="operator")
@@ -68,49 +120,50 @@ def name_localbasis(leaf_element):
     return name
 
 
-@kernel_cached
-def evaluate_basis(leaf_element, name, restriction):
-    temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1),
-                       decl_method=declare_cache_temporary(leaf_element, restriction, 'Function'))
-    cache = name_localbasis_cache(leaf_element)
-    qp = pymbolic_quadrature_position_in_cell(restriction)
-    localbasis = name_localbasis(leaf_element)
-    instruction(inames=get_backend("quad_inames")(),
-                code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
-                assignees=frozenset({name}),
-                read_variables=frozenset({get_pymbolic_basename(qp)}),
-                )
-
-
-def pymbolic_basis(leaf_element, restriction, number, context=''):
-    assert not isinstance(leaf_element, MixedElement)
-    name = "phi_{}".format(FEM_name_mangling(leaf_element))
-    name = restricted_name(name, restriction)
-    evaluate_basis(leaf_element, name, restriction)
-    inames = lfs_inames(leaf_element, restriction, number, context=context)
-
-    return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))
-
-
-@kernel_cached
-def evaluate_reference_gradient(leaf_element, name, restriction):
-    temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1, world_dimension()),
-                       decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
-    cache = name_localbasis_cache(leaf_element)
-    qp = pymbolic_quadrature_position_in_cell(restriction)
-    localbasis = name_localbasis(leaf_element)
-    instruction(inames=get_backend("quad_inames")(),
-                code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
-                assignees=frozenset({name}),
-                read_variables=frozenset({get_pymbolic_basename(qp)}),
-                )
-
-
-def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
-    assert not isinstance(leaf_element, MixedElement)
-    name = "js_{}".format(FEM_name_mangling(leaf_element))
-    name = restricted_name(name, restriction)
-    evaluate_reference_gradient(leaf_element, name, restriction)
-    inames = lfs_inames(leaf_element, restriction, number, context=context)
-
-    return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))
+# @kernel_cached
+# def evaluate_basis(leaf_element, name, restriction):
+#     temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1),
+#                        decl_method=declare_cache_temporary(leaf_element, restriction, 'Function'))
+#     cache = name_localbasis_cache(leaf_element)
+#     qp = pymbolic_quadrature_position_in_cell(restriction)
+#     localbasis = name_localbasis(leaf_element)
+#     instruction(inames=get_backend("quad_inames")(),
+#                 code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
+#                 assignees=frozenset({name}),
+#                 read_variables=frozenset({get_pymbolic_basename(qp)}),
+#                 )
+#
+#
+# def pymbolic_basis(leaf_element, restriction, number, context=''):
+#     assert not isinstance(leaf_element, MixedElement)
+#     name = "phi_{}".format(FEM_name_mangling(leaf_element))
+#     name = restricted_name(name, restriction)
+#     evaluate_basis(leaf_element, name, restriction)
+#     inames = lfs_inames(leaf_element, restriction, number, context=context)
+#
+#     return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))
+
+#
+# @kernel_cached
+# def evaluate_reference_gradient(leaf_element, name, restriction):
+#     from pudb import set_trace; set_trace()
+#     temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1, world_dimension()),
+#                        decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
+#     cache = name_localbasis_cache(leaf_element)
+#     qp = pymbolic_quadrature_position_in_cell(restriction)
+#     localbasis = name_localbasis(leaf_element)
+#     instruction(inames=get_backend("quad_inames")(),
+#                 code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
+#                 assignees=frozenset({name}),
+#                 read_variables=frozenset({get_pymbolic_basename(qp)}),
+#                 )
+#
+#
+# def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
+#     assert not isinstance(leaf_element, MixedElement)
+#     name = "js_{}".format(FEM_name_mangling(leaf_element))
+#     name = restricted_name(name, restriction)
+#     evaluate_reference_gradient(leaf_element, name, restriction)
+#     inames = lfs_inames(leaf_element, restriction, number, context=context)
+#
+#     return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))

python/dune/codegen/blockstructured/geometry.py

 from dune.codegen.pdelab.restriction import restricted_name
 
-from dune.codegen.generation import (get_backend,
+from dune.codegen.generation import (geometry_mixin,
+                                     get_backend,
                                      temporary_variable,
                                      instruction,
                                      get_global_context_value)
 from dune.codegen.tools import get_pymbolic_basename
 from dune.codegen.options import (get_form_option,
                                   option_switch, get_option)
-from dune.codegen.pdelab.geometry import (world_dimension,
+from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
+                                          EquidistantGeometryMixin,
+                                          GenericPDELabGeometryMixin,
+                                          world_dimension,
                                           local_dimension,
                                           component_iname)
 from dune.codegen.blockstructured.tools import sub_element_inames
@@ -15,6 +19,31 @@ import pymbolic.primitives as prim
 from loopy.match import Writes
 
 
+@geometry_mixin("blockstructured_multilinear")
+class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
+    def facet_jacobian_determinant(self, o):
+        raise NotImplementedError("I think this case was never really implemented")
+
+    def jacobian_inverse(self, o):
+        raise NotImplementedError("I have not yet ported this one")
+
+
+@geometry_mixin("blockstructured_axiparallel")
+class AxiparallelBlockStructuredGeometryMixin(AxiparallelGeometryMixin, BlockStructuredGeometryMixin):
+    def facet_jacobian_determinant(self, o):
+        return prim.Quotient(AxiparallelGeometryMixin.facet_jacobian_determinant(self, o),
+                             prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+
+    def jacobian_inverse(self, o):
+        return prim.Quotient(AxiparallelGeometryMixin.facet_jacobian_determinant(self, o),
+                             prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+
+
+@geometry_mixin("blockstructured_equidistant")
+class EquidistantBlockStructuredGeometryMixin(EquidistantGeometryMixin, AxiparallelBlockStructuredGeometryMixin):
+    pass
+
+
 # TODO warum muss within_inames_is_final=True gesetzt werden?
 def compute_corner_diff(first, second, additional_terms=tuple()):
     corners = name_element_corners()
@@ -258,23 +287,23 @@ def name_jacobian_inverse_transposed(restriction):
     define_jacobian_inverse_transposed(name)
     return name
 
-
-# scale Jacobian according to the order of the blockstructure
-def pymbolic_jacobian_inverse(i, j, restriction):
-    if get_form_option("constant_transformation_matrix"):
-        from dune.codegen.pdelab.geometry import name_constant_jacobian_inverse_transposed
-        name_jit = name_constant_jacobian_inverse_transposed(restriction)
-    else:
-        name_jit = name_jacobian_inverse_transposed(restriction)
-    return prim.Product((get_form_option("number_of_blocks"),
-                         prim.Subscript(prim.Variable(name_jit), (j, i))))
-
-
-# scale determinant according to the order of the blockstructure
-def pymbolic_facet_jacobian_determinant():
-    from dune.codegen.pdelab.geometry import name_facet_jacobian_determinant
-    return prim.Quotient(prim.Variable(name_facet_jacobian_determinant()),
-                         prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+#
+# # scale Jacobian according to the order of the blockstructure
+# def pymbolic_jacobian_inverse(i, j, restriction):
+#     if get_form_option("constant_transformation_matrix"):
+#         from dune.codegen.pdelab.geometry import name_constant_jacobian_inverse_transposed
+#         name_jit = name_constant_jacobian_inverse_transposed(restriction)
+#     else:
+#         name_jit = name_jacobian_inverse_transposed(restriction)
+#     return prim.Product((get_form_option("number_of_blocks"),
+#                          prim.Subscript(prim.Variable(name_jit), (j, i))))
+#
+#
+# # scale determinant according to the order of the blockstructure
+# def pymbolic_facet_jacobian_determinant():
+#     from dune.codegen.pdelab.geometry import name_facet_jacobian_determinant
+#     return prim.Quotient(prim.Variable(name_facet_jacobian_determinant()),
+#                          prim.Power(get_form_option("number_of_blocks"), local_dimension()))
 
 
 # translate a point in the micro element into macro coordinates

python/dune/codegen/blockstructured/quadrature.py

 from dune.codegen.generation import (backend,
                                      quadrature_mixin,
                                      )
-from dune.codegen.pdelab.quadrature import (GenericQuadratureMixin,
-                                            quadrature_iname)
+from dune.codegen.pdelab.quadrature import GenericQuadratureMixin
 from dune.codegen.blockstructured.geometry import name_point_in_macro
 import pymbolic.primitives as prim
 
 
 @quadrature_mixin("blockstructured")
 class BlockstructuredQuadratureMixin(GenericQuadratureMixin):
-    pass
-
-
-@backend(interface="quad_pos", name='blockstructured')
-def pymbolic_quadrature_position():
-    quad_points = name_quadrature_points()
-    quad_iname = quadrature_iname()
-
-    qp_in_micro = prim.Subscript(prim.Variable(quad_points), (prim.Variable(quad_iname),))
-
-    return prim.Variable(name_point_in_macro(qp_in_micro))
-
-
-@backend(interface="qp_in_cell", name="blockstructured")
-def pymbolic_quadrature_position_in_cell(restriction):
-    from dune.codegen.pdelab.geometry import to_cell_coordinates
-    quad_pos = pymbolic_quadrature_position()
-    return to_cell_coordinates(quad_pos, restriction)
+    def quadrature_position(self):
+        original = GenericQuadratureMixin.quadrature_position(self)
+        return prim.Variable(name_point_in_macro(original))
+
+#
+# @backend(interface="quad_pos", name='blockstructured')
+# def pymbolic_quadrature_position():
+#     quad_points = name_quadrature_points()
+#     quad_iname = quadrature_iname()
+#
+#     qp_in_micro = prim.Subscript(prim.Variable(quad_points), (prim.Variable(quad_iname),))
+#
+#     return prim.Variable(name_point_in_macro(qp_in_micro))
+#
+#
+# @backend(interface="qp_in_cell", name="blockstructured")
+# def pymbolic_quadrature_position_in_cell(restriction):
+#     from dune.codegen.pdelab.geometry import to_cell_coordinates
+#     quad_pos = pymbolic_quadrature_position()
+#     return to_cell_coordinates(quad_pos, restriction)

python/dune/codegen/options.py

@@ -180,6 +180,12 @@ def process_form_options(opt, form):
                        accumulation_mixins="sumfact",
                        )
 
+    if opt.blockstructured:
+        opt = opt.copy(accumulation_mixins="blockstructured",
+                       geometry_mixins="blockstructured",
+                       quadrature_mixins="blockstructured",
+                       )
+
     if opt.control:
         opt = opt.copy(accumulation_mixins="control")
 

test/CMakeLists.txt

@@ -7,5 +7,5 @@ add_subdirectory(stokes)
 
 add_subdirectory(sumfact)
 add_subdirectory(coeffeval)
-#add_subdirectory(blockstructured)
+add_subdirectory(blockstructured)
 add_subdirectory(adjoint)