adds to_global

python/dune/codegen/blockstructured/geometry.py

@@ -27,6 +27,20 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
     def spatial_coordinate(self, o):
         return self.to_global(self.quadrature_position_in_macro())
 
+    def to_global(self, local):
+        assert isinstance(local, prim.Expression)
+        name = get_pymbolic_basename(local) + "_global"
+
+        it = get_global_context_value("integral_type")
+
+        if it == 'cell':
+            compute_multilinear_to_global_transformation(name, local, self)
+            return prim.Variable(name)
+        elif it == 'exterior_facet' or it == 'interior_facet':
+            return GenericPDELabGeometryMixin.to_global(self, local)
+        else:
+            raise NotImplementedError
+
     def facet_jacobian_determinant(self, o):
         raise NotImplementedError("I think this case was never really implemented")
 
@@ -71,6 +85,20 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
 
 @geometry_mixin("blockstructured_axiparallel")
 class AxiparallelBlockStructuredGeometryMixin(AxiparallelGeometryMixin, BlockStructuredGeometryMixin):
+    def to_global(self, local):
+        assert isinstance(local, prim.Expression)
+        name = get_pymbolic_basename(local) + "_global"
+
+        it = get_global_context_value("integral_type")
+
+        if it == 'cell':
+            compute_axiparallel_to_global_transformation(name, local, self)
+            return prim.Variable(name)
+        elif it == 'exterior_facet' or it == 'interior_facet':
+            return AxiparallelGeometryMixin.to_global(self, local)
+        else:
+            raise NotImplementedError
+
     def facet_jacobian_determinant(self, o):
         raise NotImplementedError("I think this case was never really implemented")
 
@@ -241,17 +269,6 @@ def name_jacobian_determinant_inverse(visitor):
     return name
 
 
-# scale determinant according to the order of the blockstructuring
-def pymbolic_jacobian_determinant():
-    if get_form_option("constant_transformation_matrix"):
-        from dune.codegen.pdelab.geometry import name_jacobian_determinant
-        n_det = name_jacobian_determinant()
-    else:
-        n_det = name_jacobian_determinant_abs()
-    return prim.Quotient(prim.Variable(n_det),
-                         prim.Power(get_form_option("number_of_blocks"), local_dimension()))
-
-
 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)]
@@ -369,7 +386,7 @@ def name_point_in_macro(point_in_micro, visitor):
     return name
 
 
-def apply_default_to_global_transformation(name, local):
+def compute_multilinear_to_global_transformation(name, local, visitor):
     dim = world_dimension()
     temporary_variable(name, shape=(dim,), managed=True)
     corners = name_element_corners()
@@ -400,13 +417,13 @@ def apply_default_to_global_transformation(name, local):
     assignee = prim.Subscript(prim.Variable(name), (dim_pym,))
 
     instruction(assignee=assignee, expression=expr,
-                within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")() + (dim_pym.name,)),
+                within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames() + (dim_pym.name,)),
                 within_inames_is_final=True,
-                depends_on=frozenset({Writes(local.name), Writes(corners)})
+                depends_on=frozenset({Writes(get_pymbolic_basename(local)), Writes(corners)})
                 )
 
 
-def apply_constant_to_global_transformation(name, local):
+def compute_axiparallel_to_global_transformation(name, local, visitor):
     dim = world_dimension()
     temporary_variable(name, shape=(dim,), managed=True)
     corners = name_element_corners()
@@ -423,32 +440,12 @@ def apply_constant_to_global_transformation(name, local):
     assignee = prim.Subscript(prim.Variable(name), (dim_pym,))
 
     instruction(assignee=assignee, expression=expr,
-                within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")() + (dim_pym.name,)),
+                within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames() + (dim_pym.name,)),
                 within_inames_is_final=True,
-                depends_on=frozenset({Writes(local.name), Writes(corners)})
+                depends_on=frozenset({Writes(get_pymbolic_basename(local)), Writes(corners)})
                 )
 
 
-def to_global(local, visitor):
-    macro = name_point_in_macro(local, visitor)
-    name = macro + "_global"
-
-    it = get_global_context_value("integral_type")
-
-    if it == 'cell':
-        if get_form_option("constant_transformation_matrix"):
-            apply_constant_to_global_transformation(name, prim.Variable(macro))
-        else:
-            apply_default_to_global_transformation(name, prim.Variable(macro))
-    elif it == 'exterior_facet' or it == 'interior_facet':
-        from dune.codegen.pdelab.geometry import apply_to_global_transformation
-        apply_to_global_transformation(name, prim.Variable(macro))
-    else:
-        raise NotImplementedError
-
-    return prim.Variable(name)
-
-
 def define_element_corners(name):
     from dune.codegen.pdelab.driver import get_form
     n_corners = get_form().ufl_cell().num_vertices()