--- /dev/null
+<?xml version="1.0" encoding="utf-8" standalone="no"?>
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
+ "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<!-- Created with matplotlib (https://matplotlib.org/) -->
+<svg height="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 < 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 > 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>
#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>
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
}
+
+ 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
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)
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> &)>
+ &);
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+ }
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff
