unstructured grid working again

python/dune/codegen/blockstructured/geometry.py

-from dune.codegen.pdelab.restriction import restricted_name
-
 from dune.codegen.generation import (geometry_mixin,
                                      get_backend,
                                      temporary_variable,
                                      instruction,
-                                     get_global_context_value)
+                                     get_global_context_value,
+                                     domain)
 from dune.codegen.tools import get_pymbolic_basename
 from dune.codegen.options import (get_form_option,
                                   option_switch, get_option)
@@ -13,7 +12,11 @@ from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
                                           GenericPDELabGeometryMixin,
                                           world_dimension,
                                           local_dimension,
-                                          component_iname)
+                                          component_iname,
+                                          enforce_boundary_restriction,
+                                          restricted_name,
+                                          name_geometry,
+                                          )
 from dune.codegen.blockstructured.tools import sub_element_inames
 import pymbolic.primitives as prim
 from loopy.match import Writes
@@ -27,15 +30,49 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
     def facet_jacobian_determinant(self, o):
         raise NotImplementedError("I think this case was never really implemented")
 
+    def jacobian_determinant(self, o):
+        name = 'detjac'
+        self.define_jacobian_determinant(name)
+        return prim.Quotient(prim.Variable(name),
+                             prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+
+    def define_jacobian_determinant(self, name):
+        temporary_variable(name, shape=(), managed=True)
+
+        determinant_signed = name_jacobian_determinant_signed(self)
+
+        return instruction(expression=prim.Call(prim.Variable("abs"), (prim.Variable(determinant_signed),)),
+                           assignee=prim.Variable(name),
+                           within_inames=frozenset(sub_element_inames() + self.quadrature_inames()),
+                           depends_on=frozenset({Writes(determinant_signed)})
+                           )
+
     def jacobian_inverse(self, o):
-        raise NotImplementedError("I have not yet ported this one")
+        restriction = enforce_boundary_restriction(self)
+
+        assert(len(self.indices) == 2)
+        i, j = self.indices
+        self.indices = None
+
+        name = restricted_name("jit", restriction)
+        self.define_jacobian_inverse_transposed(name, restriction)
+
+        return prim.Product((prim.Subscript(prim.Variable(name), (j, i)),
+                             get_form_option("number_of_blocks")))
+
+    def define_jacobian_inverse_transposed(self, name, restriction):
+        temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
+
+        jacobian = name_jacobian_matrix(self)
+        det_inv = name_jacobian_determinant_inverse(self)
+
+        compute_inverse_transposed(name, det_inv, jacobian, self)
 
 
 @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()))
+        raise NotImplementedError("I think this case was never really implemented")
 
     def jacobian_determinant(self, o):
         return prim.Quotient(AxiparallelGeometryMixin.jacobian_determinant(self, o),
@@ -97,8 +134,8 @@ def bilinear_transformation_coefficients():
         raise NotImplementedError()
 
 
-def compute_jacobian(name):
-    pymbolic_pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+def compute_jacobian(name, visitor):
+    pymbolic_pos = visitor.quadrature_position_in_macro()
     jac_iname = component_iname(context="jac")
 
     coefficients = bilinear_transformation_coefficients()
@@ -137,23 +174,23 @@ def compute_jacobian(name):
     for i, expr in enumerate(expr_jac):
         instruction(expression=expr,
                     assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_iname), i)),
-                    within_inames=frozenset((jac_iname, ) + sub_element_inames() + get_backend(interface="quad_inames")()),
+                    within_inames=frozenset((jac_iname, ) + sub_element_inames() + visitor.quadrature_inames()),
                     depends_on=frozenset({Writes(get_pymbolic_basename(pymbolic_pos))} | {Writes(cd) for cd in coefficients})
                     )
 
 
-def define_jacobian_matrix(name):
+def define_jacobian_matrix(name, visitor):
     temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
-    compute_jacobian(name)
+    compute_jacobian(name, visitor)
 
 
-def name_jacobian_matrix():
+def name_jacobian_matrix(visitor):
     name = "jacobian"
-    define_jacobian_matrix(name)
+    define_jacobian_matrix(name, visitor)
     return name
 
 
-def compute_determinant(name, matrix):
+def compute_determinant(name, matrix, visitor):
     dim = world_dimension()
     matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
     if dim == 2:
@@ -171,52 +208,36 @@ def compute_determinant(name, matrix):
         raise NotImplementedError()
     instruction(expression=expr_determinant,
                 assignee=prim.Variable(name),
-                within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")()),
+                within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
                 depends_on=frozenset({Writes(matrix)})
                 )
 
 
