From 52a3bb2de8aac95014e203a9bd75495e9bbe9ed3 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Fri, 24 May 2024 10:15:05 -0500
Subject: [PATCH] Add test that shows Gauss's law for low-order hexes may need
 overintegration

---
 test/gh-339.msh     | 356 ++++++++++++++++++++++++++++++++++++++++++++
 test/test_grudge.py |  56 +++++--
 2 files changed, 403 insertions(+), 9 deletions(-)
 create mode 100644 test/gh-339.msh

diff --git a/test/gh-339.msh b/test/gh-339.msh
new file mode 100644
index 00000000..1b271835
--- /dev/null
+++ b/test/gh-339.msh
@@ -0,0 +1,356 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+147
+1 0 0 1
+2 0 0 0
+3 0 1 1
+4 0 1 0
+5 1 0 1
+6 1 0 0
+7 1 1 1
+8 1 1 0
+9 0 0 0.4999999999999999
+10 0 0.4999999999999999 1
+11 0 1 0.4999999999999999
+12 0 0.4999999999999999 0
+13 1 0 0.4999999999999999
+14 1 0.4999999999999999 1
+15 1 1 0.4999999999999999
+16 1 0.4999999999999999 0
+17 0.4999999999999999 0 0
+18 0.4999999999999999 0 1
+19 0.4999999999999999 1 0
+20 0.4999999999999999 1 1
+21 0 0.5 0.5
+22 0 0.25 0.75
+23 0 0.25 0.25
+24 0 0.75 0.75
+25 0 0.75 0.25
+26 0 0.1666666666666667 0.5
+27 0 0.5 0.8333333333333334
+28 0 0.5 0.1666666666666667
+29 0 0.8333333333333334 0.5
+30 1 0.5 0.5
+31 1 0.25 0.25
+32 1 0.25 0.75
+33 1 0.75 0.75
+34 1 0.75 0.25
+35 1 0.1666666666666667 0.5
+36 1 0.5 0.8333333333333334
+37 1 0.5 0.1666666666666667
+38 1 0.8333333333333334 0.5
+39 0.5 0 0.5
+40 0.25 0 0.25
+41 0.25 0 0.75
+42 0.75 0 0.75
+43 0.75 0 0.25
+44 0.1666666666666667 0 0.5
+45 0.5 0 0.8333333333333334
+46 0.5 0 0.1666666666666667
+47 0.8333333333333334 0 0.5
+48 0.5 1 0.5
+49 0.25 1 0.75
+50 0.25 1 0.25
+51 0.75 1 0.75
+52 0.75 1 0.25
+53 0.1666666666666667 1 0.5
+54 0.5 1 0.8333333333333334
+55 0.5 1 0.1666666666666667
+56 0.8333333333333334 1 0.5
+57 0.5 0.5 0
+58 0.25 0.75 0
+59 0.25 0.25 0
+60 0.75 0.25 0
+61 0.75 0.75 0
+62 0.1666666666666667 0.5 0
+63 0.5 0.1666666666666667 0
+64 0.5 0.8333333333333334 0
+65 0.8333333333333334 0.5 0
+66 0.5 0.5 1
+67 0.25 0.25 1
+68 0.25 0.75 1
+69 0.75 0.25 1
+70 0.75 0.75 1
+71 0.1666666666666667 0.5 1
+72 0.5 0.1666666666666667 1
+73 0.5 0.8333333333333334 1
+74 0.8333333333333334 0.5 1
+75 0.75 0.5 0.25
+76 0.5 0.5 0.5
+77 0.25 0.5 0.25
+78 0.5 0.25 0.25
+79 0.25 0.25 0.5
+80 0.75 0.25 0.5
+81 0.25 0.5 0.75
+82 0.25 0.75 0.5
+83 0.5 0.75 0.75
+84 0.75 0.5 0.75
+85 0.75 0.75 0.5
+86 0.5 0.75 0.25
+87 0.5 0.25 0.75
+88 0.5 0.5 0.3333333333333333
+89 0.6666666666666666 0.3333333333333333 0.3333333333333333
+90 0.3333333333333333 0.3333333333333333 0.3333333333333333
+91 0.5 0.3333333333333333 0.5
+92 0.5 0.375 0.375
+93 0.3333333333333333 0.6666666666666666 0.6666666666666666
+94 0.5 0.5 0.6666666666666666
+95 0.6666666666666666 0.6666666666666666 0.6666666666666666
+96 0.5 0.6666666666666666 0.5
+97 0.5 0.625 0.625
+98 0.6666666666666666 0.6666666666666666 0.3333333333333333
+99 0.3333333333333333 0.6666666666666666 0.3333333333333333
+100 0.5 0.625 0.375
+101 0.3333333333333333 0.3333333333333333 0.6666666666666666
+102 0.6666666666666666 0.3333333333333333 0.6666666666666666
+103 0.5 0.375 0.625
+104 0.8333333333333334 0.8333333333333334 0.6666666666666666
+105 0.8333333333333334 0.8333333333333334 0.3333333333333333
+106 0.875 0.875 0.5
+107 0.1666666666666667 0.8333333333333334 0.3333333333333333
+108 0.1666666666666667 0.8333333333333334 0.6666666666666666
+109 0.125 0.875 0.5
+110 0.6666666666666666 0.1666666666666667 0.8333333333333334
+111 0.3333333333333333 0.1666666666666667 0.8333333333333334
+112 0.5 0.125 0.875
+113 0.1666666666666667 0.6666666666666666 0.8333333333333334
+114 0.1666666666666667 0.3333333333333333 0.8333333333333334
+115 0.125 0.5 0.875
+116 0.3333333333333333 0.8333333333333334 0.8333333333333334
+117 0.6666666666666666 0.8333333333333334 0.8333333333333334
+118 0.5 0.875 0.875
+119 0.8333333333333334 0.6666666666666666 0.1666666666666667
+120 0.8333333333333334 0.3333333333333333 0.1666666666666667
+121 0.875 0.5 0.125
+122 0.8333333333333334 0.6666666666666666 0.8333333333333334
+123 0.8333333333333334 0.3333333333333333 0.8333333333333334
+124 0.875 0.5 0.875
+125 0.1666666666666667 0.3333333333333333 0.1666666666666667
+126 0.1666666666666667 0.6666666666666666 0.1666666666666667
+127 0.125 0.5 0.125
+128 0.1666666666666667 0.1666666666666667 0.6666666666666666
+129 0.1666666666666667 0.1666666666666667 0.3333333333333333
+130 0.125 0.125 0.5
+131 0.3333333333333333 0.8333333333333334 0.1666666666666667
+132 0.6666666666666666 0.8333333333333334 0.1666666666666667
+133 0.5 0.875 0.125
+134 0.6666666666666666 0.1666666666666667 0.1666666666666667
+135 0.3333333333333333 0.1666666666666667 0.1666666666666667
+136 0.5 0.125 0.125
+137 0.8333333333333334 0.1666666666666667 0.3333333333333333
+138 0.8333333333333334 0.1666666666666667 0.6666666666666666
+139 0.875 0.125 0.5
+140 0.25 0.25 0.75
+141 0.25 0.75 0.75
+142 0.75 0.75 0.75
+143 0.25 0.75 0.25
+144 0.75 0.75 0.25
+145 0.75 0.25 0.25
+146 0.75 0.25 0.75
+147 0.25 0.25 0.25
+$EndNodes
+$Elements
+200
+1 15 2 0 1 1
+2 15 2 0 2 2
+3 15 2 0 3 3
+4 15 2 0 4 4
+5 15 2 0 5 5
+6 15 2 0 6 6
+7 15 2 0 7 7
+8 15 2 0 8 8
+9 1 2 0 1 2 9
+10 1 2 0 1 9 1
+11 1 2 0 2 1 10
+12 1 2 0 2 10 3
+13 1 2 0 3 4 11
+14 1 2 0 3 11 3
+15 1 2 0 4 2 12
+16 1 2 0 4 12 4
+17 1 2 0 5 6 13
+18 1 2 0 5 13 5
+19 1 2 0 6 5 14
+20 1 2 0 6 14 7
+21 1 2 0 7 8 15
+22 1 2 0 7 15 7
+23 1 2 0 8 6 16
+24 1 2 0 8 16 8
+25 1 2 0 9 2 17
+26 1 2 0 9 17 6
+27 1 2 0 10 1 18
+28 1 2 0 10 18 5
+29 1 2 0 11 4 19
+30 1 2 0 11 19 8
+31 1 2 0 12 3 20
+32 1 2 0 12 20 7
+33 3 2 0 1 2 9 26 23
+34 3 2 0 1 9 1 22 26
+35 3 2 0 1 23 26 22 21
+36 3 2 0 1 1 10 27 22
+37 3 2 0 1 10 3 24 27
+38 3 2 0 1 22 27 24 21
+39 3 2 0 1 4 12 28 25
+40 3 2 0 1 12 2 23 28
+41 3 2 0 1 25 28 23 21
+42 3 2 0 1 3 11 29 24
+43 3 2 0 1 11 4 25 29
+44 3 2 0 1 24 29 25 21
+45 3 2 0 2 6 31 35 13
+46 3 2 0 2 31 30 32 35
+47 3 2 0 2 13 35 32 5
+48 3 2 0 2 5 32 36 14
+49 3 2 0 2 32 30 33 36
+50 3 2 0 2 14 36 33 7
+51 3 2 0 2 8 34 37 16
+52 3 2 0 2 34 30 31 37
+53 3 2 0 2 16 37 31 6
+54 3 2 0 2 7 33 38 15
+55 3 2 0 2 33 30 34 38
+56 3 2 0 2 15 38 34 8
+57 3 2 0 3 1 9 44 41
+58 3 2 0 3 9 2 40 44
+59 3 2 0 3 41 44 40 39
+60 3 2 0 3 5 18 45 42
+61 3 2 0 3 18 1 41 45
+62 3 2 0 3 42 45 41 39
+63 3 2 0 3 2 17 46 40
+64 3 2 0 3 17 6 43 46
+65 3 2 0 3 40 46 43 39
+66 3 2 0 3 6 13 47 43
+67 3 2 0 3 13 5 42 47
+68 3 2 0 3 43 47 42 39
+69 3 2 0 4 3 49 53 11
+70 3 2 0 4 49 48 50 53
+71 3 2 0 4 11 53 50 4
+72 3 2 0 4 7 51 54 20
+73 3 2 0 4 51 48 49 54
+74 3 2 0 4 20 54 49 3
+75 3 2 0 4 4 50 55 19
+76 3 2 0 4 50 48 52 55
+77 3 2 0 4 19 55 52 8
+78 3 2 0 4 8 52 56 15
+79 3 2 0 4 52 48 51 56
+80 3 2 0 4 15 56 51 7
+81 3 2 0 5 2 12 62 59
+82 3 2 0 5 12 4 58 62
+83 3 2 0 5 59 62 58 57
+84 3 2 0 5 6 17 63 60
+85 3 2 0 5 17 2 59 63
+86 3 2 0 5 60 63 59 57
+87 3 2 0 5 4 19 64 58
+88 3 2 0 5 19 8 61 64
+89 3 2 0 5 58 64 61 57
+90 3 2 0 5 8 16 65 61
+91 3 2 0 5 16 6 60 65
+92 3 2 0 5 61 65 60 57
+93 3 2 0 6 1 67 71 10
+94 3 2 0 6 67 66 68 71
+95 3 2 0 6 10 71 68 3
+96 3 2 0 6 5 69 72 18
+97 3 2 0 6 69 66 67 72
+98 3 2 0 6 18 72 67 1
+99 3 2 0 6 3 68 73 20
+100 3 2 0 6 68 66 70 73
+101 3 2 0 6 20 73 70 7
+102 3 2 0 6 7 70 74 14
+103 3 2 0 6 70 66 69 74
+104 3 2 0 6 14 74 69 5
+105 5 2 0 1 57 75 88 77 78 89 92 90
+106 5 2 0 1 75 30 76 88 89 80 91 92
+107 5 2 0 1 77 88 76 21 90 92 91 79
+108 5 2 0 1 39 80 89 78 79 91 92 90
+109 5 2 0 1 66 81 93 83 84 94 97 95
+110 5 2 0 1 81 21 82 93 94 76 96 97
+111 5 2 0 1 83 93 82 48 95 97 96 85
+112 5 2 0 1 30 76 94 84 85 96 97 95
+113 5 2 0 1 30 76 96 85 75 88 100 98
+114 5 2 0 1 76 21 82 96 88 77 99 100
+115 5 2 0 1 85 96 82 48 98 100 99 86
+116 5 2 0 1 57 77 88 75 86 99 100 98
+117 5 2 0 1 66 81 94 84 87 101 103 102
+118 5 2 0 1 81 21 76 94 101 79 91 103
+119 5 2 0 1 84 94 76 30 102 103 91 80
+120 5 2 0 1 39 79 101 87 80 91 103 102
+121 5 2 0 1 48 51 104 85 52 56 106 105
+122 5 2 0 1 51 7 33 104 56 15 38 106
+123 5 2 0 1 85 104 33 30 105 106 38 34
+124 5 2 0 1 8 15 56 52 34 38 106 105
+125 5 2 0 1 21 82 107 25 24 108 109 29
+126 5 2 0 1 82 48 50 107 108 49 53 109
+127 5 2 0 1 25 107 50 4 29 109 53 11
+128 5 2 0 1 3 49 108 24 11 53 109 29
+129 5 2 0 1 5 18 72 69 42 45 112 110
+130 5 2 0 1 18 1 67 72 45 41 111 112
+131 5 2 0 1 69 72 67 66 110 112 111 87
+132 5 2 0 1 39 41 45 42 87 111 112 110
+133 5 2 0 1 66 68 113 81 67 71 115 114
+134 5 2 0 1 68 3 24 113 71 10 27 115
+135 5 2 0 1 81 113 24 21 114 115 27 22
+136 5 2 0 1 1 10 71 67 22 27 115 114
+137 5 2 0 1 48 49 116 83 51 54 118 117
+138 5 2 0 1 49 3 68 116 54 20 73 118
+139 5 2 0 1 83 116 68 66 117 118 73 70
+140 5 2 0 1 7 20 54 51 70 73 118 117
+141 5 2 0 1 57 61 119 75 60 65 121 120
+142 5 2 0 1 61 8 34 119 65 16 37 121
+143 5 2 0 1 75 119 34 30 120 121 37 31
+144 5 2 0 1 6 16 65 60 31 37 121 120
+145 5 2 0 1 7 14 74 70 33 36 124 122
+146 5 2 0 1 14 5 69 74 36 32 123 124
+147 5 2 0 1 70 74 69 66 122 124 123 84
+148 5 2 0 1 30 32 36 33 84 123 124 122
+149 5 2 0 1 21 23 28 25 77 125 127 126
+150 5 2 0 1 23 2 12 28 125 59 62 127
+151 5 2 0 1 25 28 12 4 126 127 62 58
+152 5 2 0 1 57 59 125 77 58 62 127 126
+153 5 2 0 1 1 9 26 22 41 44 130 128
+154 5 2 0 1 9 2 23 26 44 40 129 130
+155 5 2 0 1 22 26 23 21 128 130 129 79
+156 5 2 0 1 39 40 44 41 79 129 130 128
+157 5 2 0 1 48 50 55 52 86 131 133 132
+158 5 2 0 1 50 4 19 55 131 58 64 133
+159 5 2 0 1 52 55 19 8 132 133 64 61
+160 5 2 0 1 57 58 131 86 61 64 133 132
+161 5 2 0 1 6 60 63 17 43 134 136 46
+162 5 2 0 1 60 57 59 63 134 78 135 136
+163 5 2 0 1 17 63 59 2 46 136 135 40
+164 5 2 0 1 39 78 134 43 40 135 136 46
+165 5 2 0 1 6 13 35 31 43 47 139 137
+166 5 2 0 1 13 5 32 35 47 42 138 139
+167 5 2 0 1 31 35 32 30 137 139 138 80
+168 5 2 0 1 39 42 47 43 80 138 139 137
+169 5 2 0 1 21 81 114 22 79 101 140 128
+170 5 2 0 1 81 66 67 114 101 87 111 140
+171 5 2 0 1 22 114 67 1 128 140 111 41
+172 5 2 0 1 39 87 101 79 41 111 140 128
+173 5 2 0 1 3 24 108 49 68 113 141 116
+174 5 2 0 1 24 21 82 108 113 81 93 141
+175 5 2 0 1 49 108 82 48 116 141 93 83
+176 5 2 0 1 66 81 113 68 83 93 141 116
+177 5 2 0 1 66 83 117 70 84 95 142 122
+178 5 2 0 1 83 48 51 117 95 85 104 142
+179 5 2 0 1 70 117 51 7 122 142 104 33
+180 5 2 0 1 30 85 95 84 33 104 142 122
+181 5 2 0 1 4 50 107 25 58 131 143 126
+182 5 2 0 1 50 48 82 107 131 86 99 143
+183 5 2 0 1 25 107 82 21 126 143 99 77
+184 5 2 0 1 57 86 131 58 77 99 143 126
+185 5 2 0 1 8 34 105 52 61 119 144 132
+186 5 2 0 1 34 30 85 105 119 75 98 144
+187 5 2 0 1 52 105 85 48 132 144 98 86
+188 5 2 0 1 57 75 119 61 86 98 144 132
+189 5 2 0 1 30 75 120 31 80 89 145 137
+190 5 2 0 1 75 57 60 120 89 78 134 145
+191 5 2 0 1 31 120 60 6 137 145 134 43
+192 5 2 0 1 39 78 89 80 43 134 145 137
+193 5 2 0 1 30 32 123 84 80 138 146 102
+194 5 2 0 1 32 5 69 123 138 42 110 146
+195 5 2 0 1 84 123 69 66 102 146 110 87
+196 5 2 0 1 39 42 138 80 87 110 146 102
+197 5 2 0 1 21 23 125 77 79 129 147 90
+198 5 2 0 1 23 2 59 125 129 40 135 147
+199 5 2 0 1 77 125 59 57 90 147 135 78
+200 5 2 0 1 39 40 129 79 78 135 147 90
+$EndElements
diff --git a/test/test_grudge.py b/test/test_grudge.py
index e7bb7bdb..3420ee95 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -23,10 +23,12 @@
 THE SOFTWARE.
 """
 
-from meshmode.mesh import TensorProductElementGroup
 import numpy as np
 import numpy.linalg as la
 
+from meshmode.discretization.poly_element import \
+    InterpolatoryEdgeClusteredGroupFactory
+from meshmode.mesh import TensorProductElementGroup
 from grudge.array_context import PytestPyOpenCLArrayContextFactory
 from arraycontext import pytest_generate_tests_for_array_contexts
 
@@ -430,12 +432,16 @@ def df(x, axis):
 
 # {{{ divergence theorem
 
-@pytest.mark.parametrize("case", ["circle", "tp_box2", "tp_box3", "gh-403"])
+@pytest.mark.parametrize(
+             "case", ["circle", "tp_box2", "tp_box3", "gh-403", "gh-339"])
 def test_gauss_theorem(actx_factory, case, visualize=False):
     """Verify Gauss's theorem explicitly on a mesh"""
 
     pytest.importorskip("meshpy")
 
+    order = 2
+    use_overint = False
+
     if case == "circle":
         from meshpy.geometry import make_circle, GeometryBuilder
         from meshpy.triangle import MeshInfo, build
@@ -455,6 +461,13 @@ def test_gauss_theorem(actx_factory, case, visualize=False):
         from meshmode.mesh.io import read_gmsh
         mesh = read_gmsh("gh-403.msh")
 
+    elif case == "gh-339":
+        # https://github.com/inducer/grudge/issues/339
+        from meshmode.mesh.io import read_gmsh
+        mesh = read_gmsh("gh-339.msh")
+        order = 1
+        use_overint = True
+
     elif case.startswith("tp_box"):
         dim = int(case[-1])
         mesh = mgen.generate_regular_rect_mesh(
@@ -467,10 +480,17 @@ def test_gauss_theorem(actx_factory, case, visualize=False):
         raise ValueError(f"unknown case: {case}")
 
     from meshmode.mesh import BTAG_ALL
+    from meshmode.discretization.poly_element import QuadratureGroupFactory
 
     actx = actx_factory()
 
-    dcoll = make_discretization_collection(actx, mesh, order=2)
+    dcoll = make_discretization_collection(
+               actx, mesh, discr_tag_to_group_factory={
+                   dof_desc.DISCR_TAG_BASE: (
+                       InterpolatoryEdgeClusteredGroupFactory(order)),
+                   dof_desc.DISCR_TAG_QUAD: (
+                       QuadratureGroupFactory(order))
+               })
     volm_disc = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
     x_volm = actx.thaw(volm_disc.nodes())
 
@@ -490,13 +510,31 @@ def f(x):
         )[:dcoll.ambient_dim]
 
     f_volm = f(x_volm)
