Fully implement axiparallel facet normals for sum factorization

python/dune/perftool/pdelab/__init__.py

@@ -19,8 +19,8 @@ from dune.perftool.pdelab.geometry import (component_iname,
                                            pymbolic_facet_jacobian_determinant,
                                            pymbolic_jacobian_determinant,
                                            pymbolic_jacobian_inverse_transposed,
-                                           name_unit_inner_normal,
-                                           name_unit_outer_normal,
+                                           pymbolic_unit_inner_normal,
+                                           pymbolic_unit_outer_normal,
                                            to_global,
                                            )
 from dune.perftool.pdelab.index import (name_index,
@@ -137,11 +137,11 @@ class PDELabInterface(object):
     def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
         return pymbolic_jacobian_inverse_transposed(i, j, restriction)
 
-    def name_unit_inner_normal(self):
-        return name_unit_inner_normal()
+    def pymbolic_unit_inner_normal(self):
+        return pymbolic_unit_inner_normal()
 
-    def name_unit_outer_normal(self):
-        return name_unit_outer_normal()
+    def pymbolic_unit_outer_normal(self):
+        return pymbolic_unit_outer_normal()
 
     def name_cell_volume(self, restriction):
         return name_cell_volume(restriction)

python/dune/perftool/pdelab/geometry.py

@@ -262,7 +262,7 @@ def declare_normal(name, shape, shape_impl):
     return "auto {} = {}.centerUnitOuterNormal();".format(name, ig)
 
 
-def name_unit_outer_normal():
+def pymbolic_unit_outer_normal():
     name = "outer_normal"
     if not get_option("diagonal_transformation_matrix"):
         temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
@@ -270,7 +270,7 @@ def name_unit_outer_normal():
     else:
         declare_normal(name, None, None)
         globalarg(name, shape=(world_dimension(),), dtype=np.float64)
-    return "outer_normal"
+    return prim.Variable(name)
 
 
 def evaluate_unit_inner_normal(name):
@@ -281,11 +281,11 @@ def evaluate_unit_inner_normal(name):
                                )
 
 
-def name_unit_inner_normal():
+def pymbolic_unit_inner_normal():
     name = "inner_normal"
     temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
     evaluate_unit_inner_normal(name)
-    return "inner_normal"
+    return prim.Variable(name)
 
 
 def type_jacobian_inverse_transposed(restriction):

python/dune/perftool/sumfact/__init__.py

@@ -76,3 +76,15 @@ class SumFactInterface(PDELabInterface):
         ret, indices = get_backend(interface="spatial_coordinate", selector=option_switch("diagonal_transformation_matrix"))(self.visitor.indices, self.visitor.do_predicates, self.visitor)
         self.visitor.indices = indices
         return ret
+
+    def pymbolic_unit_outer_normal(self):
+        from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
+        ret, indices = pymbolic_unit_outer_normal(self.visitor.indices)
+        self.visitor.indices = indices
+        return ret
+
+    def pymbolic_unit_inner_normal(self):
+        from dune.perftool.sumfact.geometry import pymbolic_unit_inner_normal
+        ret, indices = pymbolic_unit_inner_normal(self.visitor.indices)
+        self.visitor.indices = indices
+        return ret

python/dune/perftool/sumfact/geometry.py

@@ -18,7 +18,8 @@ from dune.perftool.pdelab.geometry import (local_dimension,
                                            )
 from dune.perftool.sumfact.symbolic import SumfactKernelInputBase
 from dune.perftool.sumfact.vectorization import attach_vectorization_info
-from dune.perftool.options import option_switch
+from dune.perftool.options import get_option, option_switch
+from dune.perftool.ufl.modified_terminals import Restriction
 
 from pytools import ImmutableRecord
 
@@ -159,3 +160,37 @@ def pymbolic_spatial_coordinate_axiparallel(visitor_indices, do_predicates, visi
         x = pymbolic_quadrature_position(iindex, visitor)
 
     return prim.Subscript(prim.Variable(lowcorner), (index,)) + x * prim.Subscript(prim.Variable(meshwidth), (index,)), None
+
+
+def pymbolic_unit_outer_normal(visitor_indices):
+    index, = visitor_indices
+    assert isinstance(index, int)
+    if get_option("diagonal_transformation_matrix"):
+        from dune.perftool.sumfact.switch import get_facedir, get_facemod
+        if index == get_facedir(Restriction.POSITIVE):
+            if get_facemod(Restriction.POSITIVE):
+                return 1, None
+            else:
+                return -1, None
+        else:
+            return 0, None
+    else:
+        from dune.perftool.pdelab.geometry import pymbolic_unit_outer_normal as _norm
+        return _norm(), visitor_indices
+
+
+def pymbolic_unit_outer_normal(visitor_indices):
+    index, = visitor_indices
+    assert isinstance(index, int)
+    if get_option("diagonal_transformation_matrix"):
+        from dune.perftool.sumfact.switch import get_facedir, get_facemod
+        if index == get_facedir(Restriction.NEGATIVE):
+            if get_facemod(Restriction.NEGATIVE):
+                return -1, None
+            else:
+                return 1, None
+        else:
+            return 0, None
+    else:
+        from dune.perftool.pdelab.geometry import pymbolic_unit_inner_normal as _norm
+        return _norm(), visitor_indices

python/dune/perftool/ufl/visitor.py

@@ -389,24 +389,14 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         # The normal must be restricted to be well-defined
         assert self.restriction is not Restriction.NONE
 
-        # Optimize facet normal on axiparallel grids
-        # TODO move this into the sumfact backend, it is only valid there
-        from dune.perftool.options import get_option
-        if get_option("diagonal_transformation_matrix") and get_option("sumfact"):
-            index, = self.indices
-            from dune.perftool.sumfact.switch import get_facedir
-            if isinstance(index, int) and index != get_facedir(self.restriction):
-                self.indices = None
-                return 0
-
         if self.restriction == Restriction.POSITIVE:
-            return Variable(self.interface.name_unit_outer_normal())
+            return self.interface.pymbolic_unit_outer_normal()
         if self.restriction == Restriction.NEGATIVE:
             # It is highly unnatural to have this generator function,
             # but I do run into subtle trouble with return -1*outer
             # as the indexing into the normal happens only later.
             # Not investing more time into this cornercase right now.
-            return Variable(self.interface.name_unit_inner_normal())
+            return self.interface.pymbolic_unit_inner_normal()
 
     def quadrature_weight(self, o):
         return self.interface.pymbolic_quadrature_weight()