]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add class NonMatching::FEInterfaceValues 13657/head
authorSimon Sticko <simon@sticko.se>
Fri, 29 Apr 2022 08:29:56 +0000 (10:29 +0200)
committerSimon Sticko <simon@sticko.se>
Mon, 9 May 2022 09:18:20 +0000 (11:18 +0200)
This class works the same way as NonMatching::FEValues but for
assembling interface terms on faces. This is done using the
DisceteFaceQuadratureGenerator class to create quadrature rules on
the intersected faces and then creating dealii::FEInterfaceValues
objects with these quadratures.

doc/doxygen/images/non_matching_fe_interface_values.svg [new file with mode: 0644]
include/deal.II/non_matching/fe_values.h
source/non_matching/fe_values.cc
source/non_matching/fe_values.inst.in
tests/non_matching/fe_interface_values.cc [new file with mode: 0644]
tests/non_matching/fe_interface_values.output [new file with mode: 0644]
tests/non_matching/fe_interface_values_non_standard_meshes.cc [new file with mode: 0644]
tests/non_matching/fe_interface_values_non_standard_meshes.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/non_matching_fe_interface_values.svg b/doc/doxygen/images/non_matching_fe_interface_values.svg
new file mode 100644 (file)
index 0000000..3c7de30
--- /dev/null
@@ -0,0 +1,591 @@
+<?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="329.76pt" version="1.1" viewBox="0 0 444.96 329.76" width="444.96pt" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
+ <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 329.76 
+L 444.96 329.76 
+L 444.96 0 
+L 0 0 
+z
+" style="fill:#ffffff;"/>
+  </g>
+  <g id="patch_2">
+   <path d="M 7.2 322.56 
+L 214.92 322.56 
+L 214.92 7.2 
+L 7.2 7.2 
+z
+" style="fill:#ffffff;"/>
+  </g>
+  <g id="axes_1">
+   <g id="line2d_1">
+    <path clip-path="url(#pea105021c5)" d="M 37.170217 251.389514 
+L 35.835673 171.475172 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_2">
+    <path clip-path="url(#pea105021c5)" d="M 37.170217 251.389514 
+L 94.862881 197.862086 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_3">
+    <path clip-path="url(#pea105021c5)" d="M 37.170217 251.389514 
+L 135.431156 286.45657 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_4">
+    <path clip-path="url(#pea105021c5)" d="M 35.835673 171.475172 
+L 94.557491 123.19876 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_5">
+    <path clip-path="url(#pea105021c5)" d="M 35.835673 171.475172 
+L 135.825216 203.170777 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_6">
+    <path clip-path="url(#pea105021c5)" d="M 94.862881 197.862086 
+L 94.557491 123.19876 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_7">
+    <path clip-path="url(#pea105021c5)" d="M 94.862881 197.862086 
+L 188.013319 227.73846 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_8">
+    <path clip-path="url(#pea105021c5)" d="M 94.557491 123.19876 
+L 189.259858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_9">
+    <path clip-path="url(#pea105021c5)" d="M 135.431156 286.45657 
+L 135.825216 203.170777 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_10">
+    <path clip-path="url(#pea105021c5)" d="M 135.431156 286.45657 
+L 188.013319 227.73846 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_11">
+    <path clip-path="url(#pea105021c5)" d="M 135.825216 203.170777 
+L 189.259858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_12">
+    <path clip-path="url(#pea105021c5)" d="M 188.013319 227.73846 
+L 189.259858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_13">
+    <path clip-path="url(#pea105021c5)" d="M 35.835673 171.475172 
+L 34.453864 88.730522 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_14">
+    <path clip-path="url(#pea105021c5)" d="M 35.835673 171.475172 
+L 94.557491 123.19876 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_15">
+    <path clip-path="url(#pea105021c5)" d="M 35.835673 171.475172 
+L 135.825216 203.170777 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_16">
+    <path clip-path="url(#pea105021c5)" d="M 34.453864 88.730522 
+L 94.242127 46.09662 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_17">
+    <path clip-path="url(#pea105021c5)" d="M 34.453864 88.730522 
+L 136.233946 116.784417 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_18">
+    <path clip-path="url(#pea105021c5)" d="M 94.557491 123.19876 
+L 94.242127 46.09662 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_19">
+    <path clip-path="url(#pea105021c5)" d="M 94.557491 123.19876 
+L 189.259858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_20">
+    <path clip-path="url(#pea105021c5)" d="M 94.242127 46.09662 
+L 190.549027 69.864771 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_21">
+    <path clip-path="url(#pea105021c5)" d="M 135.825216 203.170777 
+L 136.233946 116.784417 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_22">
+    <path clip-path="url(#pea105021c5)" d="M 135.825216 203.170777 
+L 189.259858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_23">
+    <path clip-path="url(#pea105021c5)" d="M 136.233946 116.784417 
+L 190.549027 69.864771 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_24">
+    <path clip-path="url(#pea105021c5)" d="M 189.259858 150.128693 
+L 190.549027 69.864771 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_25">
+    <path clip-path="url(#pea105021c5)" d="M 117.545578 196.994545 
+L 117.822478 195.772678 
+L 118.07038 194.553493 
+L 118.289403 193.337176 
+L 118.479673 192.123913 
+L 118.641318 190.913886 
+L 118.774471 189.707274 
+L 118.879267 188.504256 
+L 118.955849 187.305006 
+L 119.004359 186.109695 
+L 119.024944 184.918494 
+L 119.017757 183.731568 
+L 118.982951 182.549082 
+L 118.920684 181.371198 
+L 118.831115 180.198075 
+L 118.714409 179.029869 
+L 118.570732 177.866733 
+L 118.400254 176.70882 
+L 118.203145 175.556278 
+L 117.97958 174.409253 
+L 117.729737 173.26789 
+L 117.453794 172.132328 
+L 117.151933 171.002707 
+L 116.824337 169.879164 
+L 116.471192 168.761832 
+L 116.092687 167.650842 
+L 115.68901 166.546323 
+L 115.260352 165.448403 
+L 114.806908 164.357205 
+L 114.328872 163.272852 
+L 113.82644 162.195462 
+L 113.299811 161.125153 
+L 112.749182 160.062041 
+L 112.174755 159.006237 
+L 111.576732 157.957852 
+L 110.955316 156.916995 
+L 110.31071 155.883771 
+L 109.643119 154.858285 
+L 108.952751 153.840638 
+L 108.239811 152.83093 
+L 107.504507 151.829259 
+L 106.747049 150.83572 
+L 105.967644 149.850407 
+L 105.166504 148.87341 
+L 104.343838 147.90482 
+L 103.499857 146.944725 
+L 102.634774 145.993209 
+L 101.748799 145.050357 
+L 100.842145 144.11625 
+L 99.915026 143.190968 
+L 98.967652 142.274591 
+L 98.000238 141.367193 
+L 97.012997 140.46885 
+L 96.006143 139.579635 
+L 94.979888 138.699618 
+L 93.934447 137.82887 
+L 92.870032 136.967458 
+L 91.786858 136.115448 
+L 90.685139 135.272906 
+L 89.565087 134.439893 
+L 88.426915 133.616473 
+L 87.270838 132.802704 
+L 86.097068 131.998645 
+L 84.905818 131.204354 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:1.1;"/>
+   </g>
+   <g id="Path3DCollection_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="C0_0_966171e59e"/>
+    </defs>
+    <g clip-path="url(#pea105021c5)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="61.596309" xlink:href="#C0_0_966171e59e" y="167.122462"/>
+    </g>
+    <g clip-path="url(#pea105021c5)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="103.318521" xlink:href="#C0_0_966171e59e" y="180.089042"/>
+    </g>
+    <g clip-path="url(#pea105021c5)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="81.15241" xlink:href="#C0_0_966171e59e" y="141.062131"/>
+    </g>
+    <g clip-path="url(#pea105021c5)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="98.077145" xlink:href="#C0_0_966171e59e" y="146.055019"/>
+    </g>
+   </g>
+   <g id="text_1">
+    <!-- $\psi &lt; 0$ -->
+    <defs>
+     <path d="M 47 66.59375 
+C 47 66.6875 47.40625 67.984375 47.40625 68.09375 
+C 47.40625 69 46.59375 69 46.296875 69 
+C 45.40625 69 45.296875 68.484375 44.90625 66.984375 
+L 28.40625 1.1875 
+C 20.09375 2.1875 17 6.375 17 12.359375 
+C 17 14.5625 17 16.859375 21.703125 29.21875 
+C 23 32.8125 23.5 34.109375 23.5 35.90625 
+C 23.5 40.40625 20.296875 44 15.59375 44 
+C 6.40625 44 2.703125 29.53125 2.703125 28.734375 
+C 2.703125 28.328125 3.09375 27.828125 3.796875 27.828125 
+C 4.703125 27.828125 4.796875 28.234375 5.203125 29.625 
+C 7.59375 38.296875 11.59375 41.984375 15.296875 41.984375 
+C 16.203125 41.984375 17.90625 41.890625 17.90625 38.703125 
+C 17.90625 38.296875 17.90625 36.109375 16.203125 31.71875 
+C 10.796875 17.546875 10.796875 15.359375 10.796875 13.0625 
+C 10.796875 3.484375 18.796875 -0.203125 27.796875 -0.90625 
+C 27 -3.9375 26.296875 -7.0625 25.5 -10.09375 
+C 23.90625 -15.953125 23.203125 -18.703125 23.203125 -19.109375 
+C 23.203125 -20 24 -20 24.296875 -20 
+C 24.5 -20 24.90625 -20 25.09375 -19.609375 
+C 25.5 -19.203125 29.59375 -2.765625 30 -1 
+C 33.703125 -1 41.59375 -1 50.59375 8.28125 
+C 53.90625 11.859375 56.90625 16.453125 58.59375 20.75 
+C 59.59375 23.34375 62 32.21875 62 37.3125 
+C 62 43.296875 59 44 58 44 
+C 55.59375 44 53.40625 41.59375 53.40625 39.59375 
+C 53.40625 38.40625 54.09375 37.703125 54.5 37.3125 
+C 55.40625 36.40625 58.09375 33.71875 58.09375 28.53125 
+C 58.09375 24.9375 56.09375 17.546875 49.703125 10.375 
+C 41.296875 1 33.59375 1 30.59375 1 
+z
+" id="CMMI12-32"/>
+     <path d="M 65.90625 48.96875 
+C 67.703125 49.78125 67.90625 50.484375 67.90625 51.09375 
+C 67.90625 52.1875 67.09375 53 66 53 
+C 65.796875 53 65.703125 52.890625 64.296875 52.296875 
+L 10.203125 27.15625 
+C 8.40625 26.359375 8.203125 25.640625 8.203125 25.046875 
+C 8.203125 24.328125 8.296875 23.734375 10.203125 22.8125 
+L 64.296875 -2.3125 
+C 65.59375 -2.90625 65.796875 -3 66 -3 
+C 67.09375 -3 67.90625 -2.203125 67.90625 -1.109375 
+C 67.90625 -0.5 67.703125 0.203125 65.90625 1.015625 
+L 14.40625 24.9375 
+z
+" id="CMMI12-60"/>
+     <path d="M 42 31.84375 
+C 42 37.96875 41.90625 48.421875 37.703125 56.453125 
+C 34 63.484375 28.09375 66 22.90625 66 
+C 18.09375 66 12 63.78125 8.203125 56.5625 
+C 4.203125 49.015625 3.796875 39.671875 3.796875 31.84375 
+C 3.796875 26.109375 3.90625 17.375 7 9.734375 
+C 11.296875 -0.609375 19 -2 22.90625 -2 
+C 27.5 -2 34.5 -0.109375 38.59375 9.4375 
+C 41.59375 16.375 42 24.5 42 31.84375 
+z
+M 22.90625 -0.40625 
+C 16.5 -0.40625 12.703125 5.125 11.296875 12.75 
+C 10.203125 18.6875 10.203125 27.328125 10.203125 32.953125 
+C 10.203125 40.6875 10.203125 47.109375 11.5 53.234375 
+C 13.40625 61.78125 19 64.390625 22.90625 64.390625 
+C 27 64.390625 32.296875 61.671875 34.203125 53.4375 
+C 35.5 47.71875 35.59375 40.984375 35.59375 32.953125 
+C 35.59375 26.421875 35.59375 18.375 34.40625 12.453125 
+C 32.296875 1.5 26.40625 -0.40625 22.90625 -0.40625 
+z
+" id="CMR17-48"/>
+    </defs>
+    <g transform="translate(90.312965 26.4288)scale(0.168 -0.168)">
+     <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="patch_3">
+   <path d="M 230.04 322.56 
+L 437.76 322.56 
+L 437.76 7.2 
+L 230.04 7.2 
+z
+" style="fill:#ffffff;"/>
+  </g>
+  <g id="axes_2">
+   <g id="line2d_26">
+    <path clip-path="url(#p0a8decf98d)" d="M 260.010217 251.389514 
+L 258.675673 171.475172 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_27">
+    <path clip-path="url(#p0a8decf98d)" d="M 260.010217 251.389514 
+L 317.702881 197.862086 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_28">
+    <path clip-path="url(#p0a8decf98d)" d="M 260.010217 251.389514 
+L 358.271156 286.45657 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_29">
+    <path clip-path="url(#p0a8decf98d)" d="M 258.675673 171.475172 
+L 317.397491 123.19876 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_30">
+    <path clip-path="url(#p0a8decf98d)" d="M 258.675673 171.475172 
+L 358.665216 203.170777 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_31">
+    <path clip-path="url(#p0a8decf98d)" d="M 317.702881 197.862086 
+L 317.397491 123.19876 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_32">
+    <path clip-path="url(#p0a8decf98d)" d="M 317.702881 197.862086 
+L 410.853319 227.73846 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_33">
+    <path clip-path="url(#p0a8decf98d)" d="M 317.397491 123.19876 
+L 412.099858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_34">
+    <path clip-path="url(#p0a8decf98d)" d="M 358.271156 286.45657 
+L 358.665216 203.170777 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_35">
+    <path clip-path="url(#p0a8decf98d)" d="M 358.271156 286.45657 
+L 410.853319 227.73846 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_36">
+    <path clip-path="url(#p0a8decf98d)" d="M 358.665216 203.170777 
+L 412.099858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_37">
+    <path clip-path="url(#p0a8decf98d)" d="M 410.853319 227.73846 
+L 412.099858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_38">
+    <path clip-path="url(#p0a8decf98d)" d="M 258.675673 171.475172 
+L 257.293864 88.730522 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_39">
+    <path clip-path="url(#p0a8decf98d)" d="M 258.675673 171.475172 
+L 317.397491 123.19876 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_40">
+    <path clip-path="url(#p0a8decf98d)" d="M 258.675673 171.475172 
+L 358.665216 203.170777 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_41">
+    <path clip-path="url(#p0a8decf98d)" d="M 257.293864 88.730522 
+L 317.082127 46.09662 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_42">
+    <path clip-path="url(#p0a8decf98d)" d="M 257.293864 88.730522 
+L 359.073946 116.784417 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_43">
+    <path clip-path="url(#p0a8decf98d)" d="M 317.397491 123.19876 
+L 317.082127 46.09662 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_44">
+    <path clip-path="url(#p0a8decf98d)" d="M 317.397491 123.19876 
+L 412.099858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_45">
+    <path clip-path="url(#p0a8decf98d)" d="M 317.082127 46.09662 
+L 413.389027 69.864771 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_46">
+    <path clip-path="url(#p0a8decf98d)" d="M 358.665216 203.170777 
+L 359.073946 116.784417 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_47">
+    <path clip-path="url(#p0a8decf98d)" d="M 358.665216 203.170777 
+L 412.099858 150.128693 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_48">
+    <path clip-path="url(#p0a8decf98d)" d="M 359.073946 116.784417 
+L 413.389027 69.864771 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_49">
+    <path clip-path="url(#p0a8decf98d)" d="M 412.099858 150.128693 
+L 413.389027 69.864771 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:0.9;"/>
+   </g>
+   <g id="line2d_50">
+    <path clip-path="url(#p0a8decf98d)" d="M 340.385578 196.994545 
+L 340.662478 195.772678 
+L 340.91038 194.553493 
+L 341.129403 193.337176 
+L 341.319673 192.123913 
+L 341.481318 190.913886 
+L 341.614471 189.707274 
+L 341.719267 188.504256 
+L 341.795849 187.305006 
+L 341.844359 186.109695 
+L 341.864944 184.918494 
+L 341.857757 183.731568 
+L 341.822951 182.549082 
+L 341.760684 181.371198 
+L 341.671115 180.198075 
+L 341.554409 179.029869 
+L 341.410732 177.866733 
+L 341.240254 176.70882 
+L 341.043145 175.556278 
+L 340.81958 174.409253 
+L 340.569737 173.26789 
+L 340.293794 172.132328 
+L 339.991933 171.002707 
+L 339.664337 169.879164 
+L 339.311192 168.761832 
+L 338.932687 167.650842 
+L 338.52901 166.546323 
+L 338.100352 165.448403 
+L 337.646908 164.357205 
+L 337.168872 163.272852 
+L 336.66644 162.195462 
+L 336.139811 161.125153 
+L 335.589182 160.062041 
+L 335.014755 159.006237 
+L 334.416732 157.957852 
+L 333.795316 156.916995 
+L 333.15071 155.883771 
+L 332.483119 154.858285 
+L 331.792751 153.840638 
+L 331.079811 152.83093 
+L 330.344507 151.829259 
+L 329.587049 150.83572 
+L 328.807644 149.850407 
+L 328.006504 148.87341 
+L 327.183838 147.90482 
+L 326.339857 146.944725 
+L 325.474774 145.993209 
+L 324.588799 145.050357 
+L 323.682145 144.11625 
+L 322.755026 143.190968 
+L 321.807652 142.274591 
+L 320.840238 141.367193 
+L 319.852997 140.46885 
+L 318.846143 139.579635 
+L 317.819888 138.699618 
+L 316.774447 137.82887 
+L 315.710032 136.967458 
+L 314.626858 136.115448 
+L 313.525139 135.272906 
+L 312.405087 134.439893 
+L 311.266915 133.616473 
+L 310.110838 132.802704 
+L 308.937068 131.998645 
+L 307.745818 131.204354 
+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-width:1.1;"/>
+   </g>
+   <g id="Path3DCollection_2">
+    <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="C1_0_a2ed93b880"/>
+    </defs>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="329.216368" xlink:href="#C1_0_a2ed93b880" y="135.079927"/>
+    </g>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="384.283505" xlink:href="#C1_0_a2ed93b880" y="150.966563"/>
+    </g>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="334.715774" xlink:href="#C1_0_a2ed93b880" y="130.38843"/>
+    </g>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="389.48281" xlink:href="#C1_0_a2ed93b880" y="146.022193"/>
+    </g>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="347.462302" xlink:href="#C1_0_a2ed93b880" y="186.709908"/>
+    </g>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="362.910811" xlink:href="#C1_0_a2ed93b880" y="191.511051"/>
+    </g>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="341.025463" xlink:href="#C1_0_a2ed93b880" y="151.987082"/>
+    </g>
+    <g clip-path="url(#p0a8decf98d)">
+     <use style="fill:#1f77b4;stroke:#1f77b4;" x="379.811446" xlink:href="#C1_0_a2ed93b880" y="163.429156"/>
+    </g>
+   </g>
+   <g id="text_2">
+    <!-- $\psi &gt; 0$ -->
+    <defs>
+     <path d="M 65.90625 22.8125 
+C 67.796875 23.734375 67.90625 24.328125 67.90625 25.046875 
+C 67.90625 25.640625 67.703125 26.359375 65.90625 27.15625 
+L 11.796875 52.296875 
+C 10.5 52.890625 10.296875 53 10.09375 53 
+C 8.90625 53 8.203125 51.984375 8.203125 51.1875 
+C 8.203125 49.984375 9 49.578125 10.296875 48.96875 
+L 61.703125 25.046875 
+L 10.203125 1.015625 
+C 8.203125 0.109375 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.90625 11.796875 -2.3125 
+z
+" id="CMMI12-62"/>
+    </defs>
+    <g transform="translate(313.152965 26.4288)scale(0.168 -0.168)">
+     <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>
+ <defs>
+  <clipPath id="pea105021c5">
+   <rect height="315.36" width="207.72" x="7.2" y="7.2"/>
+  </clipPath>
+  <clipPath id="p0a8decf98d">
+   <rect height="315.36" width="207.72" x="230.04" y="7.2"/>
+  </clipPath>
+ </defs>
+</svg>
index b31ba20717ba74137f4b021cac75c4a04b2bcd59..8fda2693cee8063fccf11cf8de4dd48a35901c9b 100644 (file)
@@ -22,6 +22,7 @@
 
 #include <deal.II/dofs/dof_handler.h>
 
