From 7ca0782cba672b22c222a06c8aa32b116c8407e1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Thu, 13 Dec 2018 11:42:26 +0100
Subject: [PATCH] [skip ci] Choose more sophisticated quadrature order

---
 python/dune/codegen/sumfact/permutation.py | 37 ++++++++++++++--------
 1 file changed, 23 insertions(+), 14 deletions(-)

diff --git a/python/dune/codegen/sumfact/permutation.py b/python/dune/codegen/sumfact/permutation.py
index 0c98b368..1f775eba 100644
--- a/python/dune/codegen/sumfact/permutation.py
+++ b/python/dune/codegen/sumfact/permutation.py
@@ -4,6 +4,7 @@ import itertools
 
 from dune.codegen.options import get_option
 from dune.codegen.sumfact.switch import get_facedir, get_facemod
+from dune.codegen.ufl.modified_terminals import Restriction
 
 
 def sumfact_permutation_heuristic(permutations, stage):
@@ -116,26 +117,34 @@ def sumfact_quadrature_permutation_strategy(dim, restriction):
     # Use a simpler convention for structured grids. In this case we can always
     # go through the directions in the normal order. The same is true for 2D
     # and sum factorization on volumes.
-    if (not get_option('grid_unstructured')) or dim == 2 or restriction == 0:
+    if (not get_option('grid_unstructured')) or dim == 2 or restriction == Restriction.NONE:
         return tuple(range(dim))
     else:
+        # Draw a cube with edge orientations. We always do the normal direction
+        # first. The first case for facedir=0, facemod=0 was chosen that way,
+        # all others can be derived by rotating the cube and matching edge
+        # directions.
         def _order_on_self(restriction):
             facedir = get_facedir(restriction)
             facemod = get_facemod(restriction)
 
-            # Here we specify the convention
-            if (facedir, facemod) in [(0, 0), (1, 1), (2, 0)]:
-                return tuple(range(dim))
-            else:
-                l = list(range(dim))
-                l.reverse()
-                return tuple(l)
-
-        # On neighbor we need the reverse order
-        if restriction == 1:
+            quadrature_order = {
+                (0, 0): (0, 1, 2),
+                (0, 1): (0, 2, 1),
+                (1, 0): (1, 2, 0),
+                (1, 1): (1, 0, 2),
+                (2, 0): (2, 0, 1),
+                (2, 1): (2, 1, 0),
+            }
+
+            return quadrature_order[(facedir, facemod)]
+
+        # On neighbor we also do the normal direction first. The other two can
+        # be derived by putting a second edge oriented cube besides the first
+        # one. Then rotate and match edge directions.
+        if restriction == Restriction.POSITIVE:
             return _order_on_self(restriction)
         else:
-            assert restriction == 2
+            assert restriction == Restriction.NEGATIVE
             l = list(_order_on_self(restriction))
-            l.reverse()
-            return tuple(l)
+            return (l[0], l[2], l[1])
-- 
GitLab