-def define_jacobian_determinant(name):
+def define_jacobian_determinant(name, visitor):
     temporary_variable(name, shape=(), managed=True)
-    jacobian = name_jacobian_matrix()
-    compute_determinant(name, jacobian)
+    jacobian = name_jacobian_matrix(visitor)
+    compute_determinant(name, jacobian, visitor)
 
 
-def define_jacobian_determinant_inverse(name):
+def define_jacobian_determinant_inverse(name, visitor):
     temporary_variable(name, shape=(), managed=True)
-    determinant = name_jacobian_determinant()
+    determinant = name_jacobian_determinant_signed(visitor)
     return instruction(expression=prim.Quotient(1., prim.Variable(determinant)),
                        assignee=prim.Variable(name),
-                       within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")()),
-                       depends_on=frozenset({Writes(determinant)})
-                       )
-
-
-def define_jacobian_determinant_abs(name):
-    temporary_variable(name, shape=(), managed=True)
-    determinant = name_jacobian_determinant()
-    return instruction(expression=prim.Call(prim.Variable("abs"), (prim.Variable(determinant),)),
-                       assignee=prim.Variable(name),
-                       within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")()),
+                       within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
                        depends_on=frozenset({Writes(determinant)})
                        )
 
 
-def name_jacobian_determinant():
-    name = "detjac"
-    define_jacobian_determinant(name)
+def name_jacobian_determinant_signed(visitor):
+    name = "detjac_signed"
+    define_jacobian_determinant(name, visitor)
     return name
 
 
-def name_jacobian_determinant_inverse():
+def name_jacobian_determinant_inverse(visitor):
     name = "detjac_inverse"
-    define_jacobian_determinant_inverse(name)
-    return name
-
-
-def name_jacobian_determinant_abs():
-    name = "detjac_abs"
-    define_jacobian_determinant_abs(name)
+    define_jacobian_determinant_inverse(name, visitor)
     return name
 
 
@@ -231,7 +252,7 @@ def pymbolic_jacobian_determinant():
                          prim.Power(get_form_option("number_of_blocks"), local_dimension()))
 
 
-def compute_inverse_transposed(name, det_inv, matrix):
+def compute_inverse_transposed(name, det_inv, matrix, visitor):
     dim = world_dimension()
     matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
     assignee = [[prim.Subscript(prim.Variable(name), (i, j)) for j in range(dim)] for i in range(dim)]
@@ -278,20 +299,20 @@ def compute_inverse_transposed(name, det_inv, matrix):
         for i in range(dim):
             instruction(expression=exprs[i][j],
                         assignee=assignee[i][j],
-                        within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")()),
+                        within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
                         depends_on=frozenset({Writes(matrix), Writes(det_inv)}))
 
 
-def define_jacobian_inverse_transposed(name):
+def define_jacobian_inverse_transposed(name, visitor):
     temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
-    jacobian = name_jacobian_matrix()
-    det_inv = name_jacobian_determinant_inverse()
-    compute_inverse_transposed(name, det_inv, jacobian)
+    jacobian = name_jacobian_matrix(visitor)
+    det_inv = name_jacobian_determinant_inverse(visitor)
+    compute_inverse_transposed(name, det_inv, jacobian, visitor)
 
 
-def name_jacobian_inverse_transposed(restriction):
+def name_jacobian_inverse_transposed(restriction, visitor):
     name = restricted_name("jit", restriction)
-    define_jacobian_inverse_transposed(name)
+    define_jacobian_inverse_transposed(name, visitor)
     return name
 
 #
@@ -408,8 +429,8 @@ def apply_constant_to_global_transformation(name, local):
                 )
 
 
-def to_global(local):
-    macro = name_point_in_macro(local)
+def to_global(local, visitor):
+    macro = name_point_in_macro(local, visitor)
     name = macro + "_global"
 
     it = get_global_context_value("integral_type")

test/blockstructured/poisson/poisson_unstructured.mini

@@ -19,9 +19,10 @@ compare_l2errorsquared = 1e-7
 grid_unstructured = 1
 
 [formcompiler.r]
-matrix_free = 1
+matrix_free = 0
 blockstructured = 1
 number_of_blocks = 16, 8 | expand dimension
+geometry_mixins = blockstructured_multilinear
 
 [formcompiler.ufl_variants]
 cell = quadrilateral, hexahedron | expand dimension