+#include <deal.II/fe/fe_interface_values.h>
 #include <deal.II/fe/fe_update_flags.h>
 #include <deal.II/fe/fe_values.h>
 
@@ -368,6 +369,379 @@ namespace NonMatching
     DiscreteQuadratureGenerator<dim> quadrature_generator;
   };
 
+
+
+  /**
+   * This class is intended to facilitate assembling interface terms on faces in
+   * immersed (in the sense of cut) finite element methods. These types of terms
+   * occur mainly in cut discontinuous Galerkin methods. This class works
+   * analogously to NonMatching::FEValues. The domain is assumed to be described
+   * by a level set function,
+   * $\psi : \mathbb{R}^{dim} \to \mathbb{R}$,
+   * and this class assumes that we want to integrate over two different regions
+   * of each face, $F$:
+   * @f[
+   * N = \{x \in F : \psi(x) < 0 \}, \\
+   * P = \{x \in F : \psi(x) > 0 \},
+   * @f]
+   * which we as before refer to as the "inside" and "outside" regions of the
+   * face.
+   * @image html non_matching_fe_interface_values.svg
+   * When calling one of the reinit functions of this class, this class
+   * will create quadrature rules over these regions in the background and set
+   * up dealii::FEInterfaceValues objects with these quadratures. These objects
+   * can then be accessed through one of the functions: get_inside_fe_values()
+   * or get_outside_fe_values(). For the same reason as described in
+   * NonMatching::FEValues, the returned the FEInterfaceValues objects that
+   * get_inside/outside_fe_values() returns is wrapped in a
+   * std_cxx17::optional. Assembling using this class would then typically be
+   * done in the following way:
+   *
+   * @code
+   *  NonMatching::FEInterfaceValues<dim> fe_interface_values(...)
+   *
+   *  for (const auto &cell : dof_handler.active_cell_iterators())
+   *    {
+   *      for (unsigned int face_index : cell->face_indices())
+   *        {
+   *          if (cell->at_boundary(face_index))
+   *            fe_interface_values.reinit(cell, face_index);
+   *          else
+   *            fe_interface_values.reinit(cell,
+   *                                       face_index,
+   *                                       subface_index,
+   *                                       cell->neighbor(face_index),
+   *                                       cell->neighbor_of_neighbor(face_index),
+   *                                       neighbor_subface_index);
+   *
+   *          const std_cxx17::optional<dealii::FEInterfaceValues<dim>>
+   *            &inside_fe_values = fe_interface_values.get_inside_fe_values();
+   *
+   *          // Check if the returned optional contains a value
+   *          if (inside_fe_values)
+   *            {
+   *              // Assemble locally
+   *            }
+   *        }
+   *    }
+   * @endcode
+   *
+   * To reduce the amount of work, the reinit() function of this class uses the
+   * MeshClassifier passed to the constructor to check how the incoming cell
+   * relates to the level set function. The immersed quadrature rules are only
+   * generated if the cell is intersected. If the cell is completely inside or
+   * outside, it returns a cached FEInterfaceValues object created with a
+   * quadrature over the reference cell: $[0, 1]^{dim-1}$.
+   */
+  template <int dim>
+  class FEInterfaceValues
+  {
+  public:
+    using AdditionalData =
+      typename FaceQuadratureGenerator<dim>::AdditionalData;
+
+    /**
+     * Constructor.
+     *
+     * @param fe_collection Collection of FiniteElements to be used.
+     * @param quadrature 1-dimensional quadrature rule used to generate the
+     * immersed quadrature rules. See the QuadratureGenerator class. On the non
+     * intersected faces a tensor product of this quadrature will be used.
+     * @param mesh_classifier Object used to determine when the immersed
+     * quadrature rules need to be generated.
+     * @param region_update_flags Struct storing UpdateFlags for the
+     * inside/outside region of the cell.
+     * @param dof_handler The DoFHandler associated with the discrete level set
+     * function.
+     * @param level_set The degrees of freedom of the discrete level set function.
+     * @param additional_data Additional data passed to the QuadratureGenerator
+     * class.
+     *
+     * @note  Pointers to @p fe_collection, @p mesh_classifier, @p dof_handler,
+     * and @p level_set are stored internally, so these need to have a longer life
+     * span than the instance of this class.
+     */
+    template <class VectorType>
+    FEInterfaceValues(const hp::FECollection<dim> &fe_collection,
+                      const Quadrature<1> &        quadrature,
+                      const RegionUpdateFlags      region_update_flags,
+                      const MeshClassifier<dim> &  mesh_classifier,
+                      const DoFHandler<dim> &      dof_handler,
+                      const VectorType &           level_set,
+                      const AdditionalData &additional_data = AdditionalData());
+
+    /**
+     * Constructor.
+     *
+     * @param mapping_collection Collection of Mappings to be used.
+     * @param fe_collection Collection of FiniteElements to be used.
+     * @param q_collection Collection of Quadrature rules over $[0, 1]^{dim-1}$
+     * that should be used when a face is not intersected and we do not need to
+     * generate immersed quadrature rules.
+     * @param q_collection_1D Collection of 1-dimensional quadrature rules used
+     * to generate the immersed quadrature rules. See the QuadratureGenerator
+     * class.
+     * @param mesh_classifier Object used to determine when the immersed
+     * quadrature rules need to be generated.
+     * @param region_update_flags Struct storing UpdateFlags for the
+     * inside/outside region of the cell.
+     * @param dof_handler The DoFHandler associated with the discrete level set
+     * function.
+     * @param level_set The degrees of freedom of the discrete level set function.
+     * @param additional_data Additional data passed to the QuadratureGenerator
+     * class.
+     *
+     * @note Pointers to @p mapping_collection, @p fe_collection,
+     * @p mesh_classifier, @p dof_handler, and @p level_set are stored
+     * internally, so these need to have a longer life span than the instance of
+     * this class.
+     */
+    template <class VectorType>
+    FEInterfaceValues(const hp::MappingCollection<dim> &mapping_collection,
+                      const hp::FECollection<dim> &     fe_collection,
+                      const hp::QCollection<dim - 1> &  q_collection,
+                      const hp::QCollection<1> &        q_collection_1D,
+                      const RegionUpdateFlags           region_update_flags,
+                      const MeshClassifier<dim> &       mesh_classifier,
+                      const DoFHandler<dim> &           dof_handler,
+                      const VectorType &                level_set,
+                      const AdditionalData &additional_data = AdditionalData());
+
+    /**
+     * Reinitialize on the shared face between two neighboring cells.
+     * After calling this function, an FEInterfaceValues object can be retrieved
+     * for each region using the functions get_inside_fe_values() or
+     * get_outside_fe_values().
+     *
+     * @note Currently @p cell and @p cell_neighbor need to have the same
+     * finite element associated with them.
+     *
+     * @note Specifying sub_face_no/sub_face_no_neighbor is currently not
+     * implemented, the passed values of these must equal
+     * numbers::invalid_unsigned_int.
+     */
+    template <class CellIteratorType, class CellNeighborIteratorType>
+    void
+    reinit(const CellIteratorType &        cell,
+           const unsigned int              face_no,
+           const unsigned int              sub_face_no,
+           const CellNeighborIteratorType &cell_neighbor,
+           const unsigned int              face_no_neighbor,
+           const unsigned int              sub_face_no_neighbor);
+
+
+    /**
+     * Reinitialize on a face on the boundary of the domain,
+     * when there is no neighbor and a single face @p face_no of
+     * the @p cell. After calling this function, an
+     * FEInterfaceValues object can be retrieved for each region using the
+     * functions get_inside_fe_values() or get_outside_fe_values().
+     */
+    template <class CellIteratorType>
+    void
+    reinit(const CellIteratorType &cell, const unsigned int face_no);
+
+    /**
+     * Return an dealii::FEInterfaceValues object reinitialized with a
+     * quadrature for the inside region of the cell:
+     * $\{x \in K : \psi(x) < 0 \}$.
+     *
+     * @note If the quadrature rule over the region is empty, e.g. because the
+     * cell is completely located in the outside domain, the returned
+     * optional will not contain a value.
+     */
+    const std_cxx17::optional<dealii::FEInterfaceValues<dim>> &
+    get_inside_fe_values() const;
+
+    /**
+     * Return an dealii::FEInterfaceValues object reinitialized with a
+     * quadrature for the outside region of the cell:
+     * $\{x \in K : \psi(x) > 0 \}$.
+     *
+     * @note If the quadrature rule over the region is empty, e.g. because the
+     * cell is completely located in the inside domain, the returned
+     * optional will not contain a value.
+     */
+    const std_cxx17::optional<dealii::FEInterfaceValues<dim>> &
+    get_outside_fe_values() const;
+
+  private:
+    /**
+     * Do work common to the constructors. The incoming QCollection should be
+     * quadratures integrating over $[0, 1]^{dim-1}$. These will be used on the
+     * non-intersected cells.
+     */
+    void
+    initialize(const hp::QCollection<dim - 1> &q_collection);
+
+    /**
+     * Do work that is common two the two reinit functions of the class.
+     * @p call_reinit is a std::function that describes how we should call
+     * reinit on a single dealii::FEInterfaceValues object, which is what
+     * differs between the two reinit functions.
+     */
+    template <bool level_dof_access>
+    void
+    do_reinit(
+      const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+      const unsigned int                                               face_no,
+      const std::function<void(dealii::FEInterfaceValues<dim> &)> &call_reinit);
+
+    /**
+     * A pointer to the collection of mappings to be used.
+     */
+    const SmartPointer<const hp::MappingCollection<dim>> mapping_collection;
+
+    /**
+     * A pointer to the collection of finite elements to be used.
+     */
+    const SmartPointer<const hp::FECollection<dim>> fe_collection;
+
+    /**
+     * Collection of 1-dimensional quadrature rules that are used by
+     * QuadratureGenerator as base for generating the immersed quadrature rules.
+     */
+    const hp::QCollection<1> q_collection_1D;
+
+    /**
+     * Location of the last face that reinit was called with.
+     */
+    LocationToLevelSet current_face_location;
+
+    /**
+     * Active fe index of the last cell that reinit was called with.
+     */
+    unsigned int active_fe_index;
+
+    /**
+     * The update flags passed to the constructor.
+     */
+    const RegionUpdateFlags region_update_flags;
+
+    /**
+     * Pointer to the MeshClassifier passed to the constructor.
+     */
+    const SmartPointer<const MeshClassifier<dim>> mesh_classifier;
+
+    /**
+     * For each element in the FECollection passed to the constructor,
+     * this object contains an dealii::FEInterfaceValues object created with a
+     * quadrature rule over the full reference cell: $[0, 1]^{dim-1}$ and
+     * UpdateFlags for the inside region. Thus, these optionals should always
+     * contain a value.
+     *
+     * When LocationToLevelSet of the cell is INSIDE (and we do not need
+     * to generate an immersed quadrature), we return the
+     * dealii::FEInterfaceValues object in this container corresponding to the
+     * cell's active_fe_index.
+     *
+     * This container is a std::deque, which is compatible with the
+     * `FEInterfaceValues` class that does not have a copy-constructor.
+     */
+    std::deque<std_cxx17::optional<dealii::FEInterfaceValues<dim>>>
+      fe_values_inside_full_quadrature;
+
+    /**
+     * For each element in the FECollection passed to the constructor,
+     * this object contains an dealii::FEInterfaceValues object created with a
+     * quadrature rule over the full reference cell: $[0, 1]^{dim-1}$ and
+     * UpdateFlags for the outside region. Thus, these optionals should always
+     * contain a value.
+     *
+     * When LocationToLevelSet of the cell is OUTSIDE (and we do not need
+     * to generate an immersed quadrature), we return the dealii::FEValues
+     * object in this container corresponding to the cell's active_fe_index.
+     *
+     * This container is a std::deque, which is compatible with the
+     * `FEInterfaceValues` class that does not have a copy-constructor.
+     */
+    std::deque<std_cxx17::optional<dealii::FEInterfaceValues<dim>>>
+      fe_values_outside_full_quadrature;
+
+    /**
+     * FEInterfaceValues object created with a quadrature rule integrating over
+     * the inside region, $\{x \in B : \psi(x) < 0 \}$, that was generated in
+     * the last call to reinit(..). If the cell in the last call was not
+     * intersected or if 0 quadrature points were generated, this optional will
+     * not contain a value.
+     */
+    std_cxx17::optional<dealii::FEInterfaceValues<dim>> fe_values_inside;
+
+    /**
+     * FEInterfaceValues object created with a quadrature rule integrating over
+     * the outside region, $\{x \in B : \psi(x) > 0 \}$, that was generated in
+     * the last call to reinit(..). If the cell in the last call was not
+     * intersected or if 0 quadrature points were generated, this optional will
+     * not contain a value.
+     */
+    std_cxx17::optional<dealii::FEInterfaceValues<dim>> fe_values_outside;
+
+    /**
+     * Object that generates the immersed quadrature rules.
+     */
+    DiscreteFaceQuadratureGenerator<dim> face_quadrature_generator;
+  };
+
+
+#ifndef DOXYGEN
+
+  /*---------------------- Inline functions ---------------------*/
+
+  template <int dim>
+  template <class CellIteratorType>
+  inline void
+  FEInterfaceValues<dim>::reinit(const CellIteratorType &cell,
+                                 const unsigned int      face_no)
+  {
+    // Lambda describing how we should call reinit on a single
+    // dealii::FEInterfaceValues object.
+    const auto reinit_operation =
+      [&cell, face_no](dealii::FEInterfaceValues<dim> &fe_interface_values) {
+        fe_interface_values.reinit(cell, face_no);
+      };
+
+    do_reinit(cell, face_no, reinit_operation);
+  }
+
+
+
+  template <int dim>
+  template <class CellIteratorType, class CellNeighborIteratorType>
+  inline void
+  FEInterfaceValues<dim>::reinit(const CellIteratorType &        cell,
+                                 const unsigned int              face_no,
+                                 const unsigned int              sub_face_no,
+                                 const CellNeighborIteratorType &cell_neighbor,
+                                 const unsigned int face_no_neighbor,
+                                 const unsigned int sub_face_no_neighbor)
+  {
+    Assert(sub_face_no == numbers::invalid_unsigned_int, ExcNotImplemented());
+    Assert(sub_face_no_neighbor == numbers::invalid_unsigned_int,
+           ExcNotImplemented());
+
+    // Lambda describing how we should call reinit on a single
+    // dealii::FEInterfaceValues object.
+    const auto reinit_operation =
+      [&cell,
+       face_no,
+       sub_face_no,
+       &cell_neighbor,
+       face_no_neighbor,
+       sub_face_no_neighbor](
+        dealii::FEInterfaceValues<dim> &fe_interface_values) {
+        fe_interface_values.reinit(cell,
+                                   face_no,
+                                   sub_face_no,
+                                   cell_neighbor,
+                                   face_no_neighbor,
+                                   sub_face_no_neighbor);
+      };
+
+    do_reinit(cell, face_no, reinit_operation);
+  }
+
+#endif
+
 } // namespace NonMatching
 DEAL_II_NAMESPACE_CLOSE
 
index 75d82d252669b4fcf694badb9bee145cc49445f5..fcbb399120f387451b1b3969a5507767a8a5372c 100644 (file)
@@ -268,6 +268,217 @@ namespace NonMatching
   }
 
 
+
+  template <int dim>
+  template <class VectorType>
+  FEInterfaceValues<dim>::FEInterfaceValues(
+    const hp::FECollection<dim> &fe_collection,
+    const Quadrature<1> &        quadrature,
+    const RegionUpdateFlags      region_update_flags,
+    const MeshClassifier<dim> &  mesh_classifier,
+    const DoFHandler<dim> &      dof_handler,
+    const VectorType &           level_set,
+    const AdditionalData &       additional_data)
+    : mapping_collection(&dealii::hp::StaticMappingQ1<dim>::mapping_collection)
+    , fe_collection(&fe_collection)
+    , q_collection_1D(quadrature)
+    , region_update_flags(region_update_flags)
+    , mesh_classifier(&mesh_classifier)
+    , face_quadrature_generator(q_collection_1D,
+                                dof_handler,
+                                level_set,
+                                additional_data)
+  {
+    // Tensor products of each quadrature in q_collection_1D. Used on the
+    // non-intersected cells.
+    hp::QCollection<dim - 1> q_collection;
+    for (unsigned int i = 0; i < q_collection_1D.size(); ++i)
+      q_collection.push_back(Quadrature<dim - 1>(q_collection_1D[i]));
+
+    initialize(q_collection);
+  }
+
+
+
+  template <int dim>
+  template <class VectorType>
+  FEInterfaceValues<dim>::FEInterfaceValues(
+    const hp::MappingCollection<dim> &mapping_collection,
+    const hp::FECollection<dim> &     fe_collection,
+    const hp::QCollection<dim - 1> &  q_collection,
+    const hp::QCollection<1> &        q_collection_1D,
+    const RegionUpdateFlags           region_update_flags,
+    const MeshClassifier<dim> &       mesh_classifier,
+    const DoFHandler<dim> &           dof_handler,
+    const VectorType &                level_set,
+    const AdditionalData &            additional_data)
+    : mapping_collection(&mapping_collection)
+    , fe_collection(&fe_collection)
+    , q_collection_1D(q_collection_1D)
+    , region_update_flags(region_update_flags)
+    , mesh_classifier(&mesh_classifier)
+    , face_quadrature_generator(q_collection_1D,
+                                dof_handler,
+                                level_set,
+                                additional_data)
+  {
+    initialize(q_collection);
+  }
+
+
+
+  template <int dim>
+  void
+  FEInterfaceValues<dim>::initialize(
+    const hp::QCollection<dim - 1> &q_collection)
+  {
+    current_face_location = LocationToLevelSet::unassigned;
+    active_fe_index       = numbers::invalid_unsigned_int;
+
+    Assert(fe_collection->size() > 0,
+           ExcMessage("Incoming hp::FECollection can not be empty."));
+    Assert(
+      mapping_collection->size() == fe_collection->size() ||
+        mapping_collection->size() == 1,
+      ExcMessage(
+        "Size of hp::MappingCollection must be the same as hp::FECollection or 1."));
+    Assert(
+      q_collection.size() == fe_collection->size() || q_collection.size() == 1,
+      ExcMessage(
+        "Size of hp::QCollection<dim> must be the same as hp::FECollection or 1."));
+    Assert(
+      q_collection_1D.size() == fe_collection->size() ||
+        q_collection_1D.size() == 1,
+      ExcMessage(
+        "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
+
+    // For each element in fe_collection, create dealii::FEInterfaceValues
+    // objects to use on the non-intersected cells.
+    fe_values_inside_full_quadrature.resize(fe_collection->size());
+    fe_values_outside_full_quadrature.resize(fe_collection->size());
+    for (unsigned int fe_index = 0; fe_index < fe_collection->size();
+         ++fe_index)
+      {
+        const unsigned int mapping_index =
+          mapping_collection->size() > 1 ? fe_index : 0;
+        const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0;
+
+        fe_values_inside_full_quadrature[fe_index].emplace(
+          (*mapping_collection)[mapping_index],
+          (*fe_collection)[fe_index],
+          q_collection[q_index],
+          region_update_flags.inside);
+        fe_values_outside_full_quadrature[fe_index].emplace(
+          (*mapping_collection)[mapping_index],
+          (*fe_collection)[fe_index],
+          q_collection[q_index],
+          region_update_flags.outside);
+      }
+  }
+
+
+
+  template <int dim>
+  template <bool level_dof_access>
+  void
+  FEInterfaceValues<dim>::do_reinit(
+    const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+    const unsigned int                                               face_no,
+    const std::function<void(dealii::FEInterfaceValues<dim> &)> &call_reinit)
+  {
+    current_face_location =
+      mesh_classifier->location_to_level_set(cell, face_no);
+    active_fe_index = cell->active_fe_index();
+
+    // These objects were created with a quadrature based on the previous cell
+    // and are thus no longer valid.
+    fe_values_inside.reset();
+    fe_values_outside.reset();
+
+    switch (current_face_location)
+      {
+        case LocationToLevelSet::inside:
+          {
+            call_reinit(*fe_values_inside_full_quadrature.at(active_fe_index));
+            break;
+          }
+        case LocationToLevelSet::outside:
+          {
+            call_reinit(*fe_values_outside_full_quadrature.at(active_fe_index));
+            break;
+          }
+        case LocationToLevelSet::intersected:
+          {
+            const unsigned int mapping_index =
+              mapping_collection->size() > 1 ? active_fe_index : 0;
+            const unsigned int q1D_index =
+              q_collection_1D.size() > 1 ? active_fe_index : 0;
+
+            face_quadrature_generator.set_1D_quadrature(q1D_index);
+            face_quadrature_generator.generate(cell, face_no);
+
+            const Quadrature<dim - 1> &inside_quadrature =
+              face_quadrature_generator.get_inside_quadrature();
+            const Quadrature<dim - 1> &outside_quadrature =
+              face_quadrature_generator.get_outside_quadrature();
+
+            // Even if a cell is formally intersected the number of created
+            // quadrature points can be 0. Avoid creating an FEInterfaceValues
+            // object if that is the case.
+            if (inside_quadrature.size() > 0)
+              {
+                fe_values_inside.emplace((*mapping_collection)[mapping_index],
+                                         (*fe_collection)[active_fe_index],
+                                         inside_quadrature,
+                                         region_update_flags.inside);
+
+                call_reinit(*fe_values_inside);
+              }
+
+            if (outside_quadrature.size() > 0)
+              {
+                fe_values_outside.emplace((*mapping_collection)[mapping_index],
+                                          (*fe_collection)[active_fe_index],
+                                          outside_quadrature,
+                                          region_update_flags.outside);
+
+                call_reinit(*fe_values_outside);
+              }
+            break;
+          }
+        default:
+          {
+            Assert(false, ExcInternalError());
+            break;
+          }
+      }
+  }
+
+
+
+  template <int dim>
+  const std_cxx17::optional<dealii::FEInterfaceValues<dim>> &
+  FEInterfaceValues<dim>::get_inside_fe_values() const
+  {
+    if (current_face_location == LocationToLevelSet::inside)
+      return fe_values_inside_full_quadrature.at(active_fe_index);
+    else
+      return fe_values_inside;
+  }
+
+
+
+  template <int dim>
+  const std_cxx17::optional<dealii::FEInterfaceValues<dim>> &
+  FEInterfaceValues<dim>::get_outside_fe_values() const
+  {
+    if (current_face_location == LocationToLevelSet::outside)
+      return fe_values_outside_full_quadrature.at(active_fe_index);
+    else
+      return fe_values_outside;
+  }
+
+
 #include "fe_values.inst"
 
 } // namespace NonMatching
index e9cdf571c47bc521ecf8cb79e84e64e30dd6d395..e8bccf8e596fc18381e74a15b496a643e3a20400 100644 (file)
@@ -17,6 +17,8 @@
 for (deal_II_dimension : DIMENSIONS)
   {
     template class FEValues<deal_II_dimension>;
+
+    template class FEInterfaceValues<deal_II_dimension>;
   }
 
 for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
@@ -41,12 +43,39 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
       const DoFHandler<deal_II_dimension> &,
       const VEC &,
       const typename FEValues<deal_II_dimension>::AdditionalData &);