-    div_f = op.local_div(dcoll, f_volm)
-    int_1 = op.integral(dcoll, "vol", div_f)
 
-    prj_f = op.project(dcoll, "vol", BTAG_ALL, f_volm)
-    normal = geo.normal(actx, dcoll, BTAG_ALL)
-    f_dot_n = prj_f.dot(normal)
-    int_2 = op.integral(dcoll, BTAG_ALL, f_dot_n)
+    if not use_overint:
+        div_f = op.local_div(dcoll, f_volm)
+        int_1 = op.integral(dcoll, "vol", div_f)
+
+        prj_f = op.project(dcoll, "vol", BTAG_ALL, f_volm)
+        normal = geo.normal(actx, dcoll, BTAG_ALL)
+        f_dot_n = prj_f.dot(normal)
+        int_2 = op.integral(dcoll, BTAG_ALL, f_dot_n)
+    else:
+        dd_base_vol = dof_desc.as_dofdesc(
+                            dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_BASE)
+        dd_quad_vol = dof_desc.as_dofdesc(
+                            dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_QUAD)
+        dd_quad_bd = dof_desc.as_dofdesc(BTAG_ALL, dof_desc.DISCR_TAG_QUAD)
+
+        div_f = op.local_div(
+                 dcoll, dd_quad_vol,
+                 op.project(dcoll, dd_base_vol, dd_quad_vol, f_volm))
+        int_1 = op.integral(dcoll, dd_quad_vol, div_f)
+
+        prj_f = op.project(dcoll, "vol", dd_quad_bd, f_volm)
+        normal = geo.normal(actx, dcoll, dd_quad_bd)
+        f_dot_n = prj_f.dot(normal)
+        int_2 = op.integral(dcoll, dd_quad_bd, f_dot_n)
 
     if visualize:
         from grudge.shortcuts import make_visualizer, make_boundary_visualizer