]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add class generating high-order immersed quadrature rules 12125/head
authorSimon Sticko <simon@sticko.se>
Mon, 3 May 2021 09:46:14 +0000 (11:46 +0200)
committerSimon Sticko <simon@sticko.se>
Thu, 10 Jun 2021 05:53:13 +0000 (07:53 +0200)
Add a class NonMatching::QuadratureGenerator that generates high-order
immersed quadrature rules over a BoundingBox, B, when the domain is
described by a level set function, \psi. Thus creating quadrature rules
for the following 3 regions of the box:

{x \in B : \psi(x) < 0},
{x \in B : \psi(x) > 0},
{x \in B : \psi(x) = 0}.

24 files changed:
doc/doxygen/images/immersed_quadratures.svg [new file with mode: 0644]
include/deal.II/non_matching/quadrature_generator.h [new file with mode: 0644]
source/non_matching/CMakeLists.txt
source/non_matching/quadrature_generator.cc [new file with mode: 0644]
source/non_matching/quadrature_generator.inst.in [new file with mode: 0644]
tests/non_matching/find_best_height_direction.cc [new file with mode: 0644]
tests/non_matching/find_best_height_direction.output [new file with mode: 0644]
tests/non_matching/find_extreme_values.cc [new file with mode: 0644]
tests/non_matching/find_extreme_values.output [new file with mode: 0644]
tests/non_matching/pointwise_definiteness.cc [new file with mode: 0644]
tests/non_matching/pointwise_definiteness.output [new file with mode: 0644]
tests/non_matching/quadrature_generator.cc [new file with mode: 0644]
tests/non_matching/quadrature_generator.output [new file with mode: 0644]
tests/non_matching/quadrature_generator_clears_between_calls.cc [new file with mode: 0644]
tests/non_matching/quadrature_generator_clears_between_calls.output [new file with mode: 0644]
tests/non_matching/quadrature_generator_sphere.cc [new file with mode: 0644]
tests/non_matching/quadrature_generator_sphere.output [new file with mode: 0644]
tests/non_matching/quadrature_printing.h [new file with mode: 0644]
tests/non_matching/root_finder.cc [new file with mode: 0644]
tests/non_matching/root_finder.output [new file with mode: 0644]
tests/non_matching/tensor_point_with_1D_quadrature.cc [new file with mode: 0644]
tests/non_matching/tensor_point_with_1D_quadrature.output [new file with mode: 0644]
tests/non_matching/up_through_dimension_creator.cc [new file with mode: 0644]
tests/non_matching/up_through_dimension_creator.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/immersed_quadratures.svg b/doc/doxygen/images/immersed_quadratures.svg
new file mode 100644 (file)
index 0000000..063bb8b
--- /dev/null
@@ -0,0 +1,462 @@
+<?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 &lt; 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 &gt; 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>
diff --git a/include/deal.II/non_matching/quadrature_generator.h b/include/deal.II/non_matching/quadrature_generator.h
new file mode 100644 (file)
index 0000000..6d105f2
--- /dev/null
@@ -0,0 +1,1046 @@
+// ---------------------------------------------------------------------
+//
+// 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 */
index 8ea4a0eb48d8f22646dbaa476556540ed7955838..2d1b9a8a76f79c553ce39f552d45c2d498ceb1ce 100644 (file)
 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
   )
 
diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc
new file mode 100644 (file)
index 0000000..ea0245e
--- /dev/null
@@ -0,0 +1,1358 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/source/non_matching/quadrature_generator.inst.in b/source/non_matching/quadrature_generator.inst.in
new file mode 100644 (file)
index 0000000..736ff97
--- /dev/null
@@ -0,0 +1,74 @@
+// ---------------------------------------------------------------------
+//
+// 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
+  }
diff --git a/tests/non_matching/find_best_height_direction.cc b/tests/non_matching/find_best_height_direction.cc
new file mode 100644 (file)
index 0000000..109f221
--- /dev/null
@@ -0,0 +1,106 @@
+// ---------------------------------------------------------------------
+//
+// 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();
+}
diff --git a/tests/non_matching/find_best_height_direction.output b/tests/non_matching/find_best_height_direction.output
new file mode 100644 (file)
index 0000000..43736c2
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::test_ignores_definite_functions
+DEAL::OK
+DEAL::
+DEAL::test_find_best_height_direction
+DEAL::height direction = 1
+DEAL::min_abs_dfdx = 5.00000
diff --git a/tests/non_matching/find_extreme_values.cc b/tests/non_matching/find_extreme_values.cc
new file mode 100644 (file)
index 0000000..39fea77
--- /dev/null
@@ -0,0 +1,85 @@
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
diff --git a/tests/non_matching/find_extreme_values.output b/tests/non_matching/find_extreme_values.output
new file mode 100644 (file)
index 0000000..3c68e29
--- /dev/null
@@ -0,0 +1,8 @@
+
+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
diff --git a/tests/non_matching/pointwise_definiteness.cc b/tests/non_matching/pointwise_definiteness.cc
new file mode 100644 (file)
index 0000000..a8123fa
--- /dev/null
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// 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;
+}
diff --git a/tests/non_matching/pointwise_definiteness.output b/tests/non_matching/pointwise_definiteness.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/non_matching/quadrature_generator.cc b/tests/non_matching/quadrature_generator.cc
new file mode 100644 (file)
index 0000000..e98a3cc
--- /dev/null
@@ -0,0 +1,302 @@
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
diff --git a/tests/non_matching/quadrature_generator.output b/tests/non_matching/quadrature_generator.output
new file mode 100644 (file)
index 0000000..20f5130
--- /dev/null
@@ -0,0 +1,389 @@
+
+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
diff --git a/tests/non_matching/quadrature_generator_clears_between_calls.cc b/tests/non_matching/quadrature_generator_clears_between_calls.cc
new file mode 100644 (file)
index 0000000..ee45b23
--- /dev/null
@@ -0,0 +1,88 @@
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
diff --git a/tests/non_matching/quadrature_generator_clears_between_calls.output b/tests/non_matching/quadrature_generator_clears_between_calls.output
new file mode 100644 (file)
index 0000000..b12c146
--- /dev/null
@@ -0,0 +1,34 @@
+
+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::
diff --git a/tests/non_matching/quadrature_generator_sphere.cc b/tests/non_matching/quadrature_generator_sphere.cc
new file mode 100644 (file)
index 0000000..4ed0361
--- /dev/null
@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
diff --git a/tests/non_matching/quadrature_generator_sphere.output b/tests/non_matching/quadrature_generator_sphere.output
new file mode 100644 (file)
index 0000000..caf1f3f
--- /dev/null
@@ -0,0 +1,9 @@
+
+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::
diff --git a/tests/non_matching/quadrature_printing.h b/tests/non_matching/quadrature_printing.h
new file mode 100644 (file)
index 0000000..004c8fc
--- /dev/null
@@ -0,0 +1,70 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/tests/non_matching/root_finder.cc b/tests/non_matching/root_finder.cc
new file mode 100644 (file)
index 0000000..5b1f98c
--- /dev/null
@@ -0,0 +1,160 @@
+// ---------------------------------------------------------------------
+//
+// 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();
+}
diff --git a/tests/non_matching/root_finder.output b/tests/non_matching/root_finder.output
new file mode 100644 (file)
index 0000000..e4202fe
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::test_roots_sorted_not_duplicated
+DEAL::0.250000, 0.750000
+DEAL::
+DEAL::test_find_both_roots
+DEAL::0.250000, 0.750000
diff --git a/tests/non_matching/tensor_point_with_1D_quadrature.cc b/tests/non_matching/tensor_point_with_1D_quadrature.cc
new file mode 100644 (file)
index 0000000..3f2d3a1
--- /dev/null
@@ -0,0 +1,79 @@
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
diff --git a/tests/non_matching/tensor_point_with_1D_quadrature.output b/tests/non_matching/tensor_point_with_1D_quadrature.output
new file mode 100644 (file)
index 0000000..72077f8
--- /dev/null
@@ -0,0 +1,28 @@
+
+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::
diff --git a/tests/non_matching/up_through_dimension_creator.cc b/tests/non_matching/up_through_dimension_creator.cc
new file mode 100644 (file)
index 0000000..653be41
--- /dev/null
@@ -0,0 +1,166 @@
+// ---------------------------------------------------------------------
+//
+// 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;
+}
diff --git a/tests/non_matching/up_through_dimension_creator.output b/tests/non_matching/up_through_dimension_creator.output
new file mode 100644 (file)
index 0000000..336cbb1
--- /dev/null
@@ -0,0 +1,33 @@
+
+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::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.