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