+
+    template FEInterfaceValues<deal_II_dimension>::FEInterfaceValues(
+      const hp::FECollection<deal_II_dimension> &,
+      const Quadrature<1> &,
+      const RegionUpdateFlags,
+      const MeshClassifier<deal_II_dimension> &,
+      const DoFHandler<deal_II_dimension> &,
+      const VEC &,
+      const FEInterfaceValues<deal_II_dimension>::AdditionalData &);
+
+    template FEInterfaceValues<deal_II_dimension>::FEInterfaceValues(
+      const hp::MappingCollection<deal_II_dimension> &,
+      const hp::FECollection<deal_II_dimension> &,
+      const hp::QCollection<deal_II_dimension - 1> &,
+      const hp::QCollection<1> &,
+      const RegionUpdateFlags,
+      const MeshClassifier<deal_II_dimension> &,
+      const DoFHandler<deal_II_dimension> &,
+      const VEC &,
+      const AdditionalData &);
   }
 
-// Template reinit function
+// Template reinit functions
 for (deal_II_dimension : DIMENSIONS; lda : BOOL)
   {
     template void FEValues<deal_II_dimension>::reinit(
       const TriaIterator<
         DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &);
+
+    template void FEInterfaceValues<deal_II_dimension>::do_reinit(
+      const TriaIterator<
+        DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &,
+      const unsigned int,
+      const std::function<void(dealii::FEInterfaceValues<deal_II_dimension> &)>
+        &);
   }
