--- /dev/null
+<?xml version="1.0" encoding="utf-8" standalone="no"?>
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
+ "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<!-- Created with matplotlib (https://matplotlib.org/) -->
+<svg height="140.860469pt" version="1.1" viewBox="0 0 375.02 140.860469" width="375.02pt" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
+ <metadata>
+ <rdf:RDF xmlns:cc="http://creativecommons.org/ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
+ <cc:Work>
+ <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/>
+ <dc:date>2021-05-19T11:42:51.997656</dc:date>
+ <dc:format>image/svg+xml</dc:format>
+ <dc:creator>
+ <cc:Agent>
+ <dc:title>Matplotlib v3.3.2, https://matplotlib.org/</dc:title>
+ </cc:Agent>
+ </dc:creator>
+ </cc:Work>
+ </rdf:RDF>
+ </metadata>
+ <defs>
+ <style type="text/css">*{stroke-linecap:butt;stroke-linejoin:round;}</style>
+ </defs>
+ <g id="figure_1">
+ <g id="patch_1">
+ <path d="M 0 140.860469
+L 375.02 140.860469
+L 375.02 0
+L 0 0
+z
+" style="fill:#ffffff;"/>
+ </g>
+ <g id="axes_1">
+ <g id="patch_2">
+ <path d="M 10.7 130.160469
+L 115.735294 130.160469
+L 115.735294 25.125175
+L 10.7 25.125175
+z
+" style="fill:#ffffff;"/>
+ </g>
+ <g id="line2d_1">
+ <path clip-path="url(#p6c1199b75e)" d="M 101.433173 141.860469
+L 99.818309 137.039698
+L 97.728225 131.467643
+L 95.4292 125.978497
+L 92.924514 120.58009
+L 90.217739 115.280126
+L 87.312739 110.086166
+L 84.213658 105.00562
+L 80.924918 100.045738
+L 77.45121 95.213595
+L 73.797491 90.516087
+L 69.968974 85.959914
+L 65.971122 81.551579
+L 61.809637 77.29737
+L 57.490458 73.203357
+L 53.019748 69.275381
+L 48.403883 65.519047
+L 43.649452 61.939714
+L 38.763236 58.542489
+L 33.752207 55.332218
+L 28.623515 52.313483
+L 23.384477 49.490589
+L 18.042567 46.867565
+L 12.605409 44.448153
+L 7.080758 42.235805
+L 1.476498 40.233678
+L -1 39.426889
+L -1 39.426889
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:1.5;"/>
+ </g>
+ <g id="matplotlib.axis_1"/>
+ <g id="matplotlib.axis_2"/>
+ <g id="PathCollection_1">
+ <defs>
+ <path d="M 0 2.12132
+C 0.562581 2.12132 1.102195 1.897805 1.5 1.5
+C 1.897805 1.102195 2.12132 0.562581 2.12132 0
+C 2.12132 -0.562581 1.897805 -1.102195 1.5 -1.5
+C 1.102195 -1.897805 0.562581 -2.12132 0 -2.12132
+C -0.562581 -2.12132 -1.102195 -1.897805 -1.5 -1.5
+C -1.897805 -1.102195 -2.12132 -0.562581 -2.12132 0
+C -2.12132 0.562581 -1.897805 1.102195 -1.5 1.5
+C -1.102195 1.897805 -0.562581 2.12132 0 2.12132
+z
+" id="m93c5f5d0fb" style="stroke:#1f77b4;"/>
+ </defs>
+ <g clip-path="url(#p6c1199b75e)">
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="27.108404" xlink:href="#m93c5f5d0fb" y="111.880337"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="71.937047" xlink:href="#m93c5f5d0fb" y="111.880337"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="17.662506" xlink:href="#m93c5f5d0fb" y="61.937945"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="36.684471" xlink:href="#m93c5f5d0fb" y="61.937945"/>
+ </g>
+ </g>
+ <g id="patch_3">
+ <path d="M 10.7 130.160469
+L 10.7 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_4">
+ <path d="M 115.735294 130.160469
+L 115.735294 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_5">
+ <path d="M 10.7 130.160469
+L 115.735294 130.160469
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_6">
+ <path d="M 10.7 25.125175
+L 115.735294 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="text_1">
+ <!-- $\psi < 0$ -->
+ <g transform="translate(42.470612 19.125175)scale(0.168 -0.168)">
+ <defs>
+ <path d="M 47 67.09375
+C 47 67.203125 47.40625 68.5 47.40625 68.609375
+C 47.40625 69.5 46.59375 69.5 46.296875 69.5
+C 45.40625 69.5 45.296875 69 44.90625 67.5
+L 28.40625 1.203125
+C 20.09375 2.203125 17 6.40625 17 12.40625
+C 17 14.59375 17 16.90625 21.703125 29.296875
+C 23 32.890625 23.5 34.1875 23.5 35.984375
+C 23.5 40.484375 20.296875 44.078125 15.59375 44.078125
+C 6.40625 44.078125 2.703125 29.59375 2.703125 28.796875
+C 2.703125 28.390625 3.09375 27.890625 3.796875 27.890625
+C 4.703125 27.890625 4.796875 28.296875 5.203125 29.6875
+C 7.59375 38.328125 11.59375 42 15.296875 42
+C 16.203125 42 17.90625 42 17.90625 38.796875
+C 17.90625 38.390625 17.90625 36.1875 16.203125 31.796875
+C 10.796875 17.59375 10.796875 15.390625 10.796875 13.09375
+C 10.796875 3.484375 18.796875 -0.21875 27.796875 -0.90625
+C 27 -3.953125 26.296875 -7.078125 25.5 -10.109375
+C 23.90625 -15.96875 23.203125 -18.71875 23.203125 -19.109375
+C 23.203125 -20 24 -20 24.296875 -20
+C 24.5 -20 24.90625 -20 25.09375 -19.59375
+C 25.5 -19.21875 29.59375 -2.75 30 -1
+C 33.703125 -1 41.59375 -1 50.59375 8.296875
+C 53.90625 11.890625 56.90625 16.5 58.59375 20.796875
+C 59.59375 23.40625 62 32.296875 62 37.390625
+C 62 43.390625 59 44.078125 58 44.078125
+C 55.59375 44.078125 53.40625 41.6875 53.40625 39.6875
+C 53.40625 38.484375 54.09375 37.78125 54.5 37.390625
+C 55.40625 36.484375 58.09375 33.78125 58.09375 28.59375
+C 58.09375 25 56.09375 17.59375 49.703125 10.40625
+C 41.296875 1 33.59375 1 30.59375 1
+z
+" id="CMMI12-32"/>
+ <path d="M 65.90625 48.5
+C 67.703125 49.296875 67.90625 50 67.90625 50.59375
+C 67.90625 51.703125 67.09375 52.5 66 52.5
+C 65.796875 52.5 65.703125 52.390625 64.296875 51.796875
+L 10.203125 26.859375
+C 8.40625 26.046875 8.203125 25.34375 8.203125 24.75
+C 8.203125 24.046875 8.296875 23.453125 10.203125 22.53125
+L 64.296875 -2.40625
+C 65.59375 -3 65.796875 -3 66 -3
+C 67.09375 -3 67.90625 -2.203125 67.90625 -1.09375
+C 67.90625 -0.5 67.703125 0.203125 65.90625 1
+L 14.40625 24.703125
+z
+" id="CMMI12-60"/>
+ <path d="M 42 31.640625
+C 42 37.75 41.90625 48.125 37.703125 56.109375
+C 34 63.109375 28.09375 65.59375 22.90625 65.59375
+C 18.09375 65.59375 12 63.40625 8.203125 56.203125
+C 4.203125 48.71875 3.796875 39.4375 3.796875 31.640625
+C 3.796875 25.953125 3.90625 17.28125 7 9.671875
+C 11.296875 -0.609375 19 -2 22.90625 -2
+C 27.5 -2 34.5 -0.109375 38.59375 9.375
+C 41.59375 16.28125 42 24.359375 42 31.640625
+z
+M 22.90625 -0.40625
+C 16.5 -0.40625 12.703125 5.078125 11.296875 12.6875
+C 10.203125 18.5625 10.203125 27.15625 10.203125 32.75
+C 10.203125 40.4375 10.203125 46.828125 11.5 52.921875
+C 13.40625 61.390625 19 64 22.90625 64
+C 27 64 32.296875 61.296875 34.203125 53.125
+C 35.5 47.4375 35.59375 40.734375 35.59375 32.75
+C 35.59375 26.25 35.59375 18.265625 34.40625 12.375
+C 32.296875 1.484375 26.40625 -0.40625 22.90625 -0.40625
+z
+" id="CMR17-48"/>
+ </defs>
+ <use transform="scale(0.996264)" xlink:href="#CMMI12-32"/>
+ <use transform="translate(94.667488 0)scale(0.996264)" xlink:href="#CMMI12-60"/>
+ <use transform="translate(198.213593 0)scale(0.996264)" xlink:href="#CMR17-48"/>
+ </g>
+ </g>
+ </g>
+ <g id="axes_2">
+ <g id="patch_7">
+ <path d="M 136.742353 130.160469
+L 241.777647 130.160469
+L 241.777647 25.125175
+L 136.742353 25.125175
+z
+" style="fill:#ffffff;"/>
+ </g>
+ <g id="line2d_2">
+ <path clip-path="url(#p31b09fbec5)" d="M 227.475526 141.860469
+L 225.860662 137.039698
+L 223.770578 131.467643
+L 221.471553 125.978497
+L 218.966866 120.58009
+L 216.260092 115.280126
+L 213.355092 110.086166
+L 210.256011 105.00562
+L 206.967271 100.045738
+L 203.493563 95.213595
+L 199.839844 90.516087
+L 196.011327 85.959914
+L 192.013475 81.551579
+L 187.85199 77.29737
+L 183.532811 73.203357
+L 179.062101 69.275381
+L 174.446236 65.519047
+L 169.691805 61.939714
+L 164.805589 58.542489
+L 159.79456 55.332218
+L 154.665868 52.313483
+L 149.42683 49.490589
+L 144.08492 46.867565
+L 138.647762 44.448153
+L 133.123111 42.235805
+L 127.518851 40.233678
+L 121.842977 38.444627
+L 116.103587 36.871206
+L 110.308871 35.51566
+L 104.467095 34.379921
+L 98.586594 33.465612
+L 92.675759 32.774036
+L 86.743023 32.30618
+L 80.796851 32.062711
+L 74.845725 32.043978
+L 68.898137 32.250006
+L 62.962573 32.680502
+L 57.047501 33.334851
+L 51.161361 34.21212
+L 45.312551 35.311058
+L 39.509415 36.630096
+L 33.760233 38.167352
+L 28.073208 39.920633
+L 22.456454 41.887438
+L 16.917985 44.06496
+L 11.465702 46.450093
+L 6.107385 49.039434
+L 0.850678 51.829288
+L -1 52.874084
+L -1 52.874084
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:1.5;"/>
+ </g>
+ <g id="matplotlib.axis_3"/>
+ <g id="matplotlib.axis_4"/>
+ <g id="PathCollection_2">
+ <g clip-path="url(#p31b09fbec5)">
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="158.938936" xlink:href="#m93c5f5d0fb" y="39.741361"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="219.581064" xlink:href="#m93c5f5d0fb" y="39.741361"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="158.938936" xlink:href="#m93c5f5d0fb" y="29.041521"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="219.581064" xlink:href="#m93c5f5d0fb" y="29.041521"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="220.175878" xlink:href="#m93c5f5d0fb" y="111.880337"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="235.989467" xlink:href="#m93c5f5d0fb" y="111.880337"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="184.923303" xlink:href="#m93c5f5d0fb" y="61.937945"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="226.543538" xlink:href="#m93c5f5d0fb" y="61.937945"/>
+ </g>
+ </g>
+ <g id="patch_8">
+ <path d="M 136.742353 130.160469
+L 136.742353 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_9">
+ <path d="M 241.777647 130.160469
+L 241.777647 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_10">
+ <path d="M 136.742353 130.160469
+L 241.777647 130.160469
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_11">
+ <path d="M 136.742353 25.125175
+L 241.777647 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="text_2">
+ <!-- $\psi > 0$ -->
+ <g transform="translate(168.512965 19.125175)scale(0.168 -0.168)">
+ <defs>
+ <path d="M 65.90625 22.640625
+C 67.796875 23.546875 67.90625 24.140625 67.90625 24.84375
+C 67.90625 25.4375 67.703125 26.140625 65.90625 26.953125
+L 11.796875 51.90625
+C 10.5 52.5 10.296875 52.5 10.09375 52.5
+C 8.90625 52.5 8.203125 51.5 8.203125 50.703125
+C 8.203125 49.5 9 49.09375 10.296875 48.5
+L 61.703125 24.84375
+L 10.203125 1.015625
+C 8.203125 0.09375 8.203125 -0.609375 8.203125 -1.203125
+C 8.203125 -2 8.90625 -3 10.09375 -3
+C 10.296875 -3 10.40625 -2.890625 11.796875 -2.296875
+z
+" id="CMMI12-62"/>
+ </defs>
+ <use transform="scale(0.996264)" xlink:href="#CMMI12-32"/>
+ <use transform="translate(94.667488 0)scale(0.996264)" xlink:href="#CMMI12-62"/>
+ <use transform="translate(198.213593 0)scale(0.996264)" xlink:href="#CMR17-48"/>
+ </g>
+ </g>
+ </g>
+ <g id="axes_3">
+ <g id="patch_12">
+ <path d="M 262.784706 130.160469
+L 367.82 130.160469
+L 367.82 25.125175
+L 262.784706 25.125175
+z
+" style="fill:#ffffff;"/>
+ </g>
+ <g id="line2d_3">
+ <path clip-path="url(#p368b104918)" d="M 353.517879 141.860469
+L 351.903015 137.039698
+L 349.812931 131.467643
+L 347.513906 125.978497
+L 345.009219 120.58009
+L 342.302445 115.280126
+L 339.397445 110.086166
+L 336.298364 105.00562
+L 333.009624 100.045738
+L 329.535916 95.213595
+L 325.882197 90.516087
+L 322.05368 85.959914
+L 318.055828 81.551579
+L 313.894343 77.29737
+L 309.575164 73.203357
+L 305.104454 69.275381
+L 300.488589 65.519047
+L 295.734158 61.939714
+L 290.847941 58.542489
+L 285.836913 55.332218
+L 280.708221 52.313483
+L 275.469182 49.490589
+L 270.127273 46.867565
+L 264.690115 44.448153
+L 259.165464 42.235805
+L 253.561204 40.233678
+L 247.88533 38.444627
+L 242.14594 36.871206
+L 236.351224 35.51566
+L 230.509448 34.379921
+L 224.628947 33.465612
+L 218.718112 32.774036
+L 212.785376 32.30618
+L 206.839204 32.062711
+L 200.888078 32.043978
+L 194.94049 32.250006
+L 189.004926 32.680502
+L 183.089854 33.334851
+L 177.203714 34.21212
+L 171.354904 35.311058
+L 165.551768 36.630096
+L 159.802586 38.167352
+L 154.115561 39.920633
+L 148.498807 41.887438
+L 142.960338 44.06496
+L 137.508055 46.450093
+L 132.149738 49.039434
+L 126.893031 51.829288
+L 121.745436 54.815675
+L 116.714295 57.994334
+L 111.806788 61.360729
+L 107.029916 64.910059
+L 102.390494 68.637258
+L 97.895143 72.53701
+L 93.550274 76.603749
+L 89.362089 80.831674
+L 85.336562 85.214753
+L 81.479437 89.746731
+L 77.796216 94.421144
+L 74.292155 99.231321
+L 70.972254 104.1704
+L 67.841248 109.231334
+L 64.903606 114.406902
+L 62.163518 119.68972
+L 59.624895 125.072251
+L 57.291357 130.546815
+L 55.166234 136.105601
+L 53.215214 141.860469
+L 53.215214 141.860469
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:1.5;"/>
+ </g>
+ <g id="matplotlib.axis_5"/>
+ <g id="matplotlib.axis_6"/>
+ <g id="PathCollection_3">
+ <g clip-path="url(#p368b104918)">
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="340.430156" xlink:href="#m93c5f5d0fb" y="111.880337"/>
+ <use style="fill:#1f77b4;stroke:#1f77b4;" x="295.731652" xlink:href="#m93c5f5d0fb" y="61.937945"/>
+ </g>
+ </g>
+ <g id="patch_13">
+ <path d="M 262.784706 130.160469
+L 262.784706 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_14">
+ <path d="M 367.82 130.160469
+L 367.82 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_15">
+ <path d="M 262.784706 130.160469
+L 367.82 130.160469
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="patch_16">
+ <path d="M 262.784706 25.125175
+L 367.82 25.125175
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
+ </g>
+ <g id="text_3">
+ <!-- $\psi = 0$ -->
+ <g transform="translate(294.555318 19.125175)scale(0.168 -0.168)">
+ <defs>
+ <path d="M 64.296875 31.515625
+C 65.796875 31.515625 67.296875 31.515625 67.296875 33.203125
+C 67.296875 35 65.59375 35 63.90625 35
+L 8 35
+C 6.296875 35 4.59375 35 4.59375 33.203125
+C 4.59375 31.515625 6.09375 31.515625 7.59375 31.515625
+z
+M 63.90625 14
+C 65.59375 14 67.296875 14 67.296875 15.796875
+C 67.296875 17.484375 65.796875 17.484375 64.296875 17.484375
+L 7.59375 17.484375
+C 6.09375 17.484375 4.59375 17.484375 4.59375 15.796875
+C 4.59375 14 6.296875 14 8 14
+z
+" id="CMR17-61"/>
+ </defs>
+ <use transform="scale(0.996264)" xlink:href="#CMMI12-32"/>
+ <use transform="translate(94.667488 0)scale(0.996264)" xlink:href="#CMR17-61"/>
+ <use transform="translate(194.056286 0)scale(0.996264)" xlink:href="#CMR17-48"/>
+ </g>
+ </g>
+ </g>
+ </g>
+ <defs>
+ <clipPath id="p6c1199b75e">
+ <rect height="105.035294" width="105.035294" x="10.7" y="25.125175"/>
+ </clipPath>
+ <clipPath id="p31b09fbec5">
+ <rect height="105.035294" width="105.035294" x="136.742353" y="25.125175"/>
+ </clipPath>
+ <clipPath id="p368b104918">
+ <rect height="105.035294" width="105.035294" x="262.784706" y="25.125175"/>
+ </clipPath>
+ </defs>
+</svg>
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_non_matching_quadrature_generator_h
+#define dealii_non_matching_quadrature_generator_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_restriction.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/std_cxx17/optional.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/non_matching/immersed_surface_quadrature.h>
+
+#include <functional>
+
+DEAL_II_NAMESPACE_OPEN
+namespace NonMatching
+{
+ namespace internal
+ {
+ namespace QuadratureGeneratorImplementation
+ {
+ template <int dim, int spacedim>
+ class QGenerator;
+ }
+ } // namespace internal
+
+
+ /**
+ * Struct storing settings for the QuadratureGenerator class.
+ */
+ struct AdditionalQGeneratorData
+ {
+ /**
+ * Constructor.
+ */
+ AdditionalQGeneratorData(const unsigned int max_box_splits = 4,
+ const double lower_bound_implicit_function = 1e-11,
+ const double min_distance_between_roots = 1e-12,
+ const double limit_to_be_definite = 1e-11,
+ const double root_finder_tolerance = 1e-12,
+ const unsigned int max_root_finder_splits = 2,
+ bool split_in_half = true);
+
+ /**
+ * The number of times we are allowed to split the incoming box
+ * and recurse on each child.
+ */
+ unsigned int max_box_splits;
+
+ /**
+ * For a level set function, $\psi$, the implicit function theorem states
+ * that it is possible to write one of the coordinates $x_i$ as a function
+ * of the others if
+ *
+ * $|\frac{\partial \psi}{\partial x_i}| > 0$.
+ *
+ * In practice, it might happend the bound we have for the expression in
+ * the left-hand side is only floating-point close to zero.
+ *
+ * This constant is a safety margin, $C$, that states that the implicit
+ * function theorem can be used when
+ *
+ * $|\frac{\partial \psi}{\partial x_i}| > C$
+ *
+ * Thus this constant must be non-negative.
+ */
+ double lower_bound_implicit_function;
+
+ /**
+ * If two roots are closer to each other than this distance they are
+ * merged to one.
+ */
+ double min_distance_between_roots;
+
+ /**
+ * A constant, $C$, controlling when a level set function, $\psi$, is
+ * considered positive or negative definite:
+ *
+ * $\psi(x) > C \Rightarrow \text{Positive definite}$,
+ * $\psi(x) < -C \Rightarrow \text{Negative definite}$.
+ */
+ double limit_to_be_definite;
+
+ /**
+ * Tolerance for convergence of the underlying root finder.
+ */
+ double root_finder_tolerance;
+
+ /**
+ * The number of times the underlying rootfinder is allowed to split
+ * an interval, while trying to find multiple roots.
+ */
+ unsigned int max_root_finder_splits;
+
+ /**
+ * This determines how a box is split when this is necessary. If true, the
+ * box is split in two, if set to false the box is split into its $2^{dim}$
+ * children.
+ */
+ bool split_in_half;
+ };
+
+
+
+ /**
+ * This class creates immersed quadrature rules over a BoundingBox,
+ * $B \subset \mathbb{R}^{dim}$, when the domain is described by a level set
+ * function, $\psi$.
+ *
+ * This class creates quadrature rules for the intersections between the box
+ * and the three different regions defined by the level set function. That is,
+ * it creates quadrature rules to integrate over the following regions
+ * @f[
+ * N = \{x \in B : \psi(x) < 0 \}, \\
+ * P = \{x \in B : \psi(x) > 0 \}, \\
+ * S = \{x \in B : \psi(x) = 0 \}.
+ * @f]
+ * @image html immersed_quadratures.svg
+ *
+ * When working with level set functions, the most common is to describe a
+ * domain, $\Omega$, as
+ * @f[
+ * \Omega = \{ x \in \mathbb{R}^{dim} : \psi(x) < 0 \}.
+ * @f]
+ * Given this, we shall use the name convention that $N$ is the "inside"
+ * region (i.e. inside $\Omega$), $P$ is the "outside" region and $S$ is
+ * the "surface" region. The "inside" and "outside" quadratures will also be
+ * referred to as the "bulk"-quadratures.
+ *
+ * The underlying algorithm use a 1-dimensional quadrature rule as base for
+ * creating the immersed quadrature rules. Gauss-Legendre quadrature
+ * (QGauss) is recommended. The constructor takes an hp::QCollection<1>.
+ * One can select which 1D-quadrature in the collection should be used
+ * through the set_1D_quadrature() function. The number of quadrature points
+ * in the constructed quadratures will vary depending on the level set
+ * function. More quadrature points will be created if the intersection is
+ * "bad", for example, if the zero-contour has a high curvature compared to
+ * the size of the box. However, if the number of points in the 1D quadrature
+ * is $n$ the number of points will be proportional to $n^{dim}$ in the bulk
+ * quadratures and to $n^{dim-1}$ in the surface quadrature. For example,
+ * in the 2D-example in the above figure, there are 2 points in the
+ * 1D-quadrature. If the 1D-quadrature is a Gauss-Legendre quadrature and the
+ * grid has size $h$, the immersed quadratures typically give global errors
+ * proportional to $h^{2n}$, both for the bulk and surface integrals. If the
+ * 1D-quadrature has positive weights, the weights of the immersed quadratures
+ * will also be positive.
+ *
+ * A detailed description of the underlying algorithm can be found in
+ * "High-Order %Quadrature Methods for Implicitly Defined Surfaces and
+ * Volumes in Hyperrectangles, R. I. Saye, SIAM J. Sci. Comput., 37(2), <a
+ * href="http://www.dx.doi.org/10.1137/140966290">
+ * doi:10.1137/140966290</a>". This implementation has some modifications
+ * compared to the algorithm description in the paper. In particular, it
+ * builds the three different types of quadratures (inside, outside and
+ * surface) simultaneously. Further, the so-called "pruning" step is not yet
+ * implemented.
+ */
+ template <int dim>
+ class QuadratureGenerator
+ {
+ public:
+ using AdditionalData = AdditionalQGeneratorData;
+
+ /**
+ * Constructor. Each Quadrature<1> in @p quadratures1D can be chosen as base
+ * for generating the immersed quadrature rules.
+ *
+ * @note It is important that each 1D-quadrature rule in the
+ * hp::QCollection does not contain the points 0 and 1.
+ */
+ QuadratureGenerator(
+ const hp::QCollection<1> &quadratures1D,
+ const AdditionalData & additional_data = AdditionalData());
+
+ /**
+ * Construct immersed quadratures rules for the incoming level set
+ * function over the BoundingBox.
+ *
+ * To get the constructed quadratures, use the functions
+ * get_inside_quadrature(),
+ * get_outside_quadrature(),
+ * get_surface_quadrature().
+ *
+ * @note Both value, gradient and hessian need to be implemented on the
+ * incoming function.
+ */
+ void
+ generate(const Function<dim> &level_set, const BoundingBox<dim> &box);
+
+ /**
+ * Return the quadrature rule for the region
+ * $\{x \in B : \psi(x) < 0 \}$
+ * created in the previous call to generate().
+ * Here, $B$ is BoundingBox passed to generate().
+ */
+ const Quadrature<dim> &
+ get_inside_quadrature() const;
+
+ /**
+ * Return the quadrature rule for the region
+ * $\{x \in B : \psi(x) > 0 \}$
+ * created in the previous call to generate().
+ * Here, $B$ is BoundingBox passed to generate().
+ */
+ const Quadrature<dim> &
+ get_outside_quadrature() const;
+
+ /**
+ * Return the quadrature rule for the region
+ * $\{x \in B : \psi(x) = 0 \}$
+ * created in the previous call to generate().
+ * Here, $B$ is BoundingBox passed to generate().
+ *
+ * @note The normal at the quadrature points will be parallel to $\nabla \psi$.
+ */
+ const ImmersedSurfaceQuadrature<dim> &
+ get_surface_quadrature() const;
+
+ /**
+ * Set which 1D-quadrature in the collection passed to the constructor
+ * should be used to create the immersed quadratures.
+ */
+ void
+ set_1D_quadrature(const unsigned int q_index);
+
+ private:
+ /**
+ * QuadratureGenerator is mainly used to start up the recursive
+ * algorithm. This is the object that actually generates the quadratures.
+ */
+ internal::QuadratureGeneratorImplementation::QGenerator<dim, dim>
+ q_generator;
+ };
+
+ namespace internal
+ {
+ namespace QuadratureGeneratorImplementation
+ {
+ /**
+ * A class that attempts to find multiple distinct roots of a function,
+ * $f(x)$, over an interval, $[l, r]$. This is done as follows. If there
+ * is a sign change in function value between the interval end points,
+ * we solve for the root. If there is no sign change, we attempt to
+ * bound the function value away from zero on $[a, b]$, to conclude that
+ * no roots exist. If we can't exclude that there are roots, we split
+ * the interval in two: $[l, (r + l) / 2]$, $[(r + l) / 2, r]$, and use
+ * the same algorithm recursively on each interval. This means that we
+ * can typically find 2 distinct roots, but not 3.
+ *
+ * The bounds on the functions values are estimated using the function
+ * taylor_estimate_function_bounds, which approximates the function as a
+ * second order Taylor-polynomial around the interval midpoint.
+ * When we have a sign change on an interval, this class uses
+ * boost::math::tools::toms748_solve for finding roots .
+ */
+ class RootFinder
+ {
+ public:
+ /**
+ * Struct storing settings for the RootFinder class.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor.
+ */
+ AdditionalData(const double tolerance = 1e-12,
+ const unsigned int max_recursion_depth = 2,
+ const unsigned int max_iterations = 500);
+
+ /**
+ * The tolerance in the stopping criteria for the underlying root
+ * finding algorithm boost::math::tools::toms748_solve.
+ */
+ double tolerance;
+
+ /**
+ * The number of times we are allowed to split the interval where we
+ * seek roots.
+ */
+ unsigned int max_recursion_depth;
+
+ /**
+ * The maximum number of iterations in
+ * boost::math::tools::toms748_solve.
+ */
+ unsigned int max_iterations;
+ };
+
+
+ /**
+ * Constructor.
+ */
+ RootFinder(const AdditionalData &data = AdditionalData());
+
+ /**
+ * For each of the incoming @p functions, attempt to find the roots over
+ * the interval defined by @p interval and add these to @p roots.
+ * The returned roots will be sorted in ascending order:
+ * $x_0 < x_1 <...$ and duplicate roots (with respect to the tolerance
+ * in AdditionalData) will be removed.
+ */
+ void
+ find_roots(const std::vector<std::reference_wrapper<const Function<1>>>
+ & functions,
+ const BoundingBox<1> &interval,
+ std::vector<double> & roots);
+
+ private:
+ /**
+ * Attempt to find the roots of the @p function over the interval defined by
+ * @p interval and add these to @p roots. @p recursion_depth holds the number
+ * of times this function has been called recursively.
+ */
+ void
+ find_roots(const Function<1> & function,
+ const BoundingBox<1> &interval,
+ const unsigned int recursion_depth,
+ std::vector<double> & roots);
+
+ const AdditionalData additional_data;
+ };
+
+
+ /**
+ * This is a special Quadrature class with a push_back() method for
+ * conveniently adding a point with an associated weight.
+ *
+ * Since we build the quadrature rules in step-wise fashion,
+ * it's easier to use this class than to pass around two vectors:
+ * std::vector<Point<dim>>,
+ * std::vector<double>.
+ * Further, two std::vectors could accidentally end up with different
+ * sizes. Using push_back we make sure that the number of points and
+ * weights are the same.
+ */
+ template <int dim>
+ class ExtendableQuadrature : public Quadrature<dim>
+ {
+ public:
+ /**
+ * Constructor, creates an empty quadrature rule with no
+ * points.
+ */
+ ExtendableQuadrature() = default;
+
+ /**
+ * Constructor, copies the incoming Quadrature.
+ */
+ ExtendableQuadrature(const Quadrature<dim> &quadrature);
+
+ /**
+ * Add a point with an associated weight to the quadrature.
+ */
+ void
+ push_back(const Point<dim> &point, const double weight);
+ };
+
+
+ /**
+ * Type that describes the definiteness of a function over a region.
+ */
+ enum class Definiteness
+ {
+ negative,
+ positive,
+ indefinite
+ };
+
+
+ /**
+ * Class that stores quadrature rules to integrate over 4 different
+ * regions of a single BoundingBox, $B$. Given multiple level set
+ * functions,
+ *
+ * $\psi_i : \mathbb{R}^{dim} \rightarrow \mathbb{R}$, $i = 0, 1, ...$,
+ *
+ * the box, $B \subset \mathbb{R}^{dim}$, is partitioned into a
+ * "negative", "positive", and "indefinite" region, $B = N \cup P \cup I$,
+ * according to the signs of $\psi_i$ over each region:
+ *
+ * @f[
+ * N = \{x \in B : \psi_i(x) < 0, \forall i \}, \\
+ * P = \{x \in B : \psi_i(x) > 0, \forall i \}, \\
+ * I = B \setminus (\overline{N} \cup \overline{P}).
+ * @f]
+ *
+ * Thus, all $\psi_i$ are positive over $P$ and negative over $N$. Over
+ * $I$ the level set functions differ in sign. This class holds quadrature
+ * rules for each of these regions. In addition, when there is a single
+ * level set function, $\psi$, it holds a surface quadrature for the zero
+ * contour of $\psi$:
+ *
+ * $S = \{x \in B : \psi(x) = 0 \}$.
+ *
+ * Note that when there is a single level set function, $I$ is empty
+ * and $N$ and $P$ are the regions that one typically integrates over in
+ * an immersed finite element method.
+ */
+ template <int dim>
+ class QPartitioning
+ {
+ public:
+ /**
+ * Return a reference to the "bulk" quadrature with the same name as the
+ * member in Definiteness.
+ */
+ ExtendableQuadrature<dim> &
+ quadrature_by_definiteness(const Definiteness definiteness);
+
+ /**
+ * Quadrature for the region $\{x \in B : \psi_i(x) < 0 \forall i \}$ of
+ * the box, $B$.
+ */
+ ExtendableQuadrature<dim> negative;
+
+ /**
+ * Quadrature for the region $\{x \in B : \psi_i(x) > 0 \forall i \}$ of
+ * the box, $B$.
+ */
+ ExtendableQuadrature<dim> positive;
+
+ /**
+ * Quadrature for a region where the level set functions have different
+ * sign.
+ */
+ ExtendableQuadrature<dim> indefinite;
+
+ /**
+ * Quadrature for the region $\{x \in B : \psi(x) = 0 \}$ of the
+ * box, $B$.
+ */
+ ImmersedSurfaceQuadrature<dim> surface;
+ };
+
+
+ /**
+ * This class is responsible for creating quadrature points for
+ * the $dim$-dimensional quadrature partitioning from an
+ * $(dim - 1)$-dimensional "indefinite" quadrature (see
+ * QPartitioning documentation).
+ *
+ * To be precise, let $[L, R]$ be the extents of the box in the height
+ * function direction and let $I \subset \mathbb{R}^{dim-1}$ be the lower
+ * dimensional indefinite region. This class will create quadrature points
+ * over $I \times [L, R] \subset \mathbb{R}^{dim}$ and in the case
+ * $dim=spacedim$, points for the surface quadrature.
+ *
+ * For each lower dimensional quadrature point, $(x_I, w_I)$ in the
+ * indefinite quadrature, we create several 1D-level set functions by
+ * restricting $\psi_j$ to $x_I$. We then partition the interval $[L, R]$
+ * into $[y_0, y_1, ..., y_n]$, where $y_0 = L$, $y_n = R$, and the
+ * remaining $y_i$ are the roots of the 1D-level set functions in
+ * $[L, R]$. Since the level set functions change sign between the
+ * roots, each interval belong to different regions in the quadrature
+ * partitioning.
+ *
+ * In each interval, $[y_i, y_{i+1}]$, we distribute points
+ * according to the 1D-base quadrature, $(x_q, w_q)$ and take the
+ * cartesian product with $(x_I, w_I)$ to create the $dim$-dimensional
+ * quadrature points, $(X_q, W_q)$:
+ * $X_q = x_I \times (y_i + (y_{i+1} - y_i) x_q)$,
+ * $W_q = w_I (y_{i+1} - y_i) w_q$.
+ *
+ * When $dim=spacedim$, we have a single level set function, $\psi$. Since
+ * we have fulfilled the implicit function theorem, there is a single root
+ * $y_1 \in [L, R]$. The point, $x_s = x_I \times y_1$, will be added as a
+ * point in the surface quadrature. One can show that the correct weight
+ * of this point is
+ *
+ * $w_s = \frac{\|\nabla \psi(x_s)\|}{|\partial_i \psi(x_s)|} w_I$,
+ *
+ * where $i$ is the height function direction.
+ */
+ template <int dim, int spacedim>
+ class UpThroughDimensionCreator
+ {
+ public:
+ /**
+ * Constructor. Takes the same parameters as QuadratureGenerator.
+ */
+ UpThroughDimensionCreator(
+ const hp::QCollection<1> & q_collection1D,
+ const AdditionalQGeneratorData &additional_data);
+
+ /**
+ * Create $dim$-dimensional immersed quadratures from the incoming
+ * $(dim-1)$-dimensional quadratures and add these to
+ * @p q_partitioning.
+ */
+ void
+ generate(const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> & box,
+ const Quadrature<dim - 1> &low_dim_quadrature,
+ const unsigned int height_function_direction,
+ QPartitioning<dim> & q_partitioning);
+
+ /**
+ * Set which 1D-quadrature in the collection passed to the constructor
+ * should be used to create the immersed quadratures.
+ */
+ void
+ set_1D_quadrature(const unsigned int q_index);
+
+ private:
+ /**
+ * Create a surface quadrature point from the lower-dimensional point
+ * and add it to surface_quadrature.
+ *
+ * This function is only called when $dim=spacedim$ and there is a
+ * single level set function. At this point there should only be a
+ * single root in the interval $[L, R]$
+ */
+ void
+ create_surface_point(
+ const Point<dim - 1> &point,
+ const double weight,
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> & box,
+ const unsigned int height_function_direction,
+ ImmersedSurfaceQuadrature<dim> &surface_quadrature);
+
+ /**
+ * One dimensional quadrature rules used to create the immersed
+ * quadratures.
+ */
+ const SmartPointer<const hp::QCollection<1>> q_collection1D;
+
+ /**
+ * Stores options/settings for the algorithm.
+ */
+ const AdditionalQGeneratorData additional_data;
+
+ /**
+ * Which quadrature rule in the above collection that is used to
+ * create the immersed quadrature rules.
+ */
+ unsigned int q_index;
+
+ /**
+ * 1D-functions, that are restrictions of each dim-dimensional level set
+ * function passed to generate() to some $(dim-1)$-dimensional point.
+ */
+ std::vector<Functions::PointRestriction<dim - 1>> point_restrictions;
+
+ /**
+ * Class used to find the roots of the above 1D-restictions.
+ */
+ RootFinder root_finder;
+
+ /**
+ * The roots of the functions in point_restrictions.
+ * This will be the values of the height functions, $\{H_i(x_I)\}$ at
+ * some lower dimensional quadrature point,
+ * $x_I \in \mathbb{R}^{dim-1}}$.
+ */
+ std::vector<double> roots;
+ };
+
+
+ /**
+ * Data representing the best choice of height-function direction,
+ * which is returned by the function find_best_height_direction.
+ *
+ * This data consists of a coordinate direction
+ *
+ * $i \in \{0, ..., dim - 1 \}$,
+ *
+ * and lower bound on the absolute value of the derivative of some
+ * associated function, f, taken in the above coordinate direction. That
+ * is, a bound $C$ such that
+ *
+ * $|\frac{\partial f}{\partial x_i}| > C$,
+ *
+ * holding over some subset of $\mathbb{R}^{dim}$.
+ */
+ struct HeightDirectionData
+ {
+ /**
+ * Constructor. Initializes the direction to invalid_unsigned_int and
+ * the bound to 0.
+ */
+ HeightDirectionData();
+
+
+ /**
+ * The height-function direction, described above.
+ */
+ unsigned int direction;
+
+ /**
+ * The lower bound on $|\frac{\partial f}{\partial x_i}|$, described
+ * above.
+ */
+ double min_abs_dfdx;
+ };
+
+
+ /**
+ * Base class for the class QGenerator<dim, spacedim> and the
+ * one-dimensional specialization QGenerator<1, spacedim>.
+ */
+ template <int dim, int spacedim>
+ class QGeneratorBase
+ {
+ public:
+ QGeneratorBase(const hp::QCollection<1> & q_collection1D,
+ const AdditionalQGeneratorData &additional_data);
+
+ /**
+ * Clear the quadratures created by the previous call to generate().
+ */
+ void
+ clear_quadratures();
+
+ /**
+ * Return the created quadratures.
+ */
+ const QPartitioning<dim> &
+ get_quadratures() const;
+
+ protected:
+ /**
+ * Stores options/settings for the algorithm.
+ */
+ const AdditionalQGeneratorData additional_data;
+
+ /**
+ * Which 1D-quadrature in the collection we should use to generate
+ * the immersed quadrature.
+ */
+ unsigned int q_index;
+
+ /**
+ * Index of the quadrature in q_collection1D that should use to
+ * generate the immersed quadrature rules.
+ */
+ const SmartPointer<const hp::QCollection<1>> q_collection1D;
+
+ /**
+ * Quadratures that the derived classes create.
+ */
+ QPartitioning<dim> q_partitioning;
+ };
+
+
+ /**
+ * This class implements the Saye-algorithm cited in the documentation of
+ * the QuadratureGenerator class.
+ *
+ * The generate function takes a number of $dim$-dimensional level set
+ * functions, $\psi_i$, and a BoundingBox<dim>, and builds a partitioning
+ * of quadratures, as defined in documentation of the QPartitioning class.
+ * That is, this class builds an object of type QPartitioning<dim>.
+ *
+ * If all $\psi_i$ passed to generate can be determined to be positive or
+ * negative definite, the QPartitioning will consist of a single
+ * quadrature forming a tensor product.
+ *
+ * If this is not the case, the algorithm uses recursion over the spatial
+ * dimension. The spacedim template parameter denotes the dimension we
+ * started with and dim denotes on what level we are in the recursion.
+ * That is, we first construct a QPartitioning<dim - 1> and then
+ * build the higher dimensional quadratures from these. What we in the end
+ * actually want is a spacedim-dimensional partitioning of quadratures,
+ * for a single level set function, $\psi$.
+ *
+ * The algorithm is based on the implicit function theorem. Starting with
+ * a single level set function, $\psi$, we try to find a direction $i$,
+ * such that
+ *
+ * $|\frac{\partial \psi}{\partial x_i}| > 0$.
+ *
+ * throughout the whole box. This means that the zero-contour of the
+ * level set function can be parameterized by an implicit function
+ *
+ * $H = H(x_0, ..., x_{i-1}, x_{i+1}, ..., x_{dim-1})$,
+ *
+ * so that
+ *
+ * $\psi(..., x_{i-1}, H(..., x_{i-1}, x_{i+1}, ...), x_{i+1}, ...) = 0$,
+ *
+ * over a subset, $I \subset C \subset \mathbb{R}^{dim-1}$, of the cross
+ * section, $C$, of the box (see BoundingBox::cross_section). Here, $I$ is
+ * the "indefinite"-region defined in the QPartitioning class. To follow
+ * convention in the original paper, we will -refer to $H$ as the
+ * "height-function" and to $i$ as the "height-function direction".
+ *
+ * If a height function direction can be found, we go down in dimension by
+ * creating two new level set functions, $\{\psi_0, \psi_1\}$, which are
+ * the restriction of $\psi$ to the top and bottom faces of the box (in
+ * the height function direction). We then delegate to
+ * QGenerator<dim-1, spacedim> to create a QPartitioning<dim-1> over
+ * the cross section.
+ *
+ * When we reach the base case, $dim = 1$, the creation of
+ * QPartitioning<1> is simple. See the documentation in specialized
+ * class: QGenerator<1, spacedim>.
+ *
+ * As we go up through the dimensions and create the higher dimensional
+ * quadratures, we need to know the function value of the height
+ * functions at the lower dimensional quadrature points. Since the
+ * functions are implicit, we need to do root-finding on the level set
+ * functions to find the function values. For this we use the class
+ * UpThroughDimensionCreator, see documentation there.
+ *
+ * When we have $n$ level set functions (i.e. after having gone
+ * down in dimension), we try to find a height function direction,
+ * which works for all those $\psi_i$ which are intersected by the zero
+ * contour (i.e. those not positive or negative definite).
+ * If such a direction exist, we will have a maximum of $n$ associated
+ * implicit height functions, $H_j$. Each $H_j$ parametrize the
+ * $x_i$-coordinate of the zero-contour over a region, $I_j$. The
+ * indefinite region in the lower dimensional partitioning is the union of
+ * these $I = \cup_j I_j$.
+ *
+ * As we try to find a height function direction, we estimate bounds on
+ * the gradient components by approximating each component as a 1st-order
+ * Taylor-polynomial. If a direction can not be found, the box is split
+ * and we recurse on each smaller box. This makes an implicit function
+ * more likely to exist since we seek it over a smaller portion of the
+ * zero contour. It also makes the estimated bounds tighter since we
+ * extrapolate the Taylor-polynomial a shorter distance.
+ *
+ * Since we can not split a box forever, there is an maximum number of
+ * allowed splits on the additional data struct passed to the constructor.
+ * If this is reached, the algorithm uses the midpoint method as a last
+ * resort.
+ */
+ template <int dim, int spacedim>
+ class QGenerator : public QGeneratorBase<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor. Takes the same parameters QuadratureGenerator.
+ */
+ QGenerator(const hp::QCollection<1> & q_collection1D,
+ const AdditionalQGeneratorData &additional_data);
+
+ /**
+ * Create immersed quadrature rules over the incoming @p box and add
+ * these to the internal QPartitioning<dim> object in the base class.
+ * These quadratures can then be obtained using the
+ * get_quadratures-function.
+ *
+ * This function calls itself if the incoming box need to be split.
+ * @p n_box_splits counts the number of times this function has called
+ * itself.
+ */
+ void
+ generate(const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box,
+ const unsigned int n_box_splits);
+
+ /**
+ * Set which 1D-quadrature in the collection passed to the constructor
+ * should be used to create the immersed quadratures.
+ */
+ void
+ set_1D_quadrature(const unsigned int q_index);
+
+ private:
+ /**
+ * Restricts the incoming level set functions to the top and bottom of
+ * the incoming box (w.r.t @p height_function_direction). Then call the
+ * lower dimensional QGenerator with the cross section of the box
+ * to generate the lower dimensional immersed quadrature rules.
+ */
+ void
+ create_low_dim_quadratures(
+ const unsigned int height_function_direction,
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box,
+ const unsigned int n_box_splits);
+
+ /**
+ * Gets the $(dim - 1)$-dimensional quadratures from the lower
+ * dimensional algorithm and creates the $dim$-dimensional quadrature
+ * rules over the box from the lower dimensional ones.
+ */
+ void
+ create_high_dim_quadratures(
+ const unsigned int height_function_direction,
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box);
+
+ /**
+ * Split the incoming box and call generate() recursively with each box.
+ * The box is split in 2 or 4 parts depending on the value of
+ * AdditionalQGeneratorData::split_in_half.
+ */
+ void
+ split_box_and_recurse(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> & box,
+ const std_cxx17::optional<HeightDirectionData> &direction_data,
+ const unsigned int n_box_splits);
+
+ /**
+ * Uses the midpoint-method to create a quadrature over the box.
+ * That is, add a single quadrature point at the center of the box
+ * with weight corresponding to the volume of the box.
+ *
+ * The point is added to the region defined in QPartitioning
+ * according to the signs of the level set functions at the center of
+ * the box.
+ */
+ void
+ use_midpoint_method(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box);
+
+ /**
+ * The same algorithm as this, but creating immersed quadratures
+ * in one dimension lower.
+ */
+ QGenerator<dim - 1, spacedim> low_dim_algorithm;
+
+ /**
+ * Object responsible for creating the $dim$-dimensional quadratures
+ * from
+ */
+ UpThroughDimensionCreator<dim, spacedim> up_through_dimension_creator;
+
+ /**
+ * Stores tensor products of each of the Quadrature<1>'s in
+ * q_collection1D.
+ */
+ hp::QCollection<dim> tensor_products;
+ };
+
+
+ /**
+ * The 1D-base case of the recursive algorithm QGenerator<dim, spacedim>.
+ *
+ * Let $L$ and $R$ be the left and right bounds of the one-dimensional
+ * BoundingBox. This interval is partitioned into $[x_0, x_1, ..., x_n]$
+ * where $x_0 = L$, $x_n = R$, and the remaining $x_i$ are the roots
+ * of the level set functions in the interval $[L, R]$. In each interval,
+ * $[x_i, x_{i+1}]$, quadrature points are distributed according to a
+ * 1D-quadrature rule. These points are added to one of the regions of
+ * QPartitioning determined from the signs of the level set
+ * functions on the interval (see documentation of QPartitioning).
+ *
+ * If spacedim = 1 the points $[x_1, x_{n-1}]$ are also added as surface
+ * quadrature points to QPartitioning::surface.
+ */
+ template <int spacedim>
+ class QGenerator<1, spacedim> : public QGeneratorBase<1, spacedim>
+ {
+ public:
+ /**
+ * Constructor. Takes the same parameters QuadratureGenerator.
+ */
+ QGenerator(const hp::QCollection<1> & quadratures1D,
+ const AdditionalQGeneratorData &additional_data);
+
+ /**
+ * Creates quadrature points over the interval defined by the incoming
+ * box and adds these quadrature points to the internally stored
+ * QPartitioning. These quadratures can then be obtained using
+ * the get_quadratures-function.
+ */
+ void
+ generate(const std::vector<std::reference_wrapper<const Function<1>>>
+ & level_sets,
+ const BoundingBox<1> &box,
+ const unsigned int n_box_splits);
+
+ /**
+ * Set which 1D-quadrature in the collection passed to the constructor
+ * should be used to create the immersed quadratures.
+ */
+ void
+ set_1D_quadrature(const unsigned int q_index);
+
+ private:
+ /**
+ * Adds the point defined by coordinate to the surface quadrature of
+ * ImmersedQuadrature with unit weight.
+ */
+ void
+ create_surface_points(
+ const std::vector<std::reference_wrapper<const Function<1>>>
+ &level_sets);
+
+ /**
+ * Class used to find the roots of the functions passed to generate().
+ */
+ RootFinder root_finder;
+
+ /**
+ * Roots of the functions passed to generate().
+ */
+ std::vector<double> roots;
+
+ /**
+ * This would be the height-function direction in higher dimensions,
+ * but in 1D there is only one coordinate direction.
+ */
+ const unsigned int direction = 0;
+
+ /**
+ * To reuse the distribute_points_between_roots()-function
+ * we need a zero-dimensional quadrature point with unit weight.
+ */
+ const Point<0> zero_dim_point;
+ const double unit_weight = 1;
+ };
+
+
+ /**
+ * Take the tensor product between (point, weight) and @p quadrature1D
+ * scaled over [start, end] and add the resulting dim-dimensional
+ * quadrature points to @p quadrature.
+ *
+ * @p component_in_dim specifies which dim-dimensional coordinate
+ * quadrature1D should be written to.
+ */
+ template <int dim>
+ void
+ tensor_point_with_1D_quadrature(const Point<dim - 1> &point,
+ const double weight,
+ const Quadrature<1> & quadrature1D,
+ const double start,
+ const double end,
+ const unsigned int component_in_dim,
+ ExtendableQuadrature<dim> &quadrature);
+
+
+ /**
+ * Checks the sign of the incoming Functions at the incoming point and
+ * returns Definiteness::positive/negative if all the functions are
+ * positive/negative at the point, otherwise returns
+ * Definiteness::indefinite.
+ */
+ template <int dim>
+ Definiteness
+ pointwise_definiteness(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & functions,
+ const Point<dim> &point);
+
+
+ /**
+ * A struct storing the bounds on the function value and bounds
+ * on each component of the gradient.
+ */
+ template <int dim>
+ struct FunctionBounds
+ {
+ public:
+ /**
+ * Lower and upper bounds on the functions value.
+ */
+ std::pair<double, double> value;
+
+ /**
+ * Lower and upper bounds on each component of the gradient.
+ */
+ std::array<std::pair<double, double>, dim> gradient;
+ };
+
+
+ /**
+ * Returns the max/min bounds on the value, taken over all the entries
+ * in the incoming vector of FunctionBounds. That is, given the incoming
+ * function bounds, $[L_j, U_j]$, this function returns
+ * $[L, U]$,
+ * where $L = \min_{j} L_j$ and $U = \max_{j} U_j$.
+ */
+ template <int dim>
+ std::pair<double, double>
+ find_extreme_values(
+ const std::vector<FunctionBounds<dim>> &all_function_bounds);
+
+
+ /**
+ * Finds the best choice of height function direction, given the
+ * FunctionBounds for a number of functions $\{\psi_j\}_{j=0}^{n-1}$.
+ * Here, "best" is meant in the sense of the implicit function theorem.
+ *
+ * Let $J_I$ be the index set of the indefinite functions:
+ *
+ * $J_I = \{0,..., n - 1\} \setminus \{ j : |\psi_j| > 0 \}$.
+ *
+ * This function converts the incoming bounds to a lower bound, $L_{ij}$,
+ * on the absolute value of each component of the gradient:
+ *
+ * $|\partial_k \psi_j| > L_{jk}$.
+ *
+ * and then returns a coordindate direction, $i$, and a lower bound $L$,
+ * such that
+ *
+ * @f[
+ * i = \arg \max_{k} \min_{j \in J_I} L_{jk}, \\
+ * L = \max_{k} \min_{j \in J_I} L_{jk}.
+ * @f]
+ *
+ * This means $i$ is a coordinate direction such that all functions
+ * intersected by the zero contour (i.e. those belonging to $J_I$) fulfill
+ *
+ * $|\partial_i \psi_j| > L$.
+ *
+ * Note that the estimated lower bound, $L$, can be zero or negative. This
+ * means that no suitable height function direction exists. If all of the
+ * incoming functions are positive or negative definite the returned
+ * std::optional is non-set.
+ */
+ template <int dim>
+ std_cxx17::optional<HeightDirectionData>
+ find_best_height_direction(
+ const std::vector<FunctionBounds<dim>> &all_function_bounds);
+
+ } // namespace QuadratureGeneratorImplementation
+ } // namespace internal
+
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
+
+#endif /* dealii_non_matching_quadrature_generator_h */
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
SET(_src
+ quadrature_generator.cc
coupling.cc
immersed_surface_quadrature.cc
)
SET(_inst
+ quadrature_generator.inst.in
coupling.inst.in
)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/function_tools.h>
+
+#include <deal.II/grid/reference_cell.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <boost/math/special_functions/relative_difference.hpp>
+#include <boost/math/special_functions/sign.hpp>
+#include <boost/math/tools/roots.hpp>
+
+#include <algorithm>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+namespace NonMatching
+{
+ namespace internal
+ {
+ namespace QuadratureGeneratorImplementation
+ {
+ template <int dim>
+ void
+ tensor_point_with_1D_quadrature(const Point<dim - 1> &point,
+ const double weight,
+ const Quadrature<1> & quadrature1D,
+ const double start,
+ const double end,
+ const unsigned int component_in_dim,
+ ExtendableQuadrature<dim> &quadrature)
+ {
+ Assert(start < end,
+ ExcMessage("Interval start must be less than interval end."));
+
+ const double length = end - start;
+ for (unsigned int j = 0; j < quadrature1D.size(); ++j)
+ {
+ const double x = start + length * quadrature1D.point(j)[0];
+ quadrature.push_back(dealii::internal::create_higher_dim_point(
+ point, component_in_dim, x),
+ length * weight * quadrature1D.weight(j));
+ }
+ }
+
+
+
+ /**
+ * For each (point, weight) in lower create a dim-dimensional quadrature
+ * using tensor_point_with_1D_quadrature and add the results to @p quadrature.
+ */
+ template <int dim>
+ void
+ add_tensor_product(const Quadrature<dim - 1> &lower,
+ const Quadrature<1> & quadrature1D,
+ const double start,
+ const double end,
+ const unsigned int component_in_dim,
+ ExtendableQuadrature<dim> &quadrature)
+ {
+ for (unsigned int j = 0; j < lower.size(); ++j)
+ {
+ tensor_point_with_1D_quadrature(lower.point(j),
+ lower.weight(j),
+ quadrature1D,
+ start,
+ end,
+ component_in_dim,
+ quadrature);
+ }
+ }
+
+
+
+ template <int dim>
+ Definiteness
+ pointwise_definiteness(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & functions,
+ const Point<dim> &point)
+ {
+ Assert(functions.size() > 0,
+ ExcMessage(
+ "The incoming vector must contain at least one function."));
+
+ const int sign_of_first =
+ boost::math::sign(functions[0].get().value(point));
+
+ if (sign_of_first == 0)
+ return Definiteness::indefinite;
+
+ for (unsigned int j = 1; j < functions.size(); ++j)
+ {
+ const int sign = boost::math::sign(functions[j].get().value(point));
+
+ if (sign != sign_of_first)
+ return Definiteness::indefinite;
+ }
+ // If we got here all functions have the same sign.
+ if (sign_of_first < 0)
+ return Definiteness::negative;
+ else
+ return Definiteness::positive;
+ }
+
+
+
+ /**
+ * Given the incoming lower and upper bounds on the value of a function
+ * $[L, U]$, return the minimum/maximum of $[L, U]$ and the function
+ * values at the vertices. That is, this function returns
+ *
+ * $[\min(L, L_f), \max(U, U_f)]$,
+ *
+ * where $L_f = \min_{v} f(x_v)$, $U_f = \max_{v} f(x_v)|$,
+ * and $x_v$ is a vertex.
+ *
+ * It is assumed that the incoming function is scalar valued.
+ */
+ template <int dim>
+ void
+ take_min_max_at_vertices(const Function<dim> & function,
+ const BoundingBox<dim> & box,
+ std::pair<double, double> &value_bounds)
+ {
+ const ReferenceCell &cube = ReferenceCells::get_hypercube<dim>();
+ for (unsigned int i = 0; i < cube.n_vertices(); ++i)
+ {
+ const double vertex_value = function.value(box.vertex(i));
+
+ value_bounds.first = std::min(value_bounds.first, vertex_value);
+ value_bounds.second = std::max(value_bounds.second, vertex_value);
+ }
+ }
+
+
+
+ /**
+ * Estimate bounds on each of the functions in the incoming vector over
+ * the incoming box.
+ *
+ * Bounds on the functions value and the gradient components are first
+ * computed using FunctionTools::taylor_estimate_function_bounds.
+ * In addition, the function value is checked for min/max at the at
+ * the vertices of the box. The gradient is not checked at the box
+ * vertices.
+ */
+ template <int dim>
+ void
+ estimate_function_bounds(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & functions,
+ const BoundingBox<dim> & box,
+ std::vector<FunctionBounds<dim>> &all_function_bounds)
+ {
+ all_function_bounds.clear();
+ all_function_bounds.reserve(functions.size());
+ for (const Function<dim> &function : functions)
+ {
+ FunctionBounds<dim> bounds;
+ FunctionTools::taylor_estimate_function_bounds<dim>(
+ function, box, bounds.value, bounds.gradient);
+ take_min_max_at_vertices(function, box, bounds.value);
+
+ all_function_bounds.push_back(bounds);
+ }
+ }
+
+
+
+ template <int dim>
+ std::pair<double, double>
+ find_extreme_values(const std::vector<FunctionBounds<dim>> &bounds)
+ {
+ Assert(bounds.size() > 0, ExcMessage("The incoming vector is empty."));
+
+ std::pair<double, double> extremes = bounds[0].value;
+ for (unsigned int i = 1; i < bounds.size(); ++i)
+ {
+ extremes.first = std::min(extremes.first, bounds[i].value.first);
+ extremes.second = std::max(extremes.second, bounds[i].value.second);
+ }
+
+ return extremes;
+ }
+
+
+
+ /**
+ * Return true if the incoming function bounds correspond to a function
+ * which is indefinite, i.e., that is not negative or positive definite.
+ */
+ inline bool
+ is_indefinite(const std::pair<double, double> &function_bounds)
+ {
+ if (function_bounds.first > 0)
+ return false;
+ if (function_bounds.second < 0)
+ return false;
+ return true;
+ }
+
+
+
+ /**
+ * Return a lower bound, $L_a$, on the absolute value of a function,
+ * $f(x)$:
+ *
+ * $L_a \leq |f(x)|$,
+ *
+ * by estimating it from the incoming lower and upper bounds:
+ * $L \leq f(x) \leq U$.
+ *
+ * By rewriting the lower and upper bounds as
+ * $F - C \leq f(x) \leq F + C$,
+ * where $L = F - C$, $U = F + C$ (or $F = (U + L)/2$, $C = (U - L)/2$),
+ * we get $|f(x) - F| \leq C$.
+ * Using the inverse triangle inequality gives
+ * $|F| - |f(x)| \leq |f(x) - F| \leq C$.
+ * Thus, $L_a = |F| - C$.
+ *
+ * Note that the returned value can be negative. This is used to indicate
+ * "how far away" a function is from being definite.
+ */
+ inline double
+ lower_bound_on_abs(const std::pair<double, double> &function_bounds)
+ {
+ Assert(function_bounds.first <= function_bounds.second,
+ ExcMessage("Function bounds reversed, max < min."));
+
+ return 0.5 * (std::abs(function_bounds.second + function_bounds.first) -
+ (function_bounds.second - function_bounds.first));
+ }
+
+
+
+ HeightDirectionData::HeightDirectionData()
+ {
+ direction = numbers::invalid_unsigned_int;
+ min_abs_dfdx = 0;
+ }
+
+
+
+ template <int dim>
+ std_cxx17::optional<HeightDirectionData>
+ find_best_height_direction(
+ const std::vector<FunctionBounds<dim>> &all_function_bounds)
+ {
+ // Minimum (taken over the indefinite functions) on the lower bound on
+ // each component of the gradient.
+ std_cxx17::optional<std::array<double, dim>> min_lower_abs_grad;
+
+ for (const FunctionBounds<dim> &bounds : all_function_bounds)
+ {
+ if (is_indefinite(bounds.value))
+ {
+ // For the first indefinite function we find, we write the lower
+ // bounds on each gradient component to min_lower_abs_grad.
+ if (!min_lower_abs_grad)
+ {
+ min_lower_abs_grad.emplace();
+ for (int d = 0; d < dim; ++d)
+ {
+ (*min_lower_abs_grad)[d] =
+ lower_bound_on_abs(bounds.gradient[d]);
+ }
+ }
+ else
+ {
+ for (int d = 0; d < dim; ++d)
+ {
+ (*min_lower_abs_grad)[d] =
+ std::min((*min_lower_abs_grad)[d],
+ lower_bound_on_abs(bounds.gradient[d]));
+ }
+ }
+ }
+ }
+
+ if (min_lower_abs_grad)
+ {
+ const auto max_element =
+ std::max_element(min_lower_abs_grad->begin(),
+ min_lower_abs_grad->end());
+
+ HeightDirectionData data;
+ data.direction =
+ std::distance(min_lower_abs_grad->begin(), max_element);
+ data.min_abs_dfdx = *max_element;
+
+ return data;
+ }
+
+ return std_cxx17::optional<HeightDirectionData>();
+ }
+
+
+
+ /**
+ * Return true if there are exactly two incoming FunctionBounds and
+ * they corresponds to one function being positive definite and
+ * one being negative definite. Return false otherwise.
+ */
+ template <int dim>
+ inline bool
+ one_positive_one_negative_definite(
+ const std::vector<FunctionBounds<dim>> &all_function_bounds)
+ {
+ if (all_function_bounds.size() != 2)
+ return false;
+ else
+ {
+ const FunctionBounds<dim> &bounds0 = all_function_bounds.at(0);
+ const FunctionBounds<dim> &bounds1 = all_function_bounds.at(1);
+
+ if (bounds0.value.first > 0 && bounds1.value.second < 0)
+ return true;
+ if (bounds1.value.first > 0 && bounds0.value.second < 0)
+ return true;
+ return false;
+ }
+ }
+
+
+
+ /**
+ * Transform the points and weights of the incoming quadrature,
+ * unit_quadrature, from unit space to the incoming box and add these to
+ * quadrature.
+ *
+ * Note that unit_quadrature should be a quadrature over [0,1]^dim.
+ */
+ template <int dim>
+ void
+ map_quadrature_to_box(const Quadrature<dim> & unit_quadrature,
+ const BoundingBox<dim> & box,
+ ExtendableQuadrature<dim> &quadrature)
+ {
+ for (unsigned int i = 0; i < unit_quadrature.size(); i++)
+ {
+ const Point<dim> point = box.unit_to_real(unit_quadrature.point(i));
+ const double weight = unit_quadrature.weight(i) * box.volume();
+
+ quadrature.push_back(point, weight);
+ }
+ }
+
+
+
+ /**
+ * For each of the incoming dim-dimensional functions, create the
+ * restriction to the top and bottom of the incoming BoundingBox and add
+ * these two (dim-1)-dimensional functions to @p restrictions. Here, top and bottom is
+ * meant with respect to the incoming @p direction. For each function, the
+ * "bottom-restriction" will be added before the "top-restriction"
+ *
+ * @note @p restrictions will be cleared, so after this function
+ * restrictions.size() == 2 * functions.size().
+ */
+ template <int dim>
+ void
+ restrict_to_top_and_bottom(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & functions,
+ const BoundingBox<dim> & box,
+ const unsigned int direction,
+ std::vector<Functions::CoordinateRestriction<dim - 1>> &restrictions)
+ {
+ AssertIndexRange(direction, dim);
+
+ restrictions.clear();
+ restrictions.reserve(2 * functions.size());
+
+ const double bottom = box.lower_bound(direction);
+ const double top = box.upper_bound(direction);
+
+ for (const auto &function : functions)
+ {
+ restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
+ function, direction, bottom));
+ restrictions.push_back(Functions::CoordinateRestriction<dim - 1>(
+ function, direction, top));
+ }
+ }
+
+
+
+ /**
+ * Restrict each of the incoming @p functions to @p point,
+ * while keeping the coordinate direction @p open_direction open,
+ * and add the restriction to @p restrictions.
+ *
+ * @note @p restrictions will be cleared, so after this function
+ * restrictions.size() == functions.size().
+ */
+ template <int dim>
+ void
+ restrict_to_point(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & functions,
+ const Point<dim - 1> & point,
+ const unsigned int open_direction,
+ std::vector<Functions::PointRestriction<dim - 1>> &restrictions)
+ {
+ AssertIndexRange(open_direction, dim);
+
+ restrictions.clear();
+ restrictions.reserve(functions.size());
+ for (const auto &function : functions)
+ {
+ restrictions.push_back(Functions::PointRestriction<dim - 1>(
+ function, open_direction, point));
+ }
+ }
+
+
+
+ /**
+ * Let $\{ y_0, ..., y_{n+1} \}$ be such that $[y_0, y_{n+1}]$ is the
+ * @p interval and $\{ y_1, ..., y_n \}$ are the @p roots. In each
+ * subinterval, $[y_i, y_{i+1}]$, distribute point according to the
+ * 1D-quadrature rule $\{(x_q, w_q)\}_q$ (@p quadrature1D).
+ * Take the tensor product with the quadrature point $(x, w)$
+ * (@p point, @p weight) to create dim-dimensional quadrature points
+ * @f[
+ * X_q = x_I \times (y_i + (y_{i+1} - y_i) x_q),
+ * W_q = w_I (y_{i+1} - y_i) w_q,
+ * @f]
+ * and add these points to @p q_partitioning.
+ */
+ template <int dim>
+ void
+ distribute_points_between_roots(
+ const Quadrature<1> & quadrature1D,
+ const BoundingBox<1> & interval,
+ const std::vector<double> &roots,
+ const Point<dim - 1> & point,
+ const double weight,
+ const unsigned int height_function_direction,
+ const std::vector<std::reference_wrapper<const Function<1>>>
+ & level_sets,
+ const AdditionalQGeneratorData &additional_data,
+ QPartitioning<dim> & q_partitioning)
+ {
+ // Make this int to avoid a warning signed/unsigned comparision.
+ const int n_roots = roots.size();
+
+ // The number of intervals are roots.size() + 1
+ for (int i = -1; i < n_roots; ++i)
+ {
+ // Start and end point of the subinterval.
+ const double start = i < 0 ? interval.lower_bound(0) : roots[i];
+ const double end =
+ i + 1 < n_roots ? roots[i + 1] : interval.upper_bound(0);
+
+ const double length = end - start;
+ // It might be that the end points of the subinterval are roots.
+ // If this is the case then the subinterval has length zero.
+ // Don't distribute points on the subinterval if it is shorter than
+ // some tolerance.
+ if (length > additional_data.min_distance_between_roots)
+ {
+ // All points on the interval belong to the same region in
+ // the QPartitioning. Determine the quadrature we should add
+ // the points to.
+ const Point<1> center(start + 0.5 * length);
+ const Definiteness definiteness =
+ pointwise_definiteness(level_sets, center);
+ ExtendableQuadrature<dim> &target_quadrature =
+ q_partitioning.quadrature_by_definiteness(definiteness);
+
+ tensor_point_with_1D_quadrature(point,
+ weight,
+ quadrature1D,
+ start,
+ end,
+ height_function_direction,
+ target_quadrature);
+ }
+ }
+ }
+
+
+
+ RootFinder::AdditionalData::AdditionalData(
+ const double tolerance,
+ const unsigned int max_recursion_depth,
+ const unsigned int max_iterations)
+ : tolerance(tolerance)
+ , max_recursion_depth(max_recursion_depth)
+ , max_iterations(max_iterations)
+ {}
+
+
+
+ RootFinder::RootFinder(const AdditionalData &data)
+ : additional_data(data)
+ {}
+
+
+
+ void
+ RootFinder::find_roots(
+ const std::vector<std::reference_wrapper<const Function<1>>> &functions,
+ const BoundingBox<1> & interval,
+ std::vector<double> & roots)
+ {
+ for (const Function<1> &function : functions)
+ {
+ const unsigned int recursion_depth = 0;
+ find_roots(function, interval, recursion_depth, roots);
+ }
+ // Sort and make sure no roots are duplicated
+ std::sort(roots.begin(), roots.end());
+
+ const auto roots_are_equal = [this](const double &a, const double &b) {
+ return std::abs(a - b) < additional_data.tolerance;
+ };
+ roots.erase(unique(roots.begin(), roots.end(), roots_are_equal),
+ roots.end());
+ }
+
+
+
+ void
+ RootFinder::find_roots(const Function<1> & function,
+ const BoundingBox<1> &interval,
+ const unsigned int recursion_depth,
+ std::vector<double> & roots)
+ {
+ // Compute function values at end points.
+ const double left_value = function.value(interval.vertex(0));
+ const double right_value = function.value(interval.vertex(1));
+
+ // If we have a sign change we solve for the root.
+ if (boost::math::sign(left_value) != boost::math::sign(right_value))
+ {
+ const auto lambda = [&function](const double x) {
+ return function.value(Point<1>(x));
+ };
+
+ const auto stopping_criteria = [this](const double &a,
+ const double &b) {
+ return std::abs(a - b) < additional_data.tolerance;
+ };
+
+ boost::uintmax_t iterations = additional_data.max_iterations;
+
+ const std::pair<double, double> root_bracket =
+ boost::math::tools::toms748_solve(lambda,
+ interval.lower_bound(0),
+ interval.upper_bound(0),
+ left_value,
+ right_value,
+ stopping_criteria,
+ iterations);
+
+ const double root = .5 * (root_bracket.first + root_bracket.second);
+ roots.push_back(root);
+ }
+ else
+ {
+ // Compute bounds on the incoming function to check if there are
+ // roots. If the function is positive or negative on the whole
+ // interval we do nothing.
+ std::pair<double, double> value_bounds;
+ std::array<std::pair<double, double>, 1> gradient_bounds;
+ FunctionTools::taylor_estimate_function_bounds<1>(function,
+ interval,
+ value_bounds,
+ gradient_bounds);
+
+ // Since we already know the function values at the interval ends we
+ // might as well check these for min/max too.
+ const double function_min =
+ std::min(std::min(left_value, right_value), value_bounds.first);
+
+ // If the functions is positive there are no roots.
+ if (function_min > 0)
+ return;
+
+ const double function_max =
+ std::max(std::max(left_value, right_value), value_bounds.second);
+
+ // If the functions is negative there are no roots.
+ if (function_max < 0)
+ return;
+
+ // If we can't say that the function is strictly positive/negative
+ // we split the interval in half. We can't split forever, so if we
+ // have reached the max recursion, we stop looking for roots.
+ if (recursion_depth < additional_data.max_recursion_depth)
+ {
+ find_roots(function,
+ interval.child(0),
+ recursion_depth + 1,
+ roots);
+ find_roots(function,
+ interval.child(1),
+ recursion_depth + 1,
+ roots);
+ }
+ }
+ }
+
+
+
+ template <int dim>
+ ExtendableQuadrature<dim>::ExtendableQuadrature(
+ const Quadrature<dim> &quadrature)
+ : Quadrature<dim>(quadrature)
+ {}
+
+
+
+ template <int dim>
+ void
+ ExtendableQuadrature<dim>::push_back(const Point<dim> &point,
+ const double weight)
+ {
+ this->quadrature_points.push_back(point);
+ this->weights.push_back(weight);
+ }
+
+
+
+ template <int dim>
+ ExtendableQuadrature<dim> &
+ QPartitioning<dim>::quadrature_by_definiteness(
+ const Definiteness definiteness)
+ {
+ switch (definiteness)
+ {
+ case Definiteness::negative:
+ return negative;
+ case Definiteness::positive:
+ return positive;
+ default:
+ return indefinite;
+ }
+ }
+
+
+
+ /**
+ * Takes a (dim-1)-dimensional point from the cross-section (orthogonal
+ * to direction) of the box. Creates the two dim-dimensional points, which
+ * are the projections from the cross section to the faces of the box and
+ * returns the point closest to the zero-contour of the incoming level set
+ * function.
+ */
+ template <int dim>
+ Point<dim>
+ face_projection_closest_zero_contour(const Point<dim - 1> & point,
+ const unsigned int direction,
+ const BoundingBox<dim> &box,
+ const Function<dim> & level_set)
+ {
+ const Point<dim> bottom_point =
+ dealii::internal::create_higher_dim_point(point,
+ direction,
+ box.lower_bound(direction));
+ const double bottom_value = level_set.value(bottom_point);
+
+ const Point<dim> top_point =
+ dealii::internal::create_higher_dim_point(point,
+ direction,
+ box.upper_bound(direction));
+ const double top_value = level_set.value(top_point);
+
+ // The end point closest to the zero-contour is the one with smallest
+ // absolute value.
+ return std::abs(bottom_value) < std::abs(top_value) ? bottom_point :
+ top_point;
+ }
+
+
+
+ template <int dim, int spacedim>
+ UpThroughDimensionCreator<dim, spacedim>::UpThroughDimensionCreator(
+ const hp::QCollection<1> & q_collection1D,
+ const AdditionalQGeneratorData &additional_data)
+ : q_collection1D(&q_collection1D)
+ , additional_data(additional_data)
+ , root_finder(
+ RootFinder::AdditionalData(additional_data.root_finder_tolerance,
+ additional_data.max_root_finder_splits))
+ {
+ q_index = 0;
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ UpThroughDimensionCreator<dim, spacedim>::generate(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> & box,
+ const Quadrature<dim - 1> &low_dim_quadrature,
+ const unsigned int height_function_direction,
+ QPartitioning<dim> & q_partitioning)
+ {
+ const Quadrature<1> &quadrature1D = (*q_collection1D)[q_index];
+
+ for (unsigned int q = 0; q < low_dim_quadrature.size(); ++q)
+ {
+ const Point<dim - 1> &point = low_dim_quadrature.point(q);
+ const double weight = low_dim_quadrature.weight(q);
+ restrict_to_point(level_sets,
+ point,
+ height_function_direction,
+ point_restrictions);
+
+ // We need a vector of references to do the recursive call.
+ const std::vector<std::reference_wrapper<const Function<1>>>
+ restrictions(point_restrictions.begin(),
+ point_restrictions.end());
+
+ const BoundingBox<1> bounds_in_direction =
+ box.bounds(height_function_direction);
+
+ roots.clear();
+ root_finder.find_roots(restrictions, bounds_in_direction, roots);
+
+ distribute_points_between_roots(quadrature1D,
+ bounds_in_direction,
+ roots,
+ point,
+ weight,
+ height_function_direction,
+ restrictions,
+ additional_data,
+ q_partitioning);
+
+ if (dim == spacedim)
+ create_surface_point(point,
+ weight,
+ level_sets,
+ box,
+ height_function_direction,
+ q_partitioning.surface);
+ }
+
+ point_restrictions.clear();
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ UpThroughDimensionCreator<dim, spacedim>::create_surface_point(
+ const Point<dim - 1> &point,
+ const double weight,
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> & box,
+ const unsigned int height_function_direction,
+ ImmersedSurfaceQuadrature<dim> &surface_quadrature)
+ {
+ AssertIndexRange(roots.size(), 2);
+ Assert(level_sets.size() == 1, ExcInternalError());
+
+
+ const Function<dim> &level_set = level_sets.at(0);
+
+ Point<dim> surface_point;
+ if (roots.size() == 1)
+ {
+ surface_point = dealii::internal::create_higher_dim_point(
+ point, height_function_direction, roots[0]);
+ }
+ else
+ {
+ // If we got here, we have missed roots in the lower dimensional
+ // algorithm. This is a rare event but can happen if the
+ // zero-contour has a high curvature. The problem is that the
+ // incoming point has been incorrectly added to the indefinite
+ // quadrature in QPartitioning<dim-1>. Since we missed a root on
+ // this box, we will likely miss it on the neighboring box too. If
+ // this happens, the point will NOT be in the indefinite quadrature
+ // on the neighbor. The best thing we can do is to compute the
+ // surface point by projecting the lower dimensional point to the
+ // face closest to the zero-contour. We should add a surface point
+ // because the neighbor will not.
+ surface_point = face_projection_closest_zero_contour(
+ point, height_function_direction, box, level_set);
+ }
+
+ const Tensor<1, dim> gradient = level_set.gradient(surface_point);
+ Tensor<1, dim> normal = gradient;
+ normal *= 1. / normal.norm();
+
+ // Note that gradient[height_function_direction] is non-zero
+ // because of the implicit function theorem.
+ const double surface_weight =
+ weight * gradient.norm() /
+ std::abs(gradient[height_function_direction]);
+ surface_quadrature.push_back(surface_point, surface_weight, normal);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ UpThroughDimensionCreator<dim, spacedim>::set_1D_quadrature(
+ unsigned int q_index)
+ {
+ AssertIndexRange(q_index, q_collection1D->size());
+ this->q_index = q_index;
+ }
+
+
+
+ template <int dim, int spacedim>
+ QGeneratorBase<dim, spacedim>::QGeneratorBase(
+ const hp::QCollection<1> & q_collection1D,
+ const AdditionalQGeneratorData &additional_data)
+ : additional_data(additional_data)
+ , q_collection1D(&q_collection1D)
+ {
+ q_index = 0;
+ }
+
+
+
+ template <int dim, int spacedim>
+ QGenerator<dim, spacedim>::QGenerator(
+ const hp::QCollection<1> & q_collection1D,
+ const AdditionalQGeneratorData &additional_data)
+ : QGeneratorBase<dim, spacedim>(q_collection1D, additional_data)
+ , low_dim_algorithm(q_collection1D, additional_data)
+ , up_through_dimension_creator(q_collection1D, additional_data)
+ {
+ for (unsigned int i = 0; i < q_collection1D.size(); i++)
+ tensor_products.push_back(Quadrature<dim>(q_collection1D[i]));
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ QGeneratorBase<dim, spacedim>::clear_quadratures()
+ {
+ q_partitioning = QPartitioning<dim>();
+ }
+
+
+
+ template <int dim, int spacedim>
+ const QPartitioning<dim> &
+ QGeneratorBase<dim, spacedim>::get_quadratures() const
+ {
+ return q_partitioning;
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ QGenerator<dim, spacedim>::generate(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box,
+ const unsigned int n_box_splits)
+ {
+ std::vector<FunctionBounds<dim>> all_function_bounds;
+ estimate_function_bounds(level_sets, box, all_function_bounds);
+
+ const std::pair<double, double> extreme_values =
+ find_extreme_values(all_function_bounds);
+
+ if (extreme_values.first > this->additional_data.limit_to_be_definite)
+ {
+ map_quadrature_to_box(tensor_products[this->q_index],
+ box,
+ this->q_partitioning.positive);
+ }
+ else if (extreme_values.second <
+ -(this->additional_data.limit_to_be_definite))
+ {
+ map_quadrature_to_box(tensor_products[this->q_index],
+ box,
+ this->q_partitioning.negative);
+ }
+ else if (one_positive_one_negative_definite(all_function_bounds))
+ {
+ map_quadrature_to_box(tensor_products[this->q_index],
+ box,
+ this->q_partitioning.indefinite);
+ }
+ else
+ {
+ const std_cxx17::optional<HeightDirectionData> data =
+ find_best_height_direction(all_function_bounds);
+
+ // Check larger than a constant to avoid that min_abs_dfdx is only
+ // larger by 0 by floating point precision.
+ if (data && data->min_abs_dfdx >
+ this->additional_data.lower_bound_implicit_function)
+ {
+ create_low_dim_quadratures(data->direction,
+ level_sets,
+ box,
+ n_box_splits);
+ create_high_dim_quadratures(data->direction, level_sets, box);
+ }
+ else if (n_box_splits < this->additional_data.max_box_splits)
+ {
+ split_box_and_recurse(level_sets, box, data, n_box_splits);
+ }
+ else
+ {
+ // We can't split the box recursively forever. Use the midpoint
+ // method as a last resort.
+ use_midpoint_method(level_sets, box);
+ }
+ }
+ }
+
+
+
+ /**
+ * Return the coordinate direction of the largest side of the box.
+ * If two or more sides have the same length the returned std::optional
+ * will be non-set.
+ */
+ template <int dim>
+ std_cxx17::optional<unsigned int>
+ direction_of_largest_extent(const BoundingBox<dim> &box)
+ {
+ // Get the side lengths for each direction and sort them.
+ std::array<std::pair<double, unsigned int>, dim> side_lengths;
+ for (int i = 0; i < dim; i++)
+ {
+ side_lengths[i].first = box.side_length(i);
+ side_lengths[i].second = i;
+ }
+ // Sort is lexicographic, so this sorts based on side length first.
+ std::sort(side_lengths.begin(), side_lengths.end());
+
+ // Check if the two largest side lengths have the same length. This
+ // function isn't called in 1D, so the (dim - 2)-element exists.
+ if (boost::math::epsilon_difference(side_lengths[dim - 1].first,
+ side_lengths[dim - 2].first) < 100)
+ return std_cxx17::optional<unsigned int>();
+
+ return side_lengths.back().second;
+ }
+
+
+
+ /**
+ * Return the coordinate direction that the box should be split in,
+ * assuming that the box should be split it half.
+ *
+ * If the box is larger in one coordante direction, this direction is
+ * returned. If the box have the same extent in all directions, we choose
+ * the coordinate direction which is closest to being a height-function
+ * direction. That is, the direction $i$ that has a least negative
+ * estimate of $|\partial_i \psi_j|$. As a last resort, we choose the
+ * direction 0, if @p height_direction_data non-set.
+ */
+ template <int dim>
+ unsigned int
+ compute_split_direction(
+ const BoundingBox<dim> & box,
+ const std_cxx17::optional<HeightDirectionData> &height_direction_data)
+ {
+ const std_cxx17::optional<unsigned int> direction =
+ direction_of_largest_extent(box);
+
+ if (direction)
+ return *direction;
+
+ // This direction is closest to being a height direction, so
+ // we split in this direction.
+ if (height_direction_data)
+ return height_direction_data->direction;
+
+ // We have to choose some direction, we might aswell take 0.
+ return 0;
+ }
+
+
+
+ /**
+ * Split the incoming box in half with respect to the incoming coordinate
+ * direction and return the left half.
+ */
+ template <int dim>
+ inline BoundingBox<dim>
+ left_half(const BoundingBox<dim> &box, const unsigned int direction)
+ {
+ AssertIndexRange(direction, dim);
+
+ // Move the upper corner half a side-length to the left.
+ std::pair<Point<dim>, Point<dim>> corners = box.get_boundary_points();
+ corners.second[direction] -= .5 * box.side_length(direction);
+
+ return BoundingBox<dim>(corners);
+ }
+
+
+
+ /**
+ * Split the incoming box in half with respect to the incoming coordinate
+ * direction and return the right half.
+ */
+ template <int dim>
+ inline BoundingBox<dim>
+ right_half(const BoundingBox<dim> &box, const unsigned int direction)
+ {
+ AssertIndexRange(direction, dim);
+
+ // Move the lower corner half a side-length to the right.
+ std::pair<Point<dim>, Point<dim>> corners = box.get_boundary_points();
+ corners.first[direction] += .5 * box.side_length(direction);
+
+ return BoundingBox<dim>(corners);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ QGenerator<dim, spacedim>::split_box_and_recurse(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> & box,
+ const std_cxx17::optional<HeightDirectionData> &direction_data,
+ const unsigned int n_box_splits)
+ {
+ if (this->additional_data.split_in_half)
+ {
+ const unsigned int direction =
+ compute_split_direction(box, direction_data);
+
+ const BoundingBox<dim> left_box = left_half(box, direction);
+ const BoundingBox<dim> right_box = right_half(box, direction);
+
+ generate(level_sets, left_box, n_box_splits + 1);
+ generate(level_sets, right_box, n_box_splits + 1);
+ }
+ else
+ {
+ for (unsigned int i = 0;
+ i < GeometryInfo<dim>::max_children_per_cell;
+ ++i)
+ {
+ generate(level_sets, box.child(i), n_box_splits + 1);
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ QGenerator<dim, spacedim>::create_low_dim_quadratures(
+ const unsigned int height_function_direction,
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box,
+ const unsigned int n_box_splits)
+ {
+ std::vector<Functions::CoordinateRestriction<dim - 1>>
+ face_restrictions;
+ restrict_to_top_and_bottom(level_sets,
+ box,
+ height_function_direction,
+ face_restrictions);
+
+ // We need a vector of references to do the recursive call.
+ const std::vector<std::reference_wrapper<const Function<dim - 1>>>
+ restrictions(face_restrictions.begin(), face_restrictions.end());
+
+ const BoundingBox<dim - 1> cross_section =
+ box.cross_section(height_function_direction);
+
+ low_dim_algorithm.clear_quadratures();
+ low_dim_algorithm.generate(restrictions, cross_section, n_box_splits);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ QGenerator<dim, spacedim>::create_high_dim_quadratures(
+ const unsigned int height_function_direction,
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box)
+ {
+ const QPartitioning<dim - 1> &low_dim_quadratures =
+ low_dim_algorithm.get_quadratures();
+
+ const Quadrature<1> &quadrature1D =
+ (*this->q_collection1D)[this->q_index];
+
+ add_tensor_product(low_dim_quadratures.negative,
+ quadrature1D,
+ box.lower_bound(height_function_direction),
+ box.upper_bound(height_function_direction),
+ height_function_direction,
+ this->q_partitioning.negative);
+
+ add_tensor_product(low_dim_quadratures.positive,
+ quadrature1D,
+ box.lower_bound(height_function_direction),
+ box.upper_bound(height_function_direction),
+ height_function_direction,
+ this->q_partitioning.positive);
+
+ up_through_dimension_creator.generate(level_sets,
+ box,
+ low_dim_quadratures.indefinite,
+ height_function_direction,
+ this->q_partitioning);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ QGenerator<dim, spacedim>::use_midpoint_method(
+ const std::vector<std::reference_wrapper<const Function<dim>>>
+ & level_sets,
+ const BoundingBox<dim> &box)
+ {
+ const Point<dim> center = box.center();
+ const Definiteness definiteness =
+ pointwise_definiteness(level_sets, center);
+
+ ExtendableQuadrature<dim> &quadrature =
+ this->q_partitioning.quadrature_by_definiteness(definiteness);
+
+ quadrature.push_back(center, box.volume());
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ QGenerator<dim, spacedim>::set_1D_quadrature(const unsigned int q_index)
+ {
+ AssertIndexRange(q_index, this->q_collection1D->size());
+
+ this->q_index = q_index;
+ low_dim_algorithm.set_1D_quadrature(q_index);
+ up_through_dimension_creator.set_1D_quadrature(q_index);
+ }
+
+
+
+ template <int spacedim>
+ QGenerator<1, spacedim>::QGenerator(
+ const hp::QCollection<1> & q_collection1D,
+ const AdditionalQGeneratorData &additional_data)
+ : QGeneratorBase<1, spacedim>(q_collection1D, additional_data)
+ , root_finder(
+ RootFinder::AdditionalData(additional_data.root_finder_tolerance,
+ additional_data.max_root_finder_splits))
+ {
+ Assert(q_collection1D.size() > 0,
+ ExcMessage("Incoming quadrature collection is empty."));
+ }
+
+
+
+ template <int spacedim>
+ void
+ QGenerator<1, spacedim>::generate(
+ const std::vector<std::reference_wrapper<const Function<1>>>
+ & level_sets,
+ const BoundingBox<1> &box,
+ const unsigned int n_box_splits)
+ {
+ (void)n_box_splits;
+
+ roots.clear();
+ root_finder.find_roots(level_sets, box, roots);
+
+ const Quadrature<1> &quadrature1D =
+ (*this->q_collection1D)[this->q_index];
+
+ distribute_points_between_roots(quadrature1D,
+ box,
+ roots,
+ zero_dim_point,
+ unit_weight,
+ direction,
+ level_sets,
+ this->additional_data,
+ this->q_partitioning);
+
+ if (spacedim == 1)
+ this->create_surface_points(level_sets);
+ }
+
+
+
+ template <int spacedim>
+ void
+ QGenerator<1, spacedim>::create_surface_points(
+ const std::vector<std::reference_wrapper<const Function<1>>>
+ &level_sets)
+ {
+ Assert(level_sets.size() == 1, ExcInternalError());
+
+ for (const double root : roots)
+ {
+ // A surface integral in 1D is just a point evaluation,
+ // so the weight is always 1.
+ const double weight = 1;
+ const Point<1> point(root);
+
+ Tensor<1, 1> normal = level_sets[0].get().gradient(point);
+ const double gradient_norm = normal.norm();
+ Assert(
+ gradient_norm > 1e-11,
+ ExcMessage(
+ "The level set function has a gradient almost equal to 0."));
+ normal *= 1. / gradient_norm;
+
+ this->q_partitioning.surface.push_back(point, weight, normal);
+ }
+ }
+
+
+
+ template <int spacedim>
+ void
+ QGenerator<1, spacedim>::set_1D_quadrature(const unsigned int q_index)
+ {
+ AssertIndexRange(q_index, this->q_collection1D->size());
+ this->q_index = q_index;
+ }
+ } // namespace QuadratureGeneratorImplementation
+ } // namespace internal
+
+ using namespace internal::QuadratureGeneratorImplementation;
+
+
+
+ AdditionalQGeneratorData::AdditionalQGeneratorData(
+ const unsigned int max_box_splits,
+ const double lower_bound_implicit_function,
+ const double min_distance_between_roots,
+ const double limit_to_be_definite,
+ const double root_finder_tolerance,
+ const unsigned int max_root_finder_splits,
+ bool split_in_half)
+ : max_box_splits(max_box_splits)
+ , lower_bound_implicit_function(lower_bound_implicit_function)
+ , min_distance_between_roots(min_distance_between_roots)
+ , limit_to_be_definite(limit_to_be_definite)
+ , root_finder_tolerance(root_finder_tolerance)
+ , max_root_finder_splits(max_root_finder_splits)
+ , split_in_half(split_in_half)
+ {}
+
+
+
+ template <int dim>
+ QuadratureGenerator<dim>::QuadratureGenerator(
+ const hp::QCollection<1> &q_collection,
+ const AdditionalData & additional_data)
+ : q_generator(q_collection, additional_data)
+ {
+ Assert(q_collection.size() > 0,
+ ExcMessage("Incoming hp::QCollection<1> is empty."));
+ }
+
+
+
+ template <int dim>
+ void
+ QuadratureGenerator<dim>::generate(const Function<dim> & level_set,
+ const BoundingBox<dim> &box)
+ {
+ Assert(level_set.n_components == 1,
+ ExcMessage(
+ "The incoming function should be a scalar level set function,"
+ " it should have one component."));
+ Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume."));
+
+ q_generator.clear_quadratures();
+
+ std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
+ level_sets.push_back(level_set);
+
+ const unsigned int n_box_splits = 0;
+ q_generator.generate(level_sets, box, n_box_splits);
+
+ // With a single level set function, the "indefinite" quadrature should be
+ // zero. If you call generate() with a ZeroFunction nothing good can be
+ // done. You will end up here.
+ Assert(
+ q_generator.get_quadratures().indefinite.size() == 0,
+ ExcMessage(
+ "Generation of quadrature rules failed. This can mean that the level"
+ "set function is degenerate in some way, e.g. oscillating extremely"
+ "rapidly."));
+ }
+
+
+
+ template <int dim>
+ const Quadrature<dim> &
+ QuadratureGenerator<dim>::get_inside_quadrature() const
+ {
+ return q_generator.get_quadratures().negative;
+ }
+
+
+
+ template <int dim>
+ const Quadrature<dim> &
+ QuadratureGenerator<dim>::get_outside_quadrature() const
+ {
+ return q_generator.get_quadratures().positive;
+ }
+
+
+
+ template <int dim>
+ const ImmersedSurfaceQuadrature<dim> &
+ QuadratureGenerator<dim>::get_surface_quadrature() const
+ {
+ return q_generator.get_quadratures().surface;
+ }
+
+
+ template <int dim>
+ void
+ QuadratureGenerator<dim>::set_1D_quadrature(const unsigned int q_index)
+ {
+ q_generator.set_1D_quadrature(q_index);
+ }
+
+} // namespace NonMatching
+#include "quadrature_generator.inst"
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS)
+ {
+ namespace NonMatching
+ \{
+ template class QuadratureGenerator<deal_II_dimension>;
+
+ namespace internal
+ \{
+ namespace QuadratureGeneratorImplementation
+ \{
+ template struct FunctionBounds<deal_II_dimension>;
+
+ template std::pair<double, double>
+ find_extreme_values(
+ const std::vector<FunctionBounds<deal_II_dimension>> &);
+
+ template void
+ estimate_function_bounds(
+ const std::vector<
+ std::reference_wrapper<const Function<deal_II_dimension>>> &,
+ const BoundingBox<deal_II_dimension> &,
+ std::vector<FunctionBounds<deal_II_dimension>> &);
+
+// gcc gives a maybe-uninitialized warning in this function when dim = 1, but
+// gcc is wrong. We don't need the function when dim = 1, so we avoid
+// instantiating it.
+#if 1 < deal_II_dimension
+ template std_cxx17::optional<HeightDirectionData>
+ find_best_height_direction(
+ const std::vector<FunctionBounds<deal_II_dimension>> &);
+#endif
+
+ template void
+ map_quadrature_to_box(const Quadrature<deal_II_dimension> &,
+ const BoundingBox<deal_II_dimension> &,
+ ExtendableQuadrature<deal_II_dimension> &);
+
+ template void
+ tensor_point_with_1D_quadrature(
+ const Point<deal_II_dimension - 1> &,
+ const double,
+ const Quadrature<1> &,
+ const double,
+ const double,
+ const unsigned int,
+ ExtendableQuadrature<deal_II_dimension> &q);
+ \}
+ \}
+ \}
+ }
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
+ {
+#if 1 < deal_II_dimension
+ template class NonMatching::internal::QuadratureGeneratorImplementation::
+ UpThroughDimensionCreator<deal_II_dimension, deal_II_space_dimension>;
+#endif
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test the function find_best_height_direction in
+ * NonMatching::internal::QuadratureGeneratorImplementation.
+ */
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+using namespace NonMatching::internal::QuadratureGeneratorImplementation;
+
+
+// Return a pair with both entries equal to value.
+std::pair<double, double>
+pair_with_equal_entries(const double value)
+{
+ return std::pair<double, double>(value, value);
+}
+
+
+
+/*
+ * Test that find_best_height_direction returns an unset optional if the
+ * incoming bounds correspond to negative/positive definite functions.
+ */
+void
+test_ignores_definite_functions()
+{
+ const int dim = 2;
+ deallog << "test_ignores_definite_functions" << std::endl;
+
+ // Bounds corresponding to one negative and one positive definite function.
+ std::vector<FunctionBounds<dim>> bounds(2);
+ bounds[0].value = pair_with_equal_entries(-1);
+ bounds[1].value = pair_with_equal_entries(1);
+
+ const std_cxx17::optional<HeightDirectionData> data =
+ find_best_height_direction(bounds);
+
+ if (!data)
+ deallog << "OK" << std::endl;
+}
+
+
+
+/**
+ * Create a vector containing two FunctionBounds, set them up so that
+ * there is one height function direction that is the best. Test that this is
+ * the direction returned from find_best_height_direction().
+ */
+void
+test_find_best_height_direction()
+{
+ deallog << "test_find_best_height_direction" << std::endl;
+
+ const int dim = 2;
+
+ std::vector<FunctionBounds<dim>> bounds(2);
+ // Set up so that the bounds correspond to indefinite functions.
+ for (unsigned int i = 0; i < bounds.size(); i++)
+ {
+ bounds[i].value.first = -1;
+ bounds[i].value.second = 1;
+ }
+
+ // Set up the bounds so that the componenetwise min (over function bounds)
+ // of the gradient is [3, 5]. This makes 1 the best direction.
+ bounds[0].gradient[0] = pair_with_equal_entries(3);
+ bounds[1].gradient[0] = pair_with_equal_entries(4);
+ bounds[0].gradient[1] = pair_with_equal_entries(6);
+ bounds[1].gradient[1] = pair_with_equal_entries(5);
+
+ const std_cxx17::optional<HeightDirectionData> data =
+ find_best_height_direction(bounds);
+
+ deallog << "height direction = " << data->direction << std::endl;
+ deallog << "min_abs_dfdx = " << data->min_abs_dfdx << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+ test_ignores_definite_functions();
+ deallog << std::endl;
+ test_find_best_height_direction();
+}
--- /dev/null
+
+DEAL::test_ignores_definite_functions
+DEAL::OK
+DEAL::
+DEAL::test_find_best_height_direction
+DEAL::height direction = 1
+DEAL::min_abs_dfdx = 5.00000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test the function find_extreme_values in
+ * NonMatching::internal::QuadratureGeneratorImplementation.
+ */
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+using namespace NonMatching::internal::QuadratureGeneratorImplementation;
+
+
+/**
+ * Send in a vector with two different function bounds to
+ * find_extreme_values(). Check that what we get back are the actual extreme
+ * values.
+ */
+template <int dim>
+void
+test_extreme_values_are_found()
+{
+ deallog << "test_extreme_values_are_found" << std::endl;
+
+ std::vector<FunctionBounds<dim>> bounds(2);
+ bounds[0].value.first = 1;
+ bounds[0].value.second = 2;
+ bounds[1].value.first = -1;
+ bounds[1].value.second = 3;
+
+ const std::pair<double, double> extremes = find_extreme_values(bounds);
+
+ deallog << "min = " << extremes.first << std::endl;
+ deallog << "max = " << extremes.second << std::endl;
+}
+
+
+
+/**
+ * Since the implementation of find_extreme_values() treats the 0th entry
+ * differently, we check that we get the same entry back if we send in a
+ * vector with only one entry.
+ */
+template <int dim>
+void
+test_extreme_values_initialized_to_first()
+{
+ deallog << "test_extreme_values_initialized_to_first" << std::endl;
+
+ std::vector<FunctionBounds<dim>> bounds(1);
+ bounds[0].value.first = 1;
+ bounds[0].value.second = 2;
+
+ const std::pair<double, double> extremes = find_extreme_values(bounds);
+
+ deallog << "min = " << extremes.first << std::endl;
+ deallog << "max = " << extremes.second << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+ test_extreme_values_initialized_to_first<1>();
+ deallog << std::endl;
+ test_extreme_values_are_found<1>();
+}
--- /dev/null
+
+DEAL::test_extreme_values_initialized_to_first
+DEAL::min = 1.00000
+DEAL::max = 2.00000
+DEAL::
+DEAL::test_extreme_values_are_found
+DEAL::min = -1.00000
+DEAL::max = 3.00000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test the function pointwise_definiteness in
+ * NonMatching::internal::QuadratureGeneratorImplementation.
+ */
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/point.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+using namespace dealii;
+using namespace NonMatching::internal::QuadratureGeneratorImplementation;
+
+
+/**
+ * Call pointwise_definiteness with two positive Functions,
+ * check that it returns Definiteness::positive.
+ */
+template <int dim>
+void
+test_with_positive_functions()
+{
+ std::vector<Functions::ConstantFunction<dim>> functions;
+ functions.push_back(Functions::ConstantFunction<dim>(1));
+ functions.push_back(Functions::ConstantFunction<dim>(1));
+
+ const std::vector<std::reference_wrapper<const Function<dim>>> function_refs(
+ functions.begin(), functions.end());
+
+ const Definiteness definiteness =
+ pointwise_definiteness(function_refs, Point<dim>());
+
+ AssertThrow(definiteness == Definiteness::positive, ExcInternalError());
+}
+
+
+
+/**
+ * Call pointwise_definiteness with two negative Functions,
+ * check that it returns Definiteness::negative.
+ */
+template <int dim>
+void
+test_with_negative_functions()
+{
+ std::vector<Functions::ConstantFunction<dim>> functions;
+ functions.push_back(Functions::ConstantFunction<dim>(-1));
+ functions.push_back(Functions::ConstantFunction<dim>(-1));
+
+ const std::vector<std::reference_wrapper<const Function<dim>>> function_refs(
+ functions.begin(), functions.end());
+
+ const Definiteness definiteness =
+ pointwise_definiteness(function_refs, Point<dim>());
+
+ AssertThrow(definiteness == Definiteness::negative, ExcInternalError());
+}
+
+
+
+/**
+ * Call pointwise_definiteness with with one positive and one negative Function,
+ * check that it returns Definiteness::indefinite.
+ */
+template <int dim>
+void
+test_with_functions_of_different_sign()
+{
+ std::vector<Functions::ConstantFunction<dim>> functions;
+ functions.push_back(Functions::ConstantFunction<dim>(-1));
+ functions.push_back(Functions::ConstantFunction<dim>(1));
+
+ const std::vector<std::reference_wrapper<const Function<dim>>> function_refs(
+ functions.begin(), functions.end());
+
+ const Definiteness definiteness =
+ pointwise_definiteness(function_refs, Point<dim>());
+
+ AssertThrow(definiteness == Definiteness::indefinite, ExcInternalError());
+}
+
+
+
+/**
+ * Call pointwise_definiteness with a single Function which is zero,
+ * check that it returns Definiteness::indefinite.
+ *
+ * This is a special case in the implementation.
+ */
+template <int dim>
+void
+test_first_function_zero()
+{
+ Functions::ZeroFunction<dim> zero_function;
+
+ std::vector<std::reference_wrapper<const Function<dim>>> function_refs;
+ function_refs.push_back(zero_function);
+
+ const Definiteness definiteness =
+ pointwise_definiteness(function_refs, Point<dim>());
+
+ AssertThrow(definiteness == Definiteness::indefinite, ExcInternalError());
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+ test_with_positive_functions<dim>();
+ test_with_negative_functions<dim>();
+ test_with_functions_of_different_sign<dim>();
+ test_first_function_zero<dim>();
+}
+
+
+
+int
+main()
+{
+ initlog();
+ run_test<1>();
+ run_test<2>();
+ run_test<3>();
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test the QuadratureGenerator class, by setting up a few simple cuts over the
+ * unit box and writing the generated quadrature rules to the output file.
+ *
+ * Each function beginning with "test_" sets up a level set function and then
+ * calls the function create_and_print_quadratures() to generate the
+ * quadratures.
+ */
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_level_set.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+#include "quadrature_printing.h"
+
+using NonMatching::QuadratureGenerator;
+
+/*
+ * Create immersed quadrature rules over a unit box intersected by the
+ * incoming level set function. Use a 1D-Gauss quadrature with n_1D_points
+ * points as a base. Print the constructed quadrature rules to deallog.
+ */
+template <int dim>
+void
+create_and_print_quadratures(
+ const Function<dim> & level_set,
+ const unsigned int n_1D_points = 2,
+ const typename QuadratureGenerator<dim>::AdditionalData &additional_data =
+ typename QuadratureGenerator<dim>::AdditionalData())
+{
+ deallog << "dim=" << dim << std::endl;
+
+ hp::QCollection<1> q_collection;
+ q_collection.push_back(QGauss<1>(n_1D_points));
+
+ QuadratureGenerator<dim> quadrature_generator(q_collection, additional_data);
+
+ const BoundingBox<dim> box = create_unit_bounding_box<dim>();
+ quadrature_generator.generate(level_set, box);
+
+ deallog << "Inside quadrature" << std::endl;
+ print_quadrature(quadrature_generator.get_inside_quadrature());
+ deallog << "Outside quadrature" << std::endl;
+ print_quadrature(quadrature_generator.get_outside_quadrature());
+ deallog << "Surface quadrature" << std::endl;
+ print_surface_quadrature(quadrature_generator.get_surface_quadrature());
+}
+
+
+
+/*
+ * Construct level set with a zero contour as a plane cutting straight
+ * through the unit box. Create and print the constructed quadratures. Do this
+ * for all unit normals aligned with the coordinate directions. We expect that
+ * the constructed quadrature has equal number of points in the inside/outside
+ * region and that they are tensor products.
+ */
+template <int dim>
+void
+test_vertical_cuts_through_center()
+{
+ deallog << "test_vertical_cuts_through_center" << std::endl;
+
+ Point<dim> center;
+ for (int d = 0; d < dim; ++d)
+ {
+ center(d) = .5;
+ }
+ for (int direction = 0; direction < dim; ++direction)
+ {
+ deallog << "direction=" << direction << std::endl;
+ const Tensor<1, dim> normal = Point<dim>::unit_vector(direction);
+ const Functions::LevelSet::Plane<dim> level_set(center, normal);
+ create_and_print_quadratures(level_set);
+ }
+}
+
+
+
+/*
+ * Set up a constant negative/positive level set function. Check that
+ * the constructed quadratures only the inside/outside have points and that
+ * this quadrature is a tensor product.
+ */
+template <int dim>
+void
+test_constant_level_sets_both_signs()
+{
+ const Functions::ConstantFunction<dim> constant_positive(1);
+ const Functions::ConstantFunction<dim> constant_negative(-1);
+ deallog << std::endl;
+
+ deallog << "constant_positive" << std::endl;
+ create_and_print_quadratures(constant_positive);
+
+ deallog << std::endl;
+
+ deallog << "constant_negative" << std::endl;
+ create_and_print_quadratures(constant_negative);
+}
+
+
+
+// Set up a level set function corresponding to a plane with normal (1,1) in 2D
+// and (1,1,1) in 3D. This makes the inside region a simplex, with vertices
+// (0, 0), (0, l), (l, 0), in 2D,
+// (0, 0 ,0), (0, 0, l), (0, l, 0), (l, 0, 0), in 3D.
+// where l is the edge length.
+// This is a good test because we know that the inside weights should sum to the
+// area/volume: $V = l^{dim}/dim!$, and that the the surface weights should sum
+// to $S = \sqrt(2) l$ in 2D and $S = \sqrt(3)/2 l^2$ in 3D.
+template <int dim>
+void
+test_simplex_cut()
+{
+ deallog << "test_simplex_cut" << std::endl;
+
+ const double edge_length = 1. / std::sqrt(2);
+
+ Tensor<1, dim> normal;
+ for (int i = 0; i < dim; ++i)
+ normal[i] = 1;
+
+ Point<dim> point_in_plane;
+ point_in_plane[0] = edge_length;
+
+ const Functions::LevelSet::Plane<dim> level_set(point_in_plane, normal);
+
+ create_and_print_quadratures(level_set);
+}
+
+
+
+// Set up a level set function with a zero contour being a plane in the
+// direction (1,1) in 2D and (1,1,1) in 3D, such that it cuts the bottom corner
+// of the reference cell with a cut of size epsilon. Test that the epsilon cut
+// is ignored and we get a tensor product quadrature over the outside region.
+template <int dim>
+void
+test_epsilon_cut_at_bottom_corner()
+{
+ deallog << "test_epsilon_cut_at_bottom_corner" << std::endl;
+ const double epsilon = 1e-15;
+ Tensor<1, dim> normal;
+ Point<dim> center;
+ for (int i = 0; i < dim; ++i)
+ {
+ normal[i] = 1;
+ center(i) += epsilon;
+ }
+ const Functions::LevelSet::Plane<dim> level_set(center, normal);
+
+ create_and_print_quadratures(level_set);
+}
+
+
+
+/*
+ * Set up a spherical level set with radius R centered in (0, R) in 2D and
+ * (0, 0, R) in 3D. The result of this is that the zero contour of the level set
+ * function cuts exactly through vertex 0 of the unit box
+ *
+ * When we choose to split the cell in 4. This test case is difficult for the
+ * algorithm in 3D. We first get dim - 1 as the first height-direction. But
+ * after restricting once we get L_a = 0 (where L_a is defined by |\partial_i
+ * psi| > L_a) for both i = 1,2. Thus we can not choose a second height function
+ * direction. The results is that the cell is split several times until the
+ * maximum recursion is reached. When this happens the algorithm uses the
+ * midpoint method as a fallback.
+ */
+template <int dim>
+void
+test_sphere_cutting_corner_exactly()
+{
+ deallog << "test_sphere_cutting_corner_exactly" << std::endl;
+ const double radius = 4;
+ Point<dim> center;
+ center[dim - 1] = radius;
+ const Functions::LevelSet::Sphere<dim> level_set(center, radius);
+
+ typename QuadratureGenerator<dim>::AdditionalData data;
+ data.split_in_half = false;
+ data.max_box_splits = 2;
+
+ const unsigned int n_1D_points = 2;
+
+ create_and_print_quadratures(level_set, n_1D_points, data);
+}
+
+
+
+// A "fake" function used in test_splitting(). This function is constant 1,
+// except close to the unit box center, x_i = 0.5, where it has a very large
+// Hessian.
+template <int dim>
+class ConstantOneButLargeHessianInCenter
+ : public Functions::ConstantFunction<dim>
+{
+public:
+ ConstantOneButLargeHessianInCenter()
+ : Functions::ConstantFunction<dim>(1)
+ {
+ for (int d = 0; d < dim; ++d)
+ unit_box_center(d) = .5;
+ }
+
+ SymmetricTensor<2, dim>
+ hessian(const Point<dim> &point, const unsigned int) const override
+ {
+ SymmetricTensor<2, dim> hessian;
+
+ const double max_distance = 1e-3;
+ const double diagonal_value =
+ point.distance(unit_box_center) < max_distance ? 1E3 : 0;
+
+ for (int d = 0; d < dim; ++d)
+ hessian[d][d] = diagonal_value;
+
+ return hessian;
+ }
+
+private:
+ Point<dim> unit_box_center;
+};
+
+
+
+// Test the box splitting. Call QuadratureGenerator with a function that is
+// constant 1, but has a large Hessian close to the center of the box. This
+// should make the algorithm split the box, since the function bounds will be
+// large.
+template <int dim>
+void
+test_splitting()
+{
+ deallog << "test_splitting" << std::endl;
+
+ const ConstantOneButLargeHessianInCenter<dim> level_set;
+ create_and_print_quadratures(level_set);
+}
+
+
+
+// Some of the test cases only make sense for a given dimension,
+// so we list the cases for each dimension.
+int
+main()
+{
+ initlog();
+ // 1D
+ test_vertical_cuts_through_center<1>();
+ deallog << std::endl;
+ // 2D
+ test_vertical_cuts_through_center<2>();
+ deallog << std::endl;
+ test_constant_level_sets_both_signs<2>();
+ deallog << std::endl;
+ test_simplex_cut<2>();
+ deallog << std::endl;
+ test_epsilon_cut_at_bottom_corner<2>();
+ deallog << std::endl;
+ test_sphere_cutting_corner_exactly<2>();
+ deallog << std::endl;
+ test_splitting<2>();
+ deallog << std::endl;
+ // 3D
+ test_vertical_cuts_through_center<3>();
+ deallog << std::endl;
+ test_simplex_cut<3>();
+ deallog << std::endl;
+ test_epsilon_cut_at_bottom_corner<3>();
+ deallog << std::endl;
+ test_sphere_cutting_corner_exactly<3>();
+}
--- /dev/null
+
+DEAL::test_vertical_cuts_through_center
+DEAL::direction=0
+DEAL::dim=1
+DEAL::Inside quadrature
+DEAL::0.105662, 0.250000
+DEAL::0.394338, 0.250000
+DEAL::Outside quadrature
+DEAL::0.605662, 0.250000
+DEAL::0.894338, 0.250000
+DEAL::Surface quadrature
+DEAL::0.500000, 1.00000, 1.00000
+DEAL::
+DEAL::test_vertical_cuts_through_center
+DEAL::direction=0
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::0.105662, 0.211325, 0.125000
+DEAL::0.394338, 0.211325, 0.125000
+DEAL::0.105662, 0.788675, 0.125000
+DEAL::0.394338, 0.788675, 0.125000
+DEAL::Outside quadrature
+DEAL::0.605662, 0.211325, 0.125000
+DEAL::0.894338, 0.211325, 0.125000
+DEAL::0.605662, 0.788675, 0.125000
+DEAL::0.894338, 0.788675, 0.125000
+DEAL::Surface quadrature
+DEAL::0.500000, 0.211325, 0.500000, 1.00000, 0.00000
+DEAL::0.500000, 0.788675, 0.500000, 1.00000, 0.00000
+DEAL::direction=1
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::0.211325, 0.105662, 0.125000
+DEAL::0.211325, 0.394338, 0.125000
+DEAL::0.788675, 0.105662, 0.125000
+DEAL::0.788675, 0.394338, 0.125000
+DEAL::Outside quadrature
+DEAL::0.211325, 0.605662, 0.125000
+DEAL::0.211325, 0.894338, 0.125000
+DEAL::0.788675, 0.605662, 0.125000
+DEAL::0.788675, 0.894338, 0.125000
+DEAL::Surface quadrature
+DEAL::0.211325, 0.500000, 0.500000, 0.00000, 1.00000
+DEAL::0.788675, 0.500000, 0.500000, 0.00000, 1.00000
+DEAL::
+DEAL::
+DEAL::constant_positive
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::Outside quadrature
+DEAL::0.211325, 0.211325, 0.250000
+DEAL::0.788675, 0.211325, 0.250000
+DEAL::0.211325, 0.788675, 0.250000
+DEAL::0.788675, 0.788675, 0.250000
+DEAL::Surface quadrature
+DEAL::
+DEAL::constant_negative
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::0.211325, 0.211325, 0.250000
+DEAL::0.788675, 0.211325, 0.250000
+DEAL::0.211325, 0.788675, 0.250000
+DEAL::0.788675, 0.788675, 0.250000
+DEAL::Outside quadrature
+DEAL::Surface quadrature
+DEAL::
+DEAL::test_simplex_cut
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::0.117851, 0.149429, 0.0985844
+DEAL::0.439826, 0.149429, 0.0985844
+DEAL::0.0315781, 0.557678, 0.0264156
+DEAL::0.117851, 0.557678, 0.0264156
+DEAL::Outside quadrature
+DEAL::0.211325, 0.769002, 0.0732233
+DEAL::0.788675, 0.769002, 0.0732233
+DEAL::0.211325, 0.938104, 0.0732233
+DEAL::0.788675, 0.938104, 0.0732233
+DEAL::0.651151, 0.149429, 0.0781923
+DEAL::0.906526, 0.149429, 0.0781923
+DEAL::0.329176, 0.557678, 0.150361
+DEAL::0.820253, 0.557678, 0.150361
+DEAL::Surface quadrature
+DEAL::0.557678, 0.149429, 0.500000, 0.707107, 0.707107
+DEAL::0.149429, 0.557678, 0.500000, 0.707107, 0.707107
+DEAL::
+DEAL::test_epsilon_cut_at_bottom_corner
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::Outside quadrature
+DEAL::0.211325, 0.211325, 0.250000
+DEAL::0.788675, 0.211325, 0.250000
+DEAL::0.211325, 0.788675, 0.250000
+DEAL::0.788675, 0.788675, 0.250000
+DEAL::Surface quadrature
+DEAL::
+DEAL::test_sphere_cutting_corner_exactly
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::0.211325, 0.215731, 0.248603
+DEAL::0.211325, 0.789856, 0.248603
+DEAL::0.788675, 0.273253, 0.230370
+DEAL::0.788675, 0.805269, 0.230370
+DEAL::Outside quadrature
+DEAL::0.211325, 0.00118050, 0.00139654
+DEAL::0.211325, 0.00440568, 0.00139654
+DEAL::0.788675, 0.0165936, 0.0196304
+DEAL::0.788675, 0.0619282, 0.0196304
+DEAL::Surface quadrature
+DEAL::0.211325, 0.00558618, 0.500699, 0.0528312, -0.998603
+DEAL::0.788675, 0.0785218, 0.510012, 0.197169, -0.980370
+DEAL::
+DEAL::test_splitting
+DEAL::dim=2
+DEAL::Inside quadrature
+DEAL::Outside quadrature
+DEAL::0.105662, 0.211325, 0.125000
+DEAL::0.394338, 0.211325, 0.125000
+DEAL::0.105662, 0.788675, 0.125000
+DEAL::0.394338, 0.788675, 0.125000
+DEAL::0.605662, 0.211325, 0.125000
+DEAL::0.894338, 0.211325, 0.125000
+DEAL::0.605662, 0.788675, 0.125000
+DEAL::0.894338, 0.788675, 0.125000
+DEAL::Surface quadrature
+DEAL::
+DEAL::test_vertical_cuts_through_center
+DEAL::direction=0
+DEAL::dim=3
+DEAL::Inside quadrature
+DEAL::0.105662, 0.211325, 0.211325, 0.0625000
+DEAL::0.394338, 0.211325, 0.211325, 0.0625000
+DEAL::0.105662, 0.788675, 0.211325, 0.0625000
+DEAL::0.394338, 0.788675, 0.211325, 0.0625000
+DEAL::0.105662, 0.211325, 0.788675, 0.0625000
+DEAL::0.394338, 0.211325, 0.788675, 0.0625000
+DEAL::0.105662, 0.788675, 0.788675, 0.0625000
+DEAL::0.394338, 0.788675, 0.788675, 0.0625000
+DEAL::Outside quadrature
+DEAL::0.605662, 0.211325, 0.211325, 0.0625000
+DEAL::0.894338, 0.211325, 0.211325, 0.0625000
+DEAL::0.605662, 0.788675, 0.211325, 0.0625000
+DEAL::0.894338, 0.788675, 0.211325, 0.0625000
+DEAL::0.605662, 0.211325, 0.788675, 0.0625000
+DEAL::0.894338, 0.211325, 0.788675, 0.0625000
+DEAL::0.605662, 0.788675, 0.788675, 0.0625000
+DEAL::0.894338, 0.788675, 0.788675, 0.0625000
+DEAL::Surface quadrature
+DEAL::0.500000, 0.211325, 0.211325, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.500000, 0.788675, 0.211325, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.500000, 0.211325, 0.788675, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.500000, 0.788675, 0.788675, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::direction=1
+DEAL::dim=3
+DEAL::Inside quadrature
+DEAL::0.211325, 0.105662, 0.211325, 0.0625000
+DEAL::0.211325, 0.394338, 0.211325, 0.0625000
+DEAL::0.211325, 0.105662, 0.788675, 0.0625000
+DEAL::0.211325, 0.394338, 0.788675, 0.0625000
+DEAL::0.788675, 0.105662, 0.211325, 0.0625000
+DEAL::0.788675, 0.394338, 0.211325, 0.0625000
+DEAL::0.788675, 0.105662, 0.788675, 0.0625000
+DEAL::0.788675, 0.394338, 0.788675, 0.0625000
+DEAL::Outside quadrature
+DEAL::0.211325, 0.605662, 0.211325, 0.0625000
+DEAL::0.211325, 0.894338, 0.211325, 0.0625000
+DEAL::0.211325, 0.605662, 0.788675, 0.0625000
+DEAL::0.211325, 0.894338, 0.788675, 0.0625000
+DEAL::0.788675, 0.605662, 0.211325, 0.0625000
+DEAL::0.788675, 0.894338, 0.211325, 0.0625000
+DEAL::0.788675, 0.605662, 0.788675, 0.0625000
+DEAL::0.788675, 0.894338, 0.788675, 0.0625000
+DEAL::Surface quadrature
+DEAL::0.211325, 0.500000, 0.211325, 0.250000, 0.00000, 1.00000, 0.00000
+DEAL::0.211325, 0.500000, 0.788675, 0.250000, 0.00000, 1.00000, 0.00000
+DEAL::0.788675, 0.500000, 0.211325, 0.250000, 0.00000, 1.00000, 0.00000
+DEAL::0.788675, 0.500000, 0.788675, 0.250000, 0.00000, 1.00000, 0.00000
+DEAL::direction=2
+DEAL::dim=3
+DEAL::Inside quadrature
+DEAL::0.211325, 0.211325, 0.105662, 0.0625000
+DEAL::0.211325, 0.211325, 0.394338, 0.0625000
+DEAL::0.788675, 0.211325, 0.105662, 0.0625000
+DEAL::0.788675, 0.211325, 0.394338, 0.0625000
+DEAL::0.211325, 0.788675, 0.105662, 0.0625000
+DEAL::0.211325, 0.788675, 0.394338, 0.0625000
+DEAL::0.788675, 0.788675, 0.105662, 0.0625000
+DEAL::0.788675, 0.788675, 0.394338, 0.0625000
+DEAL::Outside quadrature
+DEAL::0.211325, 0.211325, 0.605662, 0.0625000
+DEAL::0.211325, 0.211325, 0.894338, 0.0625000
+DEAL::0.788675, 0.211325, 0.605662, 0.0625000
+DEAL::0.788675, 0.211325, 0.894338, 0.0625000
+DEAL::0.211325, 0.788675, 0.605662, 0.0625000
+DEAL::0.211325, 0.788675, 0.894338, 0.0625000
+DEAL::0.788675, 0.788675, 0.605662, 0.0625000
+DEAL::0.788675, 0.788675, 0.894338, 0.0625000
+DEAL::Surface quadrature
+DEAL::0.211325, 0.211325, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000
+DEAL::0.788675, 0.211325, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000
+DEAL::0.211325, 0.788675, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000
+DEAL::0.788675, 0.788675, 0.500000, 0.250000, 0.00000, 0.00000, 1.00000
+DEAL::
+DEAL::test_simplex_cut
+DEAL::dim=3
+DEAL::Inside quadrature
+DEAL::0.0929463, 0.117851, 0.149429, 0.0216800
+DEAL::0.346880, 0.117851, 0.149429, 0.0216800
+DEAL::0.0249049, 0.439826, 0.149429, 0.00580914
+DEAL::0.0929463, 0.439826, 0.149429, 0.00580914
+DEAL::0.0249049, 0.0315781, 0.557678, 0.00155655
+DEAL::0.0929463, 0.0315781, 0.557678, 0.00155655
+DEAL::0.00667324, 0.117851, 0.557678, 0.000417078
+DEAL::0.0249049, 0.117851, 0.557678, 0.000417078
+DEAL::Outside quadrature
+DEAL::0.211325, 0.211325, 0.769002, 0.0366117
+DEAL::0.788675, 0.211325, 0.769002, 0.0366117
+DEAL::0.211325, 0.788675, 0.769002, 0.0366117
+DEAL::0.788675, 0.788675, 0.769002, 0.0366117
+DEAL::0.211325, 0.211325, 0.938104, 0.0366117
+DEAL::0.788675, 0.211325, 0.938104, 0.0366117
+DEAL::0.211325, 0.788675, 0.938104, 0.0366117
+DEAL::0.788675, 0.788675, 0.938104, 0.0366117
+DEAL::0.211325, 0.651151, 0.149429, 0.0390962
+DEAL::0.788675, 0.651151, 0.149429, 0.0390962
+DEAL::0.211325, 0.906526, 0.149429, 0.0390962
+DEAL::0.788675, 0.906526, 0.149429, 0.0390962
+DEAL::0.211325, 0.329176, 0.557678, 0.0751805
+DEAL::0.788675, 0.329176, 0.557678, 0.0751805
+DEAL::0.211325, 0.820253, 0.557678, 0.0751805
+DEAL::0.788675, 0.820253, 0.557678, 0.0751805
+DEAL::0.558205, 0.117851, 0.149429, 0.0276122
+DEAL::0.881621, 0.117851, 0.149429, 0.0276122
+DEAL::0.304271, 0.439826, 0.149429, 0.0434831
+DEAL::0.813580, 0.439826, 0.149429, 0.0434831
+DEAL::0.304271, 0.0315781, 0.557678, 0.0116512
+DEAL::0.813580, 0.0315781, 0.557678, 0.0116512
+DEAL::0.236230, 0.117851, 0.557678, 0.0127907
+DEAL::0.795348, 0.117851, 0.557678, 0.0127907
+DEAL::Surface quadrature
+DEAL::0.439826, 0.117851, 0.149429, 0.170753, 0.577350, 0.577350, 0.577350
+DEAL::0.117851, 0.439826, 0.149429, 0.170753, 0.577350, 0.577350, 0.577350
+DEAL::0.117851, 0.0315781, 0.557678, 0.0457532, 0.577350, 0.577350, 0.577350
+DEAL::0.0315781, 0.117851, 0.557678, 0.0457532, 0.577350, 0.577350, 0.577350
+DEAL::
+DEAL::test_epsilon_cut_at_bottom_corner
+DEAL::dim=3
+DEAL::Inside quadrature
+DEAL::Outside quadrature
+DEAL::0.211325, 0.211325, 0.211325, 0.125000
+DEAL::0.788675, 0.211325, 0.211325, 0.125000
+DEAL::0.211325, 0.788675, 0.211325, 0.125000
+DEAL::0.788675, 0.788675, 0.211325, 0.125000
+DEAL::0.211325, 0.211325, 0.788675, 0.125000
+DEAL::0.788675, 0.211325, 0.788675, 0.125000
+DEAL::0.211325, 0.788675, 0.788675, 0.125000
+DEAL::0.788675, 0.788675, 0.788675, 0.125000
+DEAL::Surface quadrature
+DEAL::
+DEAL::test_sphere_cutting_corner_exactly
+DEAL::dim=3
+DEAL::Inside quadrature
+DEAL::0.125000, 0.125000, 0.214407, 0.0311279
+DEAL::0.125000, 0.125000, 0.789501, 0.0311279
+DEAL::0.302831, 0.0528312, 0.220655, 0.00772008
+DEAL::0.302831, 0.0528312, 0.791175, 0.00772008
+DEAL::0.447169, 0.0528312, 0.231377, 0.00761387
+DEAL::0.447169, 0.0528312, 0.794048, 0.00761387
+DEAL::0.302831, 0.197169, 0.224225, 0.00768472
+DEAL::0.302831, 0.197169, 0.792132, 0.00768472
+DEAL::0.447169, 0.197169, 0.234959, 0.00757839
+DEAL::0.447169, 0.197169, 0.795008, 0.00757839
+DEAL::0.0528312, 0.302831, 0.220655, 0.00772008
+DEAL::0.0528312, 0.302831, 0.791175, 0.00772008
+DEAL::0.197169, 0.302831, 0.224225, 0.00768472
+DEAL::0.197169, 0.302831, 0.792132, 0.00768472
+DEAL::0.0528312, 0.447169, 0.231377, 0.00761387
+DEAL::0.0528312, 0.447169, 0.794048, 0.00761387
+DEAL::0.197169, 0.447169, 0.234959, 0.00757839
+DEAL::0.197169, 0.447169, 0.795008, 0.00757839
+DEAL::0.302831, 0.302831, 0.229459, 0.00763287
+DEAL::0.302831, 0.302831, 0.793534, 0.00763287
+DEAL::0.447169, 0.302831, 0.240211, 0.00752636
+DEAL::0.447169, 0.302831, 0.796415, 0.00752636
+DEAL::0.302831, 0.447169, 0.240211, 0.00752636
+DEAL::0.302831, 0.447169, 0.796415, 0.00752636
+DEAL::0.447169, 0.447169, 0.251000, 0.00741948
+DEAL::0.447169, 0.447169, 0.799306, 0.00741948
+DEAL::0.605662, 0.105662, 0.248812, 0.0297646
+DEAL::0.605662, 0.105662, 0.798720, 0.0297646
+DEAL::0.894338, 0.105662, 0.292317, 0.0280408
+DEAL::0.894338, 0.105662, 0.810377, 0.0280408
+DEAL::0.605662, 0.394338, 0.263246, 0.0291927
+DEAL::0.605662, 0.394338, 0.802587, 0.0291927
+DEAL::0.894338, 0.394338, 0.306956, 0.0274608
+DEAL::0.894338, 0.394338, 0.814299, 0.0274608
+DEAL::0.105662, 0.605662, 0.248812, 0.0297646
+DEAL::0.105662, 0.605662, 0.798720, 0.0297646
+DEAL::0.394338, 0.605662, 0.263246, 0.0291927
+DEAL::0.394338, 0.605662, 0.802587, 0.0291927
+DEAL::0.105662, 0.894338, 0.292317, 0.0280408
+DEAL::0.105662, 0.894338, 0.810377, 0.0280408
+DEAL::0.394338, 0.894338, 0.306956, 0.0274608
+DEAL::0.394338, 0.894338, 0.814299, 0.0274608
+DEAL::0.605662, 0.605662, 0.284500, 0.0283505
+DEAL::0.605662, 0.605662, 0.808282, 0.0283505
+DEAL::0.894338, 0.605662, 0.328517, 0.0266065
+DEAL::0.894338, 0.605662, 0.820077, 0.0266065
+DEAL::0.605662, 0.894338, 0.328517, 0.0266065
+DEAL::0.605662, 0.894338, 0.820077, 0.0266065
+DEAL::0.894338, 0.894338, 0.373180, 0.0248367
+DEAL::0.894338, 0.894338, 0.832044, 0.0248367
+DEAL::Outside quadrature
+DEAL::0.125000, 0.125000, 0.000825891, 0.000122130
+DEAL::0.125000, 0.125000, 0.00308227, 0.000122130
+DEAL::0.302831, 0.0528312, 0.00249992, 9.24198e-05
+DEAL::0.302831, 0.0528312, 0.00932981, 9.24198e-05
+DEAL::0.447169, 0.0528312, 0.00537287, 0.000198630
+DEAL::0.447169, 0.0528312, 0.0200518, 0.000198630
+DEAL::0.302831, 0.197169, 0.00345648, 0.000127783
+DEAL::0.302831, 0.197169, 0.0128997, 0.000127783
+DEAL::0.447169, 0.197169, 0.00633270, 0.000234115
+DEAL::0.447169, 0.197169, 0.0236340, 0.000234115
+DEAL::0.0528312, 0.302831, 0.00249992, 9.24198e-05
+DEAL::0.0528312, 0.302831, 0.00932981, 9.24198e-05
+DEAL::0.197169, 0.302831, 0.00345648, 0.000127783
+DEAL::0.197169, 0.302831, 0.0128997, 0.000127783
+DEAL::0.0528312, 0.447169, 0.00537287, 0.000198630
+DEAL::0.0528312, 0.447169, 0.0200518, 0.000198630
+DEAL::0.197169, 0.447169, 0.00633270, 0.000234115
+DEAL::0.197169, 0.447169, 0.0236340, 0.000234115
+DEAL::0.302831, 0.302831, 0.00485894, 0.000179631
+DEAL::0.302831, 0.302831, 0.0181338, 0.000179631
+DEAL::0.447169, 0.302831, 0.00773999, 0.000286141
+DEAL::0.447169, 0.302831, 0.0288860, 0.000286141
+DEAL::0.302831, 0.447169, 0.00773999, 0.000286141
+DEAL::0.302831, 0.447169, 0.0288860, 0.000286141
+DEAL::0.447169, 0.447169, 0.0106310, 0.000393018
+DEAL::0.447169, 0.447169, 0.0396753, 0.000393018
+DEAL::0.605662, 0.105662, 0.0100446, 0.00148535
+DEAL::0.605662, 0.105662, 0.0374868, 0.00148535
+DEAL::0.894338, 0.105662, 0.0217017, 0.00320918
+DEAL::0.894338, 0.105662, 0.0809920, 0.00320918
+DEAL::0.605662, 0.394338, 0.0139121, 0.00205728
+DEAL::0.605662, 0.394338, 0.0519208, 0.00205728
+DEAL::0.894338, 0.394338, 0.0256243, 0.00378924
+DEAL::0.894338, 0.394338, 0.0956313, 0.00378924
+DEAL::0.105662, 0.605662, 0.0100446, 0.00148535
+DEAL::0.105662, 0.605662, 0.0374868, 0.00148535
+DEAL::0.394338, 0.605662, 0.0139121, 0.00205728
+DEAL::0.394338, 0.605662, 0.0519208, 0.00205728
+DEAL::0.105662, 0.894338, 0.0217017, 0.00320918
+DEAL::0.105662, 0.894338, 0.0809920, 0.00320918
+DEAL::0.394338, 0.894338, 0.0256243, 0.00378924
+DEAL::0.394338, 0.894338, 0.0956313, 0.00378924
+DEAL::0.605662, 0.605662, 0.0196073, 0.00289946
+DEAL::0.605662, 0.605662, 0.0731755, 0.00289946
+DEAL::0.894338, 0.605662, 0.0314015, 0.00464354
+DEAL::0.894338, 0.605662, 0.117192, 0.00464354
+DEAL::0.605662, 0.894338, 0.0314015, 0.00464354
+DEAL::0.605662, 0.894338, 0.117192, 0.00464354
+DEAL::0.894338, 0.894338, 0.0433691, 0.00641327
+DEAL::0.894338, 0.894338, 0.161856, 0.00641327
+DEAL::Surface quadrature
+DEAL::0.125000, 0.125000, 0.00390816, 0.0625611, 0.0312500, 0.0312500, -0.999023
+DEAL::0.302831, 0.0528312, 0.0118297, 0.0156713, 0.0757078, 0.0132078, -0.997043
+DEAL::0.447169, 0.0528312, 0.0254247, 0.0157250, 0.111792, 0.0132078, -0.993644
+DEAL::0.302831, 0.197169, 0.0163562, 0.0156892, 0.0757078, 0.0492922, -0.995911
+DEAL::0.447169, 0.197169, 0.0299667, 0.0157429, 0.111792, 0.0492922, -0.992508
+DEAL::0.0528312, 0.302831, 0.0118297, 0.0156713, 0.0132078, 0.0757078, -0.997043
+DEAL::0.197169, 0.302831, 0.0163562, 0.0156892, 0.0492922, 0.0757078, -0.995911
+DEAL::0.0528312, 0.447169, 0.0254247, 0.0157250, 0.0132078, 0.111792, -0.993644
+DEAL::0.197169, 0.447169, 0.0299667, 0.0157429, 0.0492922, 0.111792, -0.992508
+DEAL::0.302831, 0.302831, 0.0229928, 0.0157153, 0.0757078, 0.0757078, -0.994252
+DEAL::0.447169, 0.302831, 0.0366260, 0.0157694, 0.111792, 0.0757078, -0.990843
+DEAL::0.302831, 0.447169, 0.0366260, 0.0157694, 0.0757078, 0.111792, -0.990843
+DEAL::0.447169, 0.447169, 0.0503063, 0.0158240, 0.111792, 0.111792, -0.987423
+DEAL::0.605662, 0.105662, 0.0475313, 0.0632516, 0.151416, 0.0264156, -0.988117
+DEAL::0.894338, 0.105662, 0.102694, 0.0641469, 0.223584, 0.0264156, -0.974327
+DEAL::0.605662, 0.394338, 0.0658329, 0.0635459, 0.151416, 0.0985844, -0.983542
+DEAL::0.894338, 0.394338, 0.121256, 0.0644538, 0.223584, 0.0985844, -0.969686
+DEAL::0.105662, 0.605662, 0.0475313, 0.0632516, 0.0264156, 0.151416, -0.988117
+DEAL::0.394338, 0.605662, 0.0658329, 0.0635459, 0.0985844, 0.151416, -0.983542
+DEAL::0.105662, 0.894338, 0.102694, 0.0641469, 0.0264156, 0.223584, -0.974327
+DEAL::0.394338, 0.894338, 0.121256, 0.0644538, 0.0985844, 0.223584, -0.969686
+DEAL::0.605662, 0.605662, 0.0927828, 0.0639842, 0.151416, 0.151416, -0.976804
+DEAL::0.894338, 0.605662, 0.148593, 0.0649113, 0.223584, 0.151416, -0.962852
+DEAL::0.605662, 0.894338, 0.148593, 0.0649113, 0.151416, 0.223584, -0.962852
+DEAL::0.894338, 0.894338, 0.205225, 0.0658801, 0.223584, 0.223584, -0.948694
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test that QuadratureGenerator clears its previously created quadratures
+ * when we call generate() again.
+ */
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include "deal.II/non_matching/quadrature_generator.h"
+
+#include "../tests.h"
+
+
+// Print the sizes of all the quadratures that QuadratureGenerator creates to
+// deallog.
+template <int dim>
+void
+print_n_quadrature_points(
+ const NonMatching::QuadratureGenerator<dim> &quadrature_generator)
+{
+ deallog << "inside " << quadrature_generator.get_inside_quadrature().size()
+ << std::endl;
+
+ deallog << "outside " << quadrature_generator.get_outside_quadrature().size()
+ << std::endl;
+
+ deallog << "surface " << quadrature_generator.get_surface_quadrature().size()
+ << std::endl;
+
+ deallog << std::endl;
+}
+
+
+
+// Call the QuadratureGenerator::generate with the same level set function
+// twice. Make sure that the sizes of the constructed quadratures are the same
+// both times. The purpose is to make sure that the previously created
+// quadratures have been cleared before we create the new ones.
+template <int dim>
+void
+test()
+{
+ deallog << "dim = " << dim << std::endl;
+
+ const hp::QCollection<1> q_collection(QGauss<1>(1));
+ NonMatching::QuadratureGenerator<dim> quadrature_generator(q_collection);
+
+ const Functions::ConstantFunction<dim> level_set(1);
+
+ const BoundingBox<dim> box = create_unit_bounding_box<dim>();
+
+ deallog << "quadrature sizes first call" << std::endl;
+ quadrature_generator.generate(level_set, box);
+ print_n_quadrature_points(quadrature_generator);
+
+ deallog << "quadrature sizes second call" << std::endl;
+ quadrature_generator.generate(level_set, box);
+ print_n_quadrature_points(quadrature_generator);
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ test<1>();
+ test<2>();
+ test<3>();
+}
--- /dev/null
+
+DEAL::dim = 1
+DEAL::quadrature sizes first call
+DEAL::inside 0
+DEAL::outside 1
+DEAL::surface 0
+DEAL::
+DEAL::quadrature sizes second call
+DEAL::inside 0
+DEAL::outside 1
+DEAL::surface 0
+DEAL::
+DEAL::dim = 2
+DEAL::quadrature sizes first call
+DEAL::inside 0
+DEAL::outside 1
+DEAL::surface 0
+DEAL::
+DEAL::quadrature sizes second call
+DEAL::inside 0
+DEAL::outside 1
+DEAL::surface 0
+DEAL::
+DEAL::dim = 3
+DEAL::quadrature sizes first call
+DEAL::inside 0
+DEAL::outside 1
+DEAL::surface 0
+DEAL::
+DEAL::quadrature sizes second call
+DEAL::inside 0
+DEAL::outside 1
+DEAL::surface 0
+DEAL::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/function_level_set.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <cmath>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+/**
+ * Compute the volume and surface area of a ball/sphere by setting up a
+ * level set function immersed in a background mesh, generating the
+ * quadrature rules, and summing the weights.
+ */
+template <int dim>
+void
+calculate_volume_and_surface_area()
+{
+ // Set up a background mesh
+ Triangulation<dim> triangulation;
+ const int n_subdivisions = 12;
+ const double gridsize = 2.07;
+ GridGenerator::subdivided_hyper_cube(triangulation,
+ n_subdivisions,
+ -gridsize / 2,
+ gridsize / 2);
+
+
+ // Description of the immersed domain.
+ const Functions::LevelSet::Sphere<dim> level_set;
+
+ // Create a quadrature generator.
+ const hp::QCollection<1> q_collection1D(QGauss<1>(2));
+ NonMatching::QuadratureGenerator<dim> quadrature_generator(q_collection1D);
+
+ // Go over all cells and compute the volume and surface area.
+ double surface_area = 0, volume = 0;
+ for (const auto cell : triangulation.active_cell_iterators())
+ {
+ // Create a box corresponding to the cell.
+ std::pair<Point<dim>, Point<dim>> lower_upper_corner;
+ lower_upper_corner.first = cell->vertex(0);
+ lower_upper_corner.second =
+ cell->vertex(GeometryInfo<dim>::vertices_per_cell - 1);
+ const BoundingBox<dim> box(lower_upper_corner);
+
+ // Generate immersed quadrature rules.
+ quadrature_generator.generate(level_set, box);
+
+ // Get the quadrature rules.
+ const Quadrature<dim> &inside_quadrature =
+ quadrature_generator.get_inside_quadrature();
+ const NonMatching::ImmersedSurfaceQuadrature<dim> &surface_quadrature =
+ quadrature_generator.get_surface_quadrature();
+
+ // Sum the weights to get the area/volume of the sphere.
+ for (unsigned int i = 0; i < inside_quadrature.size(); ++i)
+ volume += inside_quadrature.weight(i);
+
+ // Sum the weights to get the circumference/surface area of the sphere.
+ for (unsigned int i = 0; i < surface_quadrature.size(); ++i)
+ surface_area += surface_quadrature.weight(i);
+ }
+
+ deallog << "dim = " << dim << std::endl;
+ deallog << (2 == dim ? "area = " : "volume = ");
+ deallog << volume / M_PI << " * pi" << std::endl;
+ deallog << (2 == dim ? "circumference = " : "surface area = ");
+ deallog << surface_area / M_PI << " * pi" << std::endl;
+ deallog << std::endl;
+}
+
+
+int
+main()
+{
+ initlog();
+ calculate_volume_and_surface_area<2>();
+ calculate_volume_and_surface_area<3>();
+}
--- /dev/null
+
+DEAL::dim = 2
+DEAL::area = 1.00000 * pi
+DEAL::circumference = 1.99997 * pi
+DEAL::
+DEAL::dim = 3
+DEAL::volume = 1.33334 * pi
+DEAL::surface area = 3.99996 * pi
+DEAL::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_quadrature_printing_h_
+#define dealii_quadrature_printing_h_
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/non_matching/immersed_surface_quadrature.h>
+
+using namespace dealii;
+
+/*
+ * Print the incoming quadrature to deallog as comma separated values:
+ * point[0], ..., point[dim-1], weight
+ */
+template <int dim>
+void
+print_quadrature(const Quadrature<dim> &quadrature)
+{
+ for (unsigned int i = 0; i < quadrature.size(); ++i)
+ {
+ const Point<dim> &point = quadrature.point(i);
+ for (int d = 0; d < dim; d++)
+ deallog << point[d] << ", ";
+
+ deallog << quadrature.weight(i) << std::endl;
+ }
+}
+
+
+
+/*
+ * Print the incoming surface quadrature to deallog as comma separated values:
+ * p[0], ..., p[dim-1], weight, normal[0], ..., normal[dim-1]
+ */
+template <int dim>
+void
+print_surface_quadrature(
+ const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature)
+{
+ for (unsigned int i = 0; i < quadrature.size(); ++i)
+ {
+ const Point<dim> &point = quadrature.point(i);
+ for (int d = 0; d < dim; d++)
+ deallog << point[d] << ", ";
+
+ deallog << quadrature.weight(i);
+
+ const Tensor<1, dim> &normal = quadrature.normal_vector(i);
+ for (int d = 0; d < dim; d++)
+ deallog << ", " << normal[d];
+ deallog << std::endl;
+ }
+}
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test the RootFinder class in internal::QuadratureGeneratorImplementation.
+ */
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/function_level_set.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+using namespace NonMatching::internal::QuadratureGeneratorImplementation;
+
+// Use RootFinder to find the roots of the incoming functions over the interval
+// [0, 1]. Print the roots to deallog.
+void
+find_and_print_roots(
+ const std::vector<std::reference_wrapper<const Function<1>>> &functions)
+{
+ const BoundingBox<1> interval = create_unit_bounding_box<1>();
+
+ std::vector<double> roots;
+ RootFinder root_finder;
+ root_finder.find_roots(functions, interval, roots);
+
+ for (unsigned int i = 0; i < roots.size(); i++)
+ {
+ deallog << roots[i];
+ if (i < roots.size() - 1)
+ deallog << ", ";
+ }
+ deallog << std::endl;
+}
+
+
+
+// Test that the roots we get back from RootFinder are sorted and
+// not duplicated.
+//
+// Call find_roots with 3 linear functions f_i(x) = x - x_i,
+// where x_0 = 0.75, x_1 = 0.25, x_2 = 0.25
+// and check that RootFinder gives back the vector {0.25, 0.75}.
+void
+test_roots_sorted_not_duplicated()
+{
+ deallog << "test_roots_sorted_not_duplicated" << std::endl;
+
+ std::vector<Functions::LevelSet::Plane<1>> linear_functions;
+
+ const std::vector<double> roots = {.75, .25, .25};
+ for (unsigned int i = 0; i < roots.size(); ++i)
+ {
+ Tensor<1, 1> normal;
+ normal[0] = 1;
+ const Point<1> point(roots.at(i));
+ linear_functions.push_back(Functions::LevelSet::Plane<1>(point, normal));
+ }
+
+ const std::vector<std::reference_wrapper<const Function<1>>> functions(
+ linear_functions.begin(), linear_functions.end());
+
+ find_and_print_roots(functions);
+}
+
+
+
+/*
+ * The function:
+ * f(x) = C(x - x_0)^2 + y_0
+ */
+class QuadraticFunction : public Function<1>
+{
+public:
+ QuadraticFunction(const double C, const double x_0, const double y_0)
+ : C(C)
+ , x_0(x_0)
+ , y_0(y_0)
+ {}
+
+ double
+ value(const Point<1> &point, const unsigned int component = 0) const override
+ {
+ return C * std::pow(point(0) - x_0, 2) + y_0;
+ };
+
+ Tensor<1, 1>
+ gradient(const Point<1> & point,
+ const unsigned int component = 0) const override
+ {
+ Tensor<1, 1> grad;
+ grad[0] = 2 * C * (point(0) - x_0);
+
+ return grad;
+ };
+
+ SymmetricTensor<2, 1>
+ hessian(const Point<1> & point,
+ const unsigned int component = 0) const override
+ {
+ SymmetricTensor<2, 1> grad;
+ grad[0][0] = 2 * C;
+
+ return grad;
+ };
+
+private:
+ const double C;
+ const double x_0;
+ const double y_0;
+};
+
+
+
+// Test that RootFinder can find both roots of the function
+// f(x) = 4(x-0.5)^2 - 0.25
+// which are x_0 = 0.25 and x_1 = 0.75.
+void
+test_find_both_roots()
+{
+ deallog << "test_find_both_roots" << std::endl;
+
+ const QuadraticFunction function(4, 0.5, -0.25);
+
+ std::vector<std::reference_wrapper<const Function<1>>> functions;
+ functions.push_back(function);
+
+ find_and_print_roots(functions);
+}
+
+
+
+int
+main()
+{
+ initlog();
+ test_roots_sorted_not_duplicated();
+ deallog << std::endl;
+ test_find_both_roots();
+}
--- /dev/null
+
+DEAL::test_roots_sorted_not_duplicated
+DEAL::0.250000, 0.750000
+DEAL::
+DEAL::test_find_both_roots
+DEAL::0.250000, 0.750000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/**
+ * Test the function tensor_point_with_1D_quadrature()
+ * in NonMatching::internal::QuadratureGeneratorImplementation.
+ */
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include "../tests.h"
+
+#include "quadrature_printing.h"
+
+using namespace dealii;
+using namespace NonMatching::internal::QuadratureGeneratorImplementation;
+
+
+/*
+ * Set up a (dim-1)-dimensional point and a 1D-quadrature. Call
+ * tensor_point_with_1D_quadrature for each possible coordinate direction
+ * and print the resulting dim-dimensional quadrature.
+ */
+template <int dim>
+void
+create_and_output_quadrature_for_each_direction()
+{
+ deallog << "dim=" << dim << std::endl;
+
+ const unsigned int n_quadrature_points = 2;
+ const QGaussLobatto<1> quadrature1D(n_quadrature_points);
+ // Choose the points coordinates to something
+ // easily distinguished.
+ Point<dim - 1> point;
+ for (int i = 0; i < dim - 1; ++i)
+ {
+ point(i) = 10 * (i + 1);
+ }
+ // Both points in the 1D-quadrature have weight 1/2 so
+ // this should also be the weight of the points in the final
+ // quadrature.
+ const double weight = 5;
+ const double start = -1, end = 1;
+
+ for (int direction = 0; direction < dim; ++direction)
+ {
+ deallog << "direction=" << direction << std::endl;
+ ExtendableQuadrature<dim> result;
+ tensor_point_with_1D_quadrature(
+ point, weight, quadrature1D, start, end, direction, result);
+ print_quadrature(result);
+ deallog << std::endl;
+ }
+}
+
+
+
+int
+main()
+{
+ initlog();
+ create_and_output_quadrature_for_each_direction<1>();
+ create_and_output_quadrature_for_each_direction<2>();
+ create_and_output_quadrature_for_each_direction<3>();
+}
--- /dev/null
+
+DEAL::dim=1
+DEAL::direction=0
+DEAL::-1.00000, 5.00000
+DEAL::1.00000, 5.00000
+DEAL::
+DEAL::dim=2
+DEAL::direction=0
+DEAL::-1.00000, 10.0000, 5.00000
+DEAL::1.00000, 10.0000, 5.00000
+DEAL::
+DEAL::direction=1
+DEAL::10.0000, -1.00000, 5.00000
+DEAL::10.0000, 1.00000, 5.00000
+DEAL::
+DEAL::dim=3
+DEAL::direction=0
+DEAL::-1.00000, 10.0000, 20.0000, 5.00000
+DEAL::1.00000, 10.0000, 20.0000, 5.00000
+DEAL::
+DEAL::direction=1
+DEAL::20.0000, -1.00000, 10.0000, 5.00000
+DEAL::20.0000, 1.00000, 10.0000, 5.00000
+DEAL::
+DEAL::direction=2
+DEAL::10.0000, 20.0000, -1.00000, 5.00000
+DEAL::10.0000, 20.0000, 1.00000, 5.00000
+DEAL::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/**
+ * Test the class UpThroughDimensionCreator
+ * in NonMatching::internal::QuadratureGeneratorImplementation.
+ */
+
+#include <deal.II/base/function_level_set.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include "../tests.h"
+
+#include "quadrature_printing.h"
+
+using namespace dealii;
+using namespace NonMatching::internal::QuadratureGeneratorImplementation;
+
+
+// Print each quadrature in the incoming QPartitioning.
+template <int dim>
+void
+print(const QPartitioning<dim> &q_partitioning)
+{
+ deallog << "Negative" << std::endl;
+ print_quadrature(q_partitioning.negative);
+ deallog << "Positive" << std::endl;
+ print_quadrature(q_partitioning.positive);
+ deallog << "Indefinite" << std::endl;
+ print_quadrature(q_partitioning.indefinite);
+ deallog << "Surface" << std::endl;
+ print_surface_quadrature(q_partitioning.surface);
+}
+
+
+
+// Let the height function direction be dim - 1 and let the lower dimensional
+// quadrature contain a single point. Call UpThroughDimensionCreator with the
+// incoming level set function over the unit box to generate a QPartitioning.
+// Print the quadratures in the partitioning to make sure it is correct.
+template <int dim>
+void
+create_and_print_partitioning(const Function<dim> &level_set)
+{
+ const hp::QCollection<1> q_collection1D(QGauss<1>(2));
+ const NonMatching::AdditionalQGeneratorData additional_data;
+
+ UpThroughDimensionCreator<dim, dim> up_through_dimension_creator(
+ q_collection1D, additional_data);
+
+ std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
+ level_sets.push_back(level_set);
+ const BoundingBox<dim> box = create_unit_bounding_box<dim>();
+ const QMidpoint<dim - 1> low_dim_quadrature;
+ const unsigned int height_function_direction = dim - 1;
+
+ QPartitioning<dim> q_partitioning;
+ up_through_dimension_creator.generate(level_sets,
+ box,
+ low_dim_quadrature,
+ height_function_direction,
+ q_partitioning);
+ print(q_partitioning);
+}
+
+
+
+// Set up a level set function with the zero contour along, x_{dim-1} = 0.5
+// Call create_and_print_partitioning to test that the points are added as
+// expected:
+// "negative" points should have x_{dim-1} \in (0, 0.5)
+// "positive" points should have x_{dim-1} \in (0.5, 1)
+// surface points should have x_{dim-1} = 0.5
+template <int dim>
+void
+test_cut_through_center()
+{
+ deallog << "test_cut_through_center" << std::endl;
+ deallog << std::endl;
+
+ Point<dim> point_through_plane = .5 * Point<dim>::unit_vector(dim - 1);
+ Tensor<1, dim> plane_normal = Point<dim>::unit_vector(dim - 1);
+ const Functions::LevelSet::Plane<dim> level_set(point_through_plane,
+ plane_normal);
+
+ create_and_print_partitioning(level_set);
+}
+
+
+
+// Fabricate the case when we have missed roots when creating the quadrature in
+// the lower dimensions. See the comment in the implementation of
+// UpThroughDimensionCreator::create_surface_point(..).
+//
+// In this test, the zero contour goes outside the cell but close to the "bottom
+// face" at x_{dim-1} = 0. Check that the surface quadrature points gets placed
+// on x_{dim-1} = 0. This is the "least bad" option we have.
+template <int dim>
+void
+test_missed_roots_on_bottom_face()
+{
+ deallog << "test_missed_roots_on_bottom_face" << std::endl;
+ deallog << std::endl;
+
+ const Tensor<1, dim> plane_normal = Point<dim>::unit_vector(dim - 1);
+ Point<dim> point_in_plane;
+ point_in_plane[dim - 1] = -.1;
+ const Functions::LevelSet::Plane<dim> level_set(point_in_plane, plane_normal);
+
+ create_and_print_partitioning(level_set);
+}
+
+
+
+// Same test as above, but the zero contour just outside the "top face" at
+// x_{dim-1} = 1. Check that the surface quadrature points are placed on
+// x_{dim-1} = 1.
+template <int dim>
+void
+test_missed_roots_on_top_face()
+{
+ deallog << "test_missed_roots_on_top_face" << std::endl;
+ deallog << std::endl;
+
+ const Tensor<1, dim> plane_normal = Point<dim>::unit_vector(dim - 1);
+ Point<dim> point_in_plane;
+ point_in_plane[dim - 1] = 1.1;
+ const Functions::LevelSet::Plane<dim> level_set(point_in_plane, plane_normal);
+
+ create_and_print_partitioning(level_set);
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ const int dim = 2;
+
+ test_cut_through_center<dim>();
+ deallog << std::endl;
+
+ test_missed_roots_on_bottom_face<dim>();
+ deallog << std::endl;
+
+ test_missed_roots_on_top_face<dim>();
+ deallog << std::endl;
+}
--- /dev/null
+
+DEAL::test_cut_through_center
+DEAL::
+DEAL::Negative
+DEAL::0.500000, 0.105662, 0.250000
+DEAL::0.500000, 0.394338, 0.250000
+DEAL::Positive
+DEAL::0.500000, 0.605662, 0.250000
+DEAL::0.500000, 0.894338, 0.250000
+DEAL::Indefinite
+DEAL::Surface
+DEAL::0.500000, 0.500000, 1.00000, 0.00000, 1.00000
+DEAL::
+DEAL::test_missed_roots_on_bottom_face
+DEAL::
+DEAL::Negative
+DEAL::Positive
+DEAL::0.500000, 0.211325, 0.500000
+DEAL::0.500000, 0.788675, 0.500000
+DEAL::Indefinite
+DEAL::Surface
+DEAL::0.500000, 0.00000, 1.00000, 0.00000, 1.00000
+DEAL::
+DEAL::test_missed_roots_on_top_face
+DEAL::
+DEAL::Negative
+DEAL::0.500000, 0.211325, 0.500000
+DEAL::0.500000, 0.788675, 0.500000
+DEAL::Positive
+DEAL::Indefinite
+DEAL::Surface
+DEAL::0.500000, 1.00000, 1.00000, 0.00000, 1.00000
+DEAL::
projects / dealii.git / commitdiff
