From 90ea8f0bde386a71139c8e46f86fef83f20e5942 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 Apr 2024 13:15:36 -0500 Subject: [PATCH] Also test Gauss theorem on hex mesh from @anderson2981 x-ref: https://github.com/inducer/meshmode/issues/403 --- test/gh-403.msh | 305 ++++++++++++++++++++++++++++++++++++++++++++ test/test_grudge.py | 79 +++++++++--- 2 files changed, 364 insertions(+), 20 deletions(-) create mode 100644 test/gh-403.msh diff --git a/test/gh-403.msh b/test/gh-403.msh new file mode 100644 index 000000000..336a749cb --- /dev/null +++ b/test/gh-403.msh @@ -0,0 +1,305 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +8 +2 3 "inflow" +2 4 "outflow" +2 5 "flow" +2 6 "isothermal_wall" +2 7 "wall_interface" +2 8 "wall_farfield" +3 1 "fluid" +3 2 "wall_insert" +$EndPhysicalNames +$Nodes +117 +1 0 -0.01 0 +2 0.1 -0.01 0 +3 0.1 0.01 0 +4 0 0.01 0 +5 0.12 -0.01 0 +6 0.12 0.01 0 +7 0 -0.01 0.02 +8 0 0.01 0.02 +9 0.1 0.01 0.02 +10 0.1 -0.01 0.02 +11 0.12 0.01 0.02 +12 0.12 -0.01 0.02 +13 0.009999999999982485 -0.01 0 +14 0.01999999999995601 -0.01 0 +15 0.02999999999992731 -0.01 0 +16 0.0399999999998959 -0.01 0 +17 0.04999999999986855 -0.01 0 +18 0.0599999999998943 -0.01 0 +19 0.06999999999991978 -0.01 0 +20 0.07999999999994742 -0.01 0 +21 0.08999999999997288 -0.01 0 +22 0.1 -2.65987232239695e-14 0 +23 0.0900000000000228 0.01 0 +24 0.08000000000004992 0.01 0 +25 0.07000000000008301 0.01 0 +26 0.06000000000011067 0.01 0 +27 0.05000000000013697 0.01 0 +28 0.04000000000011013 0.01 0 +29 0.0300000000000834 0.01 0 +30 0.02000000000005492 0.01 0 +31 0.01000000000002775 0.01 0 +32 0 2.65987232239695e-14 0 +33 0.1099999999999756 0.01 0 +34 0.1099999999999756 -0.01 0 +35 0.12 -2.65987232239695e-14 0 +36 0 -2.65987232239695e-14 0.02 +37 0.009999999999982485 0.01 0.02 +38 0.01999999999995601 0.01 0.02 +39 0.02999999999992731 0.01 0.02 +40 0.0399999999998959 0.01 0.02 +41 0.04999999999986855 0.01 0.02 +42 0.0599999999998943 0.01 0.02 +43 0.06999999999991978 0.01 0.02 +44 0.07999999999994742 0.01 0.02 +45 0.08999999999997288 0.01 0.02 +46 0.1 2.65987232239695e-14 0.02 +47 0.0900000000000228 -0.01 0.02 +48 0.08000000000004992 -0.01 0.02 +49 0.07000000000008301 -0.01 0.02 +50 0.06000000000011067 -0.01 0.02 +51 0.05000000000013697 -0.01 0.02 +52 0.04000000000011013 -0.01 0.02 +53 0.0300000000000834 -0.01 0.02 +54 0.02000000000005492 -0.01 0.02 +55 0.01000000000002775 -0.01 0.02 +56 0 -0.01 0.009999999999973929 +57 0 0.01 0.009999999999973929 +58 0.1 0.01 0.009999999999973929 +59 0.1 -0.01 0.009999999999973929 +60 0.1099999999999756 0.01 0.02 +61 0.12 2.65987232239695e-14 0.02 +62 0.1100000000000244 -0.01 0.02 +63 0.12 0.01 0.009999999999973929 +64 0.12 -0.01 0.009999999999973929 +65 0.01000000000000377 7.12693062460239e-15 0 +66 0.02000000000000486 1.908999274446981e-15 0 +67 0.03000000000000474 5.090664731859039e-16 0 +68 0.04000000000000338 1.272666182964841e-16 0 +69 0.05000000000000278 -1.762611085103303e-29 0 +70 0.0600000000000022 -1.272666182965546e-16 0 +71 0.07000000000000106 -5.090664731861571e-16 0 +72 0.07999999999999929 -1.908999274447778e-15 0 +73 0.08999999999999873 -7.126930624604066e-15 0 +74 0.1099999999999878 -1.329936161198474e-14 0 +75 0 0 0.009999999999986964 +76 0.09000000000000011 0.01 0.009999999999979494 +77 0.08000000000000039 0.01 0.009999999999983416 +78 0.07000000000000083 0.01 0.009999999999986017 +79 0.06000000000000119 0.01 0.009999999999987496 +80 0.05000000000000147 0.01 0.009999999999987975 +81 0.04000000000000164 0.01 0.009999999999987496 +82 0.03000000000000169 0.01 0.009999999999986017 +83 0.02000000000000146 0.01 0.009999999999983416 +84 0.0100000000000009 0.01 0.009999999999979494 +85 0.1 0 0.009999999999986964 +86 0.01000000000000092 -0.01 0.009999999999979494 +87 0.02000000000000149 -0.01 0.009999999999983416 +88 0.03000000000000173 -0.01 0.009999999999986017 +89 0.04000000000000167 -0.01 0.009999999999987496 +90 0.05000000000000149 -0.01 0.009999999999987975 +91 0.06000000000000121 -0.01 0.009999999999987496 +92 0.07000000000000083 -0.01 0.009999999999986017 +93 0.0800000000000004 -0.01 0.009999999999983416 +94 0.09000000000000012 -0.01 0.009999999999979494 +95 0.01000000000000377 -7.12693062460239e-15 0.02 +96 0.02000000000000486 -1.908999274446981e-15 0.02 +97 0.03000000000000474 -5.090664731859039e-16 0.02 +98 0.04000000000000338 -1.272666182964841e-16 0.02 +99 0.05000000000000278 1.762611085103303e-29 0.02 +100 0.0600000000000022 1.272666182965546e-16 0.02 +101 0.07000000000000106 5.090664731861571e-16 0.02 +102 0.07999999999999929 1.908999274447778e-15 0.02 +103 0.08999999999999873 7.126930624604066e-15 0.02 +104 0.1099999999999878 0.01 0.009999999999986964 +105 0.12 0 0.009999999999986964 +106 0.11 -0.01 0.009999999999986964 +107 0.11 1.329936161198474e-14 0.02 +108 0.00999999999999956 -1.42247325030107e-19 0.009999999999992529 +109 0.02000000000000087 6.21898366137663e-19 0.009999999999996447 +110 0.0300000000000011 7.528699885738973e-19 0.009999999999999053 +111 0.04000000000000203 1.942890293093493e-19 0.01000000000000053 +112 0.05000000000000147 1.734723475976137e-19 0.01000000000000102 +113 0.06000000000000093 3.608224830031225e-19 0.01000000000000053 +114 0.07000000000000051 7.476658181459641e-19 0.009999999999999053 +115 0.08000000000000101 3.365363543394754e-19 0.009999999999996451 +116 0.09000000000000097 6.019490461639405e-19 0.00999999999999253 +117 0.1100000000000001 -4.770489558935599e-19 0.01 +$EndNodes +$Elements +168 +1 3 2 6 1 4 31 65 32 +2 3 2 6 1 31 30 66 65 +3 3 2 6 1 30 29 67 66 +4 3 2 6 1 29 28 68 67 +5 3 2 6 1 28 27 69 68 +6 3 2 6 1 27 26 70 69 +7 3 2 6 1 26 25 71 70 +8 3 2 6 1 25 24 72 71 +9 3 2 6 1 24 23 73 72 +10 3 2 6 1 23 3 22 73 +11 3 2 6 1 32 65 13 1 +12 3 2 6 1 65 66 14 13 +13 3 2 6 1 66 67 15 14 +14 3 2 6 1 67 68 16 15 +15 3 2 6 1 68 69 17 16 +16 3 2 6 1 69 70 18 17 +17 3 2 6 1 70 71 19 18 +18 3 2 6 1 71 72 20 19 +19 3 2 6 1 72 73 21 20 +20 3 2 6 1 73 22 2 21 +21 3 2 8 2 2 22 74 34 +22 3 2 8 2 34 74 35 5 +23 3 2 8 2 22 3 33 74 +24 3 2 8 2 74 33 6 35 +25 3 2 3 16 4 57 75 32 +26 3 2 5 16 4 57 75 32 +27 3 2 3 16 57 8 36 75 +28 3 2 5 16 57 8 36 75 +29 3 2 3 16 32 75 56 1 +30 3 2 5 16 32 75 56 1 +31 3 2 3 16 75 36 7 56 +32 3 2 5 16 75 36 7 56 +33 3 2 6 20 3 58 76 23 +34 3 2 6 20 58 9 45 76 +35 3 2 6 20 23 76 77 24 +36 3 2 6 20 76 45 44 77 +37 3 2 6 20 24 77 78 25 +38 3 2 6 20 77 44 43 78 +39 3 2 6 20 25 78 79 26 +40 3 2 6 20 78 43 42 79 +41 3 2 6 20 26 79 80 27 +42 3 2 6 20 79 42 41 80 +43 3 2 6 20 27 80 81 28 +44 3 2 6 20 80 41 40 81 +45 3 2 6 20 28 81 82 29 +46 3 2 6 20 81 40 39 82 +47 3 2 6 20 29 82 83 30 +48 3 2 6 20 82 39 38 83 +49 3 2 6 20 30 83 84 31 +50 3 2 6 20 83 38 37 84 +51 3 2 6 20 31 84 57 4 +52 3 2 6 20 84 37 8 57 +53 3 2 4 24 2 59 85 22 +54 3 2 5 24 2 59 85 22 +55 3 2 7 24 2 59 85 22 +56 3 2 4 24 59 10 46 85 +57 3 2 5 24 59 10 46 85 +58 3 2 7 24 59 10 46 85 +59 3 2 4 24 22 85 58 3 +60 3 2 5 24 22 85 58 3 +61 3 2 7 24 22 85 58 3 +62 3 2 4 24 85 46 9 58 +63 3 2 5 24 85 46 9 58 +64 3 2 7 24 85 46 9 58 +65 3 2 6 28 1 56 86 13 +66 3 2 6 28 56 7 55 86 +67 3 2 6 28 13 86 87 14 +68 3 2 6 28 86 55 54 87 +69 3 2 6 28 14 87 88 15 +70 3 2 6 28 87 54 53 88 +71 3 2 6 28 15 88 89 16 +72 3 2 6 28 88 53 52 89 +73 3 2 6 28 16 89 90 17 +74 3 2 6 28 89 52 51 90 +75 3 2 6 28 17 90 91 18 +76 3 2 6 28 90 51 50 91 +77 3 2 6 28 18 91 92 19 +78 3 2 6 28 91 50 49 92 +79 3 2 6 28 19 92 93 20 +80 3 2 6 28 92 49 48 93 +81 3 2 6 28 20 93 94 21 +82 3 2 6 28 93 48 47 94 +83 3 2 6 28 21 94 59 2 +84 3 2 6 28 94 47 10 59 +85 3 2 6 29 7 36 95 55 +86 3 2 6 29 55 95 96 54 +87 3 2 6 29 54 96 97 53 +88 3 2 6 29 53 97 98 52 +89 3 2 6 29 52 98 99 51 +90 3 2 6 29 51 99 100 50 +91 3 2 6 29 50 100 101 49 +92 3 2 6 29 49 101 102 48 +93 3 2 6 29 48 102 103 47 +94 3 2 6 29 47 103 46 10 +95 3 2 6 29 36 8 37 95 +96 3 2 6 29 95 37 38 96 +97 3 2 6 29 96 38 39 97 +98 3 2 6 29 97 39 40 98 +99 3 2 6 29 98 40 41 99 +100 3 2 6 29 99 41 42 100 +101 3 2 6 29 100 42 43 101 +102 3 2 6 29 101 43 44 102 +103 3 2 6 29 102 44 45 103 +104 3 2 6 29 103 45 9 46 +105 3 2 8 42 3 33 104 58 +106 3 2 8 42 58 104 60 9 +107 3 2 8 42 33 6 63 104 +108 3 2 8 42 104 63 11 60 +109 3 2 8 46 5 64 105 35 +110 3 2 8 46 64 12 61 105 +111 3 2 8 46 35 105 63 6 +112 3 2 8 46 105 61 11 63 +113 3 2 8 50 2 59 106 34 +114 3 2 8 50 59 10 62 106 +115 3 2 8 50 34 106 64 5 +116 3 2 8 50 106 62 12 64 +117 3 2 8 51 9 60 107 46 +118 3 2 8 51 60 11 61 107 +119 3 2 8 51 46 107 62 10 +120 3 2 8 51 107 61 12 62 +121 5 2 1 1 4 32 65 31 57 75 108 84 +122 5 2 1 1 57 75 108 84 8 36 95 37 +123 5 2 1 1 31 65 66 30 84 108 109 83 +124 5 2 1 1 84 108 109 83 37 95 96 38 +125 5 2 1 1 30 66 67 29 83 109 110 82 +126 5 2 1 1 83 109 110 82 38 96 97 39 +127 5 2 1 1 29 67 68 28 82 110 111 81 +128 5 2 1 1 82 110 111 81 39 97 98 40 +129 5 2 1 1 28 68 69 27 81 111 112 80 +130 5 2 1 1 81 111 112 80 40 98 99 41 +131 5 2 1 1 27 69 70 26 80 112 113 79 +132 5 2 1 1 80 112 113 79 41 99 100 42 +133 5 2 1 1 26 70 71 25 79 113 114 78 +134 5 2 1 1 79 113 114 78 42 100 101 43 +135 5 2 1 1 25 71 72 24 78 114 115 77 +136 5 2 1 1 78 114 115 77 43 101 102 44 +137 5 2 1 1 24 72 73 23 77 115 116 76 +138 5 2 1 1 77 115 116 76 44 102 103 45 +139 5 2 1 1 23 73 22 3 76 116 85 58 +140 5 2 1 1 76 116 85 58 45 103 46 9 +141 5 2 1 1 32 1 13 65 75 56 86 108 +142 5 2 1 1 75 56 86 108 36 7 55 95 +143 5 2 1 1 65 13 14 66 108 86 87 109 +144 5 2 1 1 108 86 87 109 95 55 54 96 +145 5 2 1 1 66 14 15 67 109 87 88 110 +146 5 2 1 1 109 87 88 110 96 54 53 97 +147 5 2 1 1 67 15 16 68 110 88 89 111 +148 5 2 1 1 110 88 89 111 97 53 52 98 +149 5 2 1 1 68 16 17 69 111 89 90 112 +150 5 2 1 1 111 89 90 112 98 52 51 99 +151 5 2 1 1 69 17 18 70 112 90 91 113 +152 5 2 1 1 112 90 91 113 99 51 50 100 +153 5 2 1 1 70 18 19 71 113 91 92 114 +154 5 2 1 1 113 91 92 114 100 50 49 101 +155 5 2 1 1 71 19 20 72 114 92 93 115 +156 5 2 1 1 114 92 93 115 101 49 48 102 +157 5 2 1 1 72 20 21 73 115 93 94 116 +158 5 2 1 1 115 93 94 116 102 48 47 103 +159 5 2 1 1 73 21 2 22 116 94 59 85 +160 5 2 1 1 116 94 59 85 103 47 10 46 +161 5 2 2 2 74 22 2 34 117 85 59 106 +162 5 2 2 2 117 85 59 106 107 46 10 62 +163 5 2 2 2 35 74 34 5 105 117 106 64 +164 5 2 2 2 105 117 106 64 61 107 62 12 +165 5 2 2 2 33 3 22 74 104 58 85 117 +166 5 2 2 2 104 58 85 117 60 9 46 107 +167 5 2 2 2 6 33 74 35 63 104 117 105 +168 5 2 2 2 63 104 117 105 11 60 107 61 +$EndElements \ No newline at end of file diff --git a/test/test_grudge.py b/test/test_grudge.py index 547ce8a2c..b902ee890 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -23,11 +23,13 @@ THE SOFTWARE. """ +from meshmode.mesh import TensorProductElementGroup import numpy as np import numpy.linalg as la from grudge.array_context import PytestPyOpenCLArrayContextFactory from arraycontext import pytest_generate_tests_for_array_contexts + pytest_generate_tests = pytest_generate_tests_for_array_contexts( [PytestPyOpenCLArrayContextFactory]) @@ -428,44 +430,81 @@ def df(x, axis): # {{{ divergence theorem -def test_2d_gauss_theorem(actx_factory): +@pytest.mark.parametrize("case", ["circle", "tp_box2", "tp_box3", "gh-403"]) +def test_gauss_theorem(actx_factory, case, visualize=False): """Verify Gauss's theorem explicitly on a mesh""" pytest.importorskip("meshpy") - from meshpy.geometry import make_circle, GeometryBuilder - from meshpy.triangle import MeshInfo, build - - geob = GeometryBuilder() - geob.add_geometry(*make_circle(1)) - mesh_info = MeshInfo() - geob.set(mesh_info) - - mesh_info = build(mesh_info) + if case == "circle": + from meshpy.geometry import make_circle, GeometryBuilder + from meshpy.triangle import MeshInfo, build + + geob = GeometryBuilder() + geob.add_geometry(*make_circle(1)) + mesh_info = MeshInfo() + geob.set(mesh_info) + + mesh_info = build(mesh_info) + + from meshmode.mesh.io import from_meshpy + mesh = from_meshpy(mesh_info, order=1) + elif case == "gh-403": + # https://github.com/inducer/meshmode/issues/403 + from meshmode.mesh.io import read_gmsh + mesh = read_gmsh("gh-403.msh") + elif case.startswith("tp_box"): + dim = int(case[-1]) + mesh = mgen.generate_regular_rect_mesh( + a=(-0.5,)*dim, + b=(0.5,)*dim, + nelements_per_axis=(4,)*dim, + group_cls=TensorProductElementGroup) - from meshmode.mesh.io import from_meshpy from meshmode.mesh import BTAG_ALL - mesh = from_meshpy(mesh_info, order=1) - actx = actx_factory() - dcoll = DiscretizationCollection(actx, mesh, order=2) - volm_disc = dcoll.discr_from_dd(dof_desc.DD_VOLUME) + dcoll = make_discretization_collection(actx, mesh, order=2) + volm_disc = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL) x_volm = actx.thaw(volm_disc.nodes()) def f(x): + if len(x) == 3: + x0, x1, x2 = x + elif len(x) == 2: + x0, x1 = x + x2 = 0 + else: + raise ValueError("unsupported dimensionality") + return flat_obj_array( - actx.np.sin(3*x[0]) + actx.np.cos(3*x[1]), - actx.np.sin(2*x[0]) + actx.np.cos(x[1]) - ) + actx.np.sin(3*x0) + actx.np.cos(3*x1) + 2*actx.np.cos(2*x2), + actx.np.sin(2*x0) + actx.np.cos(x1) + 4*actx.np.cos(0.5*x2), + actx.np.sin(1*x0) + actx.np.cos(2*x1) + 3*actx.np.cos(0.8*x2), + )[:dcoll.ambient_dim] f_volm = f(x_volm) - int_1 = op.integral(dcoll, "vol", op.local_div(dcoll, f_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) - int_2 = op.integral(dcoll, BTAG_ALL, prj_f.dot(normal)) + f_dot_n = prj_f.dot(normal) + int_2 = op.integral(dcoll, BTAG_ALL, f_dot_n) + + if visualize: + from grudge.shortcuts import make_visualizer, make_boundary_visualizer + vis = make_visualizer(dcoll) + bvis = make_boundary_visualizer(dcoll) + + vis.write_vtk_file( + f"gauss-thm-{case}-vol.vtu", [("div_f", div_f),]) + bvis.write_vtk_file( + f"gauss-thm-{case}-bdry.vtu", [ + ("f_dot_n", f_dot_n), + ("normal", normal), + ]) assert abs(int_1 - int_2) < 1e-13