diff --git a/tests/non_matching/fe_interface_values.cc b/tests/non_matching/fe_interface_values.cc
new file mode 100644 (file)
index 0000000..6677f11
--- /dev/null
@@ -0,0 +1,227 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function_signed_distance.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx17/optional.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_cartesian.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/non_matching/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+// Set up a triangulation with 2 elements in a row: |-0-|-1-|, over
+// [0,2]x[0,1]^{dim-1}, and a level set  function with a cut in the plane
+// $x_[dim-1] = 0.5$,  such that it cuts through the middle of both cells in 2D
+// and 3D and thought the middle of cell 0 in 1D. Test that we can create an
+// object of type NonMatching::FEInterfaceValues. Call reinit on all the faces
+// of cell 0 and check that the dealii::FEInterfaceValues objects we get from
+// get_inside/outside_fe_values are set up correctly. That is, the
+// std::optionals that should contain a value do and the ones that should not
+// contain a value do not.
+template <int dim>
+class Test
+{
+public:
+  Test();
+
+  void
+  run();
+
+private:
+  void
+  setup_mesh();
+
+  void
+  setup_discrete_level_set();
+
+  void
+  print_which_optionals_have_values_on_cell_0(
+    NonMatching::FEInterfaceValues<dim> &fe_values);
+
+  Triangulation<dim>    triangulation;
+  hp::FECollection<dim> fe_collection;
+  DoFHandler<dim>       dof_handler;
+
+  hp::MappingCollection<dim>           mapping_collection;
+  hp::QCollection<dim - 1>             q_collection;
+  hp::QCollection<1>                   q_collection1D;
+  const NonMatching::RegionUpdateFlags region_update_flags;
+
+  Vector<double>                   level_set;
+  NonMatching::MeshClassifier<dim> mesh_classifier;
+};
+
+
+
+template <int dim>
+Test<dim>::Test()
+  : dof_handler(triangulation)
+  , mesh_classifier(dof_handler, level_set)
+{
+  fe_collection.push_back(FE_Q<dim>(1));
+  mapping_collection.push_back(MappingCartesian<dim>());
+
+  const unsigned int n_quadrature_points = 1;
+  q_collection.push_back(QGauss<dim - 1>(n_quadrature_points));
+  q_collection1D.push_back(QGauss<1>(n_quadrature_points));
+}
+
+
+
+template <int dim>
+void
+Test<dim>::run()
+{
+  setup_mesh();
+  setup_discrete_level_set();
+  mesh_classifier.reclassify();
+
+  // Test to create an object with each constructor.
+  {
+    deallog << "Simple constructor:" << std::endl;
+    NonMatching::FEInterfaceValues<dim> fe_values(fe_collection,
+                                                  q_collection1D[0],
+                                                  region_update_flags,
+                                                  mesh_classifier,
+                                                  dof_handler,
+                                                  level_set);
+    print_which_optionals_have_values_on_cell_0(fe_values);
+  }
+  {
+    deallog << "Advanced constructor:" << std::endl;
+    NonMatching::FEInterfaceValues<dim> fe_values(mapping_collection,
+                                                  fe_collection,
+                                                  q_collection,
+                                                  q_collection1D,
+                                                  region_update_flags,
+                                                  mesh_classifier,
+                                                  dof_handler,
+                                                  level_set);
+    print_which_optionals_have_values_on_cell_0(fe_values);
+  }
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_mesh()
+{
+  const Point<dim> lower_left;
+  Point<dim>       upper_right;
+
+  std::vector<unsigned int> repetitions;
+  upper_right[0] = 2;
+  repetitions.push_back(2);
+  for (unsigned int d = 1; d < dim; ++d)
+    {
+      upper_right[d] = 1;
+      repetitions.push_back(1);
+    }
+
+  GridGenerator::subdivided_hyper_rectangle(triangulation,
+                                            repetitions,
+                                            lower_left,
+                                            upper_right);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_discrete_level_set()
+{
+  dof_handler.distribute_dofs(fe_collection);
+  level_set.reinit(dof_handler.n_dofs());
+
+  Point<dim> point_on_zero_contour;
+  point_on_zero_contour[dim - 1] = 0.5;
+  const Functions::SignedDistance::Plane<dim> analytical_levelset(
+    point_on_zero_contour, Point<dim>::unit_vector(dim - 1));
+
+  VectorTools::interpolate(dof_handler, analytical_levelset, level_set);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::print_which_optionals_have_values_on_cell_0(
+  NonMatching::FEInterfaceValues<dim> &fe_values)
+{
+  const auto cell = dof_handler.begin_active();
+
+  for (unsigned int face_index : cell->face_indices())
+    {
+      deallog << "face " << face_index << std::endl;
+
+      if (cell->at_boundary(face_index))
+        {
+          fe_values.reinit(cell, face_index);
+        }
+      else
+        {
+          const unsigned int invalid_subface =
+            dealii::numbers::invalid_unsigned_int;
+
+          fe_values.reinit(cell,
+                           face_index,
+                           invalid_subface,
+                           cell->neighbor(face_index),
+                           cell->neighbor_of_neighbor(face_index),
+                           invalid_subface);
+        }
+
+      if (fe_values.get_inside_fe_values())
+        deallog << "inside has value" << std::endl;
+
+      if (fe_values.get_outside_fe_values())
+        deallog << "outside has value" << std::endl;
+    }
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+  deallog << "dim = " << dim << std::endl;
+  Test<dim> test;
+  test.run();
+  deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  run_test<1>();
+  run_test<2>();
+  run_test<3>();
+}
diff --git a/tests/non_matching/fe_interface_values.output b/tests/non_matching/fe_interface_values.output
new file mode 100644 (file)
index 0000000..52bca2f
--- /dev/null
@@ -0,0 +1,73 @@
+
+DEAL::dim = 1
+DEAL::Simple constructor:
+DEAL::face 0
+DEAL::inside has value
+DEAL::face 1
+DEAL::outside has value
+DEAL::Advanced constructor:
+DEAL::face 0
+DEAL::inside has value
+DEAL::face 1
+DEAL::outside has value
+DEAL::
+DEAL::dim = 2
+DEAL::Simple constructor:
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::face 3
+DEAL::outside has value
+DEAL::Advanced constructor:
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::face 3
+DEAL::outside has value
+DEAL::
+DEAL::dim = 3
+DEAL::Simple constructor:
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 3
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 4
+DEAL::inside has value
+DEAL::face 5
+DEAL::outside has value
+DEAL::Advanced constructor:
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 3
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 4
+DEAL::inside has value
+DEAL::face 5
+DEAL::outside has value
+DEAL::
diff --git a/tests/non_matching/fe_interface_values_non_standard_meshes.cc b/tests/non_matching/fe_interface_values_non_standard_meshes.cc
new file mode 100644 (file)
index 0000000..fcfa25c
--- /dev/null
@@ -0,0 +1,206 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// For all combinations of 2-cell triangulations that can be created by
+// GridGenerator::non_standard_orientation_mesh, set up a discrete level set
+// function corresponding to
+// (y  - 0.25) + (z - 0.25) = 0
+// and use NonMatching::FEInterfaceValues to generate the face quadratures on
+// the shared face at x = 0. Check that the quadrature points on both sides of
+// the face are the same.
+
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function_signed_distance.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx17/optional.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/non_matching/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <bitset>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class Test
+{
+public:
+  Test(const Triangulation<dim> &triangulation);
+
+  void
+  run();
+
+private:
+  void
+  setup_discrete_level_set();
+
+  /*
+   * For the incoming FEInterfaceValues object, print the distance between the
+   * real space quadrature points of the underlying FEFaceValues objects, for
+   * all quadrature points.
+   */
+  void
+  print_quadrature_point_difference_on_both_sides(
+    const dealii::FEInterfaceValues<dim> &fe_interface_values);
+
+  hp::FECollection<dim> fe_collection;
+  DoFHandler<dim>       dof_handler;
+  Vector<double>        level_set;
+
+  NonMatching::MeshClassifier<dim> mesh_classifier;
+
+  QGauss<1> quadrature_1D;
+};
+
+
+
+template <int dim>
+Test<dim>::Test(const Triangulation<dim> &triangulation)
+  : dof_handler(triangulation)
+  , mesh_classifier(dof_handler, level_set)
+  , quadrature_1D(2)
+{
+  fe_collection.push_back(FE_Q<dim>(1));
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_discrete_level_set()
+{
+  dof_handler.distribute_dofs(fe_collection);
+  level_set.reinit(dof_handler.n_dofs());
+
+  Point<dim>     point_on_zero_contour;
+  Tensor<1, dim> plane_normal;
+  for (unsigned int i = 1; i < dim; i++)
+    {
+      point_on_zero_contour(i) = 0.25;
+      plane_normal[i]          = 1;
+    }
+
+  const Functions::SignedDistance::Plane<dim> analytical_levelset(
+    point_on_zero_contour, plane_normal);
+
+  VectorTools::interpolate(dof_handler, analytical_levelset, level_set);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::print_quadrature_point_difference_on_both_sides(
+  const dealii::FEInterfaceValues<dim> &fe_interface_values)
+{
+  for (unsigned int q : fe_interface_values.quadrature_point_indices())
+    {
+      const Point<dim> &point_on_side_0 =
+        fe_interface_values.get_fe_face_values(0).quadrature_point(q);
+      const Point<dim> &point_on_side_1 =
+        fe_interface_values.get_fe_face_values(1).quadrature_point(q);
+
+      deallog << point_on_side_0.distance(point_on_side_1) << std::endl;
+    }
+}
+
+
+
+template <int dim>
+void
+Test<dim>::run()
+{
+  setup_discrete_level_set();
+  mesh_classifier.reclassify();
+
+  NonMatching::RegionUpdateFlags region_update_flags;
+  region_update_flags.inside  = update_quadrature_points;
+  region_update_flags.outside = update_quadrature_points;
+
+  NonMatching::FEInterfaceValues<dim> fe_values(fe_collection,
+                                                quadrature_1D,
+                                                region_update_flags,
+                                                mesh_classifier,
+                                                dof_handler,
+                                                level_set);
+
+  const auto cell = dof_handler.begin_active();
+
+  // Test that the points agree only on the interface between the two cells,
+  // which face this is depends on the incoming triangulation.
+  for (unsigned int face : cell->face_indices())
+    if (!cell->face(face)->at_boundary())
+      {
+        const unsigned int invalid_subface =
+          dealii::numbers::invalid_unsigned_int;
+
+        fe_values.reinit(cell,
+                         face,
+                         invalid_subface,
+                         cell->neighbor(face),
+                         cell->neighbor_of_neighbor(face),
+                         invalid_subface);
+
+        Assert(fe_values.get_inside_fe_values(), ExcInternalError());
+        Assert(fe_values.get_outside_fe_values(), ExcInternalError());
+
+        print_quadrature_point_difference_on_both_sides(
+          *fe_values.get_inside_fe_values());
+        print_quadrature_point_difference_on_both_sides(
+          *fe_values.get_outside_fe_values());
+      }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  const int dim = 3;
+
+  // non_standard_orientation_mesh takes 4 boolean parameters. Test all
+  // combinations.
+  const unsigned int n_arguments    = 4;
+  const unsigned int n_combinations = std::pow(2, 4);
+  for (unsigned int i = 0; i < n_combinations; i++)
+    {
+      const std::bitset<n_arguments> bits(i);
+
+      const bool face_orientation     = bits[0];
+      const bool face_flip            = bits[1];
+      const bool face_rotation        = bits[2];
+      const bool manipulate_left_cube = bits[3];
+
+      Triangulation<dim> triangulation;
+      GridGenerator::non_standard_orientation_mesh(triangulation,
+                                                   face_orientation,
+                                                   face_flip,
+                                                   face_rotation,
+                                                   manipulate_left_cube);
+      Test<dim> test(triangulation);
+      test.run();
+    }
+}
diff --git a/tests/non_matching/fe_interface_values_non_standard_meshes.output b/tests/non_matching/fe_interface_values_non_standard_meshes.output
new file mode 100644 (file)
index 0000000..f8b0136
--- /dev/null
@@ -0,0 +1,193 @@
+
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0

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.