]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a CutFEM tutorial as step-85
authorSimon Sticko <simon@sticko.se>
Sun, 2 Jan 2022 10:27:50 +0000 (11:27 +0100)
committerSimon Sticko <simon@sticko.se>
Sat, 29 Jan 2022 13:25:33 +0000 (14:25 +0100)
This step shows how the NonMatching classes: MeshClassifier,
FEImmersedSurfaceValues and NonMatching::FEValues are intended to work
together. This is done by solving the Poisson equation on a hypercube
mesh where the domain is described by a level set function

14 files changed:
doc/doxygen/images/step-85-active-mesh.svg [new file with mode: 0644]
doc/doxygen/images/step-85-background-mesh.svg [new file with mode: 0644]
doc/doxygen/images/step-85-ghost-faces.svg [new file with mode: 0644]
doc/doxygen/images/step-85-level-set.png [new file with mode: 0644]
doc/doxygen/images/step-85-solution.png [new file with mode: 0644]
doc/doxygen/references.bib
doc/doxygen/tutorial/tutorial.h.in
examples/step-85/CMakeLists.txt [new file with mode: 0644]
examples/step-85/doc/builds-on [new file with mode: 0644]
examples/step-85/doc/intro.dox [new file with mode: 0644]
examples/step-85/doc/kind [new file with mode: 0644]
examples/step-85/doc/results.dox [new file with mode: 0644]
examples/step-85/doc/tooltip [new file with mode: 0644]
examples/step-85/step-85.cc [new file with mode: 0644]

diff --git a/doc/doxygen/images/step-85-active-mesh.svg b/doc/doxygen/images/step-85-active-mesh.svg
new file mode 100644 (file)
index 0000000..05f00ca
--- /dev/null
@@ -0,0 +1,377 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+   xmlns:dc="http://purl.org/dc/elements/1.1/"
+   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+   xmlns:svg="http://www.w3.org/2000/svg"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   width="57.979237mm"
+   height="68.900131mm"
+   viewBox="0 0 57.979238 68.900131"
+   version="1.1"
+   id="svg8"
+   inkscape:version="0.92.5 (2060ec1f9f, 2020-04-08)"
+   sodipodi:docname="step-85-submesh.svg">
+  <defs
+     id="defs2" />
+  <sodipodi:namedview
+     id="base"
+     pagecolor="#ffffff"
+     bordercolor="#666666"
+     borderopacity="1.0"
+     inkscape:pageopacity="0.0"
+     inkscape:pageshadow="2"
+     inkscape:zoom="1.979899"
+     inkscape:cx="179.60987"
+     inkscape:cy="191.3565"
+     inkscape:document-units="mm"
+     inkscape:current-layer="layer1"
+     showgrid="false"
+     inkscape:window-width="1366"
+     inkscape:window-height="713"
+     inkscape:window-x="0"
+     inkscape:window-y="0"
+     inkscape:window-maximized="1"
+     inkscape:snap-global="false"
+     fit-margin-top="0"
+     fit-margin-left="0"
+     fit-margin-right="0"
+     fit-margin-bottom="0"
+     showguides="false" />
+  <metadata
+     id="metadata5">
+    <rdf:RDF>
+      <cc:Work
+         rdf:about="">
+        <dc:format>image/svg+xml</dc:format>
+        <dc:type
+           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
+        <dc:title></dc:title>
+      </cc:Work>
+    </rdf:RDF>
+  </metadata>
+  <g
+     inkscape:label="Layer 1"
+     inkscape:groupmode="layer"
+     id="layer1"
+     transform="translate(334.16873,-203.02181)">
+    <g
+       style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;letter-spacing:normal;word-spacing:normal;text-anchor:start;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10.43299961;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+       id="g1099"
+       transform="matrix(0.13432553,0,0,-0.13432553,-370.4005,324.63128)"
+       xml:space="preserve"
+       stroke-miterlimit="10.433"
+       font-style="normal"
+       font-variant="normal"
+       font-weight="normal"
+       font-stretch="normal"
+       font-size-adjust="none"
+       letter-spacing="normal"
+       word-spacing="normal"><path
+         style="fill:#000000;stroke-width:0"
+         inkscape:connector-curvature="0"
+         d="M 479.74,437.14 H 488 c 6.58,0 9.07,3.78 9.07,4.43 0,0.25 -0.19,0.35 -0.44,0.35 -0.56,0 -1.84,-0.5 -3.09,-1.39 -0.41,-0.3 -0.5,-0.3 -1.46,-0.3 h -21.76 c -6.18,0 -13.11,-5.33 -13.11,-9.42 0,-0.3 0,-0.6 0.4,-0.6 1.19,0 4.33,1.65 4.44,2.85 0.05,0.84 0.14,1.78 1.14,3.08 0.5,0.61 0.89,1 3.1,1 h 9.9 c -0.84,-0.55 -1.05,-1.3 -1.14,-1.8 l -5.69,-22.56 c -0.59,-2.44 -1.64,-5.58 -3.07,-8.82 -0.3,-0.6 -0.3,-0.7 -0.3,-0.81 0,-0.29 0.19,-0.34 0.44,-0.34 1,0 3.43,1.34 4.18,2.39 0.3,0.5 2.5,5.58 3.55,9.72 z"
+         id="path1093" /><path
+         style="fill:#000000;stroke-width:0"
+         inkscape:connector-curvature="0"
+         d="m 506.27,447.34 v 0 0.01 0 0.01 0 l 0.01,0.01 v 0 0.01 0.01 h 0.01 v 0.02 l 0.01,0.02 v 0.02 l 0.01,0.02 v 0.02 l 0.01,0.02 v 0.02 l 0.01,0.02 0.01,0.03 0.01,0.04 0.01,0.05 v 0.02 l 0.01,0.03 v 0.02 l 0.01,0.02 v 0.02 l 0.01,0.02 v 0.01 l 0.01,0.02 v 0.01 0 0.01 0.01 0 0.01 0 0.01 h 0.01 v 0 0.01 0 c 0,0.18 -0.14,0.5 -0.57,0.5 -0.68,0 -3.57,-0.28 -4.45,-0.36 -0.28,-0.03 -0.77,-0.07 -0.77,-0.8 0,-0.48 0.49,-0.48 0.91,-0.48 1.67,0 1.67,-0.25 1.67,-0.54 0,-0.23 -0.08,-0.45 -0.14,-0.76 l -4.92,-19.72 c -0.17,-0.62 -0.17,-0.69 -0.17,-0.77 0,-0.51 0.42,-1.1 1.19,-1.1 0.39,0 1.04,0.17 1.43,0.9 0.1,0.2 0.41,1.45 0.58,2.19 l 0.81,3.11 c 0.1,0.51 0.46,1.84 0.55,2.36 0.36,1.33 0.36,1.36 1.05,2.48 1.12,1.71 2.86,3.69 5.58,3.69 1.95,0 2.06,-1.61 2.06,-2.44 0,-2.09 -1.5,-5.95 -2.06,-7.42 -0.38,-0.98 -0.52,-1.3 -0.52,-1.89 0,-1.84 1.53,-2.98 3.3,-2.98 3.5,0 5.03,4.79 5.03,5.32 0,0.46 -0.45,0.46 -0.56,0.46 -0.49,0 -0.53,-0.21 -0.66,-0.6 -0.81,-2.78 -2.34,-4.22 -3.7,-4.22 -0.74,0 -0.88,0.49 -0.88,1.22 0,0.8 0.19,1.25 0.81,2.83 0.41,1.08 1.85,4.77 1.85,6.72 0,0.56 0,2.03 -1.3,3.03 -0.59,0.45 -1.59,0.94 -3.23,0.94 -2.55,0 -4.39,-1.39 -5.68,-2.89 z"
+         id="path1095" /><path
+         style="fill:#000000;stroke-width:0"
+         inkscape:connector-curvature="0"
+         d="m 509.02,398.26 h -1.19 c -0.03,-0.11 -0.35,-1.75 -0.63,-2.55 -0.2,-0.56 -0.31,-0.84 -2.47,-0.84 h -2.96 c 0.59,1.95 1.61,3.31 3.56,5.92 1.7,2.31 3.23,4.61 3.23,7.42 0,4.7 -5.04,8.69 -11.73,8.69 -6.71,0 -11.75,-3.99 -11.75,-8.69 0,-3.06 1.67,-5.33 4.11,-8.59 1.01,-1.33 2.12,-2.89 2.69,-4.75 h -2.97 c -2.08,0 -2.22,0.28 -2.39,0.67 -0.29,0.66 -0.71,2.64 -0.71,2.72 h -1.19 l 1.26,-5.86 h 6.28 c 0.73,0 0.96,0 0.96,0.76 0,2.52 -1.28,5.47 -2.01,7.18 -1.33,3.07 -2.3,5.34 -2.3,7.92 0,4.94 3.94,7.65 8.02,7.65 4.11,0 8.01,-2.75 8.01,-7.65 0,-2.63 -1.04,-5.02 -2.23,-7.74 -0.77,-1.75 -2.09,-4.81 -2.09,-7.32 0,-0.8 0.2,-0.8 1.01,-0.8 h 6.24 z"
+         id="path1097" /></g>    <path
+       inkscape:connector-curvature="0"
+       id="path74-6"
+       style="fill:none;stroke:#0000ed;stroke-width:0.37742284;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -279.07075,230.2073 c 0,14.30493 -11.5962,25.90111 -25.90113,25.90111 -14.3045,0 -25.90113,-11.59618 -25.90113,-25.90111 0,-14.3049 11.59663,-25.9011 25.90113,-25.9011 14.30493,0 25.90113,11.5962 25.90113,25.9011 z" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path14"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -334.01776,217.63504 h 7.06046 m 0,0 v 7.06049 m 0,0 h -7.06046 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path16"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -334.01776,224.86563 h 7.06046 m 0,0 v 7.06049 m 0,0 h -7.06046 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path18"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -334.01776,232.09676 h 7.06046 m 0,0 v 7.06049 m 0,0 h -7.06046 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path20"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -334.01776,239.32789 h 7.06046 m 0,0 v 7.06049 m 0,0 h -7.06046 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path24"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -326.78666,210.40391 h 7.0605 m 0,0 v 7.06049 m 0,0 h -7.0605 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path26"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -326.78666,217.63504 h 7.0605 m 0,0 v 7.06049 m 0,0 h -7.0605 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path28"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -326.78666,224.86563 h 7.0605 m 0,0 v 7.06049 m 0,0 h -7.0605 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path30"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -326.78666,232.09676 h 7.0605 m 0,0 v 7.06049 m 0,0 h -7.0605 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path32"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -326.78666,239.32789 h 7.0605 m 0,0 v 7.06049 m 0,0 h -7.0605 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path34"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -326.78666,246.55891 h 7.0605 m 0,0 v 7.06059 m 0,0 h -7.0605 m 0,0 v -7.06059" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path36"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -283.57159,210.40389 h 7.0605 m 0,0 v 7.06051 m 0,0 h -7.0605 m 0,0 v -7.06051" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path38"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,203.17278 h 7.06048 m 0,0 v 7.06049 m 0,0 h -7.06048 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path40"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,210.40391 h 7.06048 m 0,0 v 7.06049 m 0,0 h -7.06048 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path42"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,217.63504 h 7.06048 m 0,0 v 7.06049 m 0,0 h -7.06048 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path44"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,224.86563 h 7.06048 m 0,0 v 7.06049 m 0,0 h -7.06048 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path46"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,232.09676 h 7.06048 m 0,0 v 7.06049 m 0,0 h -7.06048 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path48"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,239.32789 h 7.06048 m 0,0 v 7.06049 m 0,0 h -7.06048 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path50"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,246.55891 h 7.06048 m 0,0 v 7.06059 m 0,0 h -7.06048 m 0,0 v -7.06059" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path52"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -319.55605,253.78959 h 7.06048 m 0,0 v 7.06051 m 0,0 h -7.06048 m 0,0 v -7.06051" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path54"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,203.17278 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path56"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,210.40391 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path58"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,217.63504 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path60"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,224.86563 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path62"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,232.09676 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path64"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,239.32789 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path66"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,246.55891 h 7.06049 m 0,0 v 7.06059 m 0,0 h -7.06049 m 0,0 v -7.06059" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path68"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -312.32493,253.78959 h 7.06049 m 0,0 v 7.06051 m 0,0 h -7.06049 m 0,0 v -7.06051" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path70"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,203.17278 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path72"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,210.40391 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path74"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,217.63504 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path76"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,224.86563 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path78"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,232.09676 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path80"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,239.32789 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path82"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,246.55891 h 7.06049 m 0,0 v 7.06059 m 0,0 h -7.06049 m 0,0 v -7.06059" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path84"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -305.0938,253.78959 h 7.06049 m 0,0 v 7.06051 m 0,0 h -7.06049 m 0,0 v -7.06051" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path86"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,203.17278 h 7.0606 m 0,0 v 7.06049 m 0,0 h -7.0606 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path88"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,210.40391 h 7.0606 m 0,0 v 7.06049 m 0,0 h -7.0606 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path90"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,217.63504 h 7.0606 m 0,0 v 7.06049 m 0,0 h -7.0606 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path92"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,224.86563 h 7.0606 m 0,0 v 7.06049 m 0,0 h -7.0606 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path94"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,232.09676 h 7.0606 m 0,0 v 7.06049 m 0,0 h -7.0606 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path96"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,239.32789 h 7.0606 m 0,0 v 7.06049 m 0,0 h -7.0606 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path98"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,246.55891 h 7.0606 m 0,0 v 7.06059 m 0,0 h -7.0606 m 0,0 v -7.06059" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path100"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -297.86278,253.78959 h 7.0606 m 0,0 v 7.06051 m 0,0 h -7.0606 m 0,0 v -7.06051" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path102"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -290.63208,210.40391 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path104"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -290.63208,217.63504 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path106"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -290.63208,224.86563 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path108"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -290.63208,232.09676 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path110"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -290.63208,239.32789 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path112"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -290.63208,246.55891 h 7.06049 m 0,0 v 7.06059 m 0,0 h -7.06049 m 0,0 v -7.06059" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path116"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -283.40095,217.63504 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path118"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -283.40095,224.86563 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path120"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -283.40095,232.09676 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path122"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -283.40095,239.32789 h 7.06049 m 0,0 v 7.06049 m 0,0 h -7.06049 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path36-3"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -290.63209,203.34338 h 7.0605 m 0,0 v 7.06051 m 0,0 h -7.0605 m 0,0 v -7.06051" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path14-7"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -333.96524,210.40391 h 7.06046 m 0,0 v 7.06049 m 0,0 h -7.06046 m 0,0 v -7.06049" />
+    <path
+       inkscape:connector-curvature="0"
+       id="path38-5"
+       style="fill:none;stroke:#000000;stroke-width:0.30193827;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+       d="m -326.82914,203.17278 h 7.06048 m 0,0 v 7.06049 m 0,0 h -7.06048 m 0,0 v -7.06049" />
+  </g>
+  <g
+     inkscape:groupmode="layer"
+     id="layer2"
+     inkscape:label="Layer 2"
+     transform="translate(334.16873,-203.02181)" />
+</svg>
diff --git a/doc/doxygen/images/step-85-background-mesh.svg b/doc/doxygen/images/step-85-background-mesh.svg
new file mode 100644 (file)
index 0000000..8387329
--- /dev/null
@@ -0,0 +1,411 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+   xmlns:dc="http://purl.org/dc/elements/1.1/"
+   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+   xmlns:svg="http://www.w3.org/2000/svg"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   width="57.979233mm"
+   height="67.967903mm"
+   viewBox="0 0 57.979234 67.967903"
+   version="1.1"
+   id="svg8"
+   inkscape:version="0.92.5 (2060ec1f9f, 2020-04-08)"
+   sodipodi:docname="step-85-background-mesh.svg">
+  <defs
+     id="defs2" />
+  <sodipodi:namedview
+     id="base"
+     pagecolor="#ffffff"
+     bordercolor="#666666"
+     borderopacity="1.0"
+     inkscape:pageopacity="0.0"
+     inkscape:pageshadow="2"
+     inkscape:zoom="0.98994948"
+     inkscape:cx="225.07702"
+     inkscape:cy="163.38754"
+     inkscape:document-units="mm"
+     inkscape:current-layer="layer1"
+     showgrid="false"
+     inkscape:window-width="1366"
+     inkscape:window-height="713"
+     inkscape:window-x="0"
+     inkscape:window-y="0"
+     inkscape:window-maximized="1"
+     inkscape:snap-global="false"
+     fit-margin-top="0"
+     fit-margin-left="0"
+     fit-margin-right="0"
+     fit-margin-bottom="0"
+     showguides="false" />
+  <metadata
+     id="metadata5">
+    <rdf:RDF>
+      <cc:Work
+         rdf:about="">
+        <dc:format>image/svg+xml</dc:format>
+        <dc:type
+           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
+        <dc:title></dc:title>
+      </cc:Work>
+    </rdf:RDF>
+  </metadata>
+  <g
+     inkscape:label="Layer 1"
+     inkscape:groupmode="layer"
+     id="layer1"
+     transform="translate(404.44732,-203.25495)">
+    <g
+       id="g2639"
+       transform="matrix(0.30193828,0,0,0.30193828,-282.32919,189.8183)">
+      <g
+         transform="matrix(1.0709114,0,0,1.0709114,187.29951,194.48832)"
+         id="g2113">
+        <path
+           id="path1093-6"
+           d="m -465.3209,55.882894 h 3.43136 c 2.73346,0 3.76786,-1.570286 3.76786,-1.840308 0,-0.103855 -0.0789,-0.145397 -0.18279,-0.145397 -0.23263,0 -0.76437,0.20771 -1.28364,0.577433 -0.17033,0.124626 -0.20771,0.124626 -0.60652,0.124626 h -9.03953 c -2.56729,0 -5.44614,2.214186 -5.44614,3.913252 0,0.124625 0,0.249251 0.16616,0.249251 0.49435,0 1.79877,-0.685442 1.84447,-1.183945 0.0208,-0.348953 0.0582,-0.739447 0.47357,-1.279492 0.20771,-0.253406 0.36973,-0.41542 1.2878,-0.41542 h 4.11266 c -0.34896,0.228481 -0.43619,0.540045 -0.47358,0.747755 l -2.36374,9.371864 c -0.2451,1.013623 -0.68129,2.318041 -1.27534,3.664 -0.12462,0.249252 -0.12462,0.290793 -0.12462,0.33649 0,0.120471 0.0789,0.141242 0.18278,0.141242 0.41542,0 1.42489,-0.556662 1.73646,-0.992852 0.12462,-0.20771 1.03854,-2.318041 1.47474,-4.037878 z"
+           inkscape:connector-curvature="0"
+           style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;letter-spacing:normal;word-spacing:normal;text-anchor:start;fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10.43299961;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+        <path
+           id="path1095-7"
+           d="m -454.29982,51.645615 v 0 -0.0042 0 -0.0042 0 l 0.004,-0.0042 v 0 -0.0042 -0.0042 h 0.004 v -0.0083 l 0.004,-0.0083 v -0.0083 l 0.004,-0.0083 v -0.0083 l 0.004,-0.0083 v -0.0083 l 0.004,-0.0083 0.004,-0.01246 0.004,-0.01662 0.004,-0.02077 v -0.0083 l 0.004,-0.01246 v -0.0083 l 0.004,-0.0083 v -0.0083 l 0.004,-0.0083 v -0.0042 l 0.004,-0.0083 v -0.0042 0 -0.0042 -0.0042 0 -0.0042 0 -0.0042 h 0.004 v 0 -0.0042 0 c 0,-0.07478 -0.0582,-0.20771 -0.23679,-0.20771 -0.28248,0 -1.48304,0.116317 -1.84861,0.149551 -0.11632,0.01246 -0.31987,0.02908 -0.31987,0.332335 0,0.199402 0.20355,0.199402 0.37803,0.199402 0.69375,0 0.69375,0.103855 0.69375,0.224326 0,0.09555 -0.0332,0.186939 -0.0582,0.315719 l -2.04387,8.192073 c -0.0706,0.25756 -0.0706,0.286639 -0.0706,0.319873 0,0.211863 0.17448,0.456961 0.49435,0.456961 0.16202,0 0.43204,-0.07062 0.59405,-0.373878 0.0415,-0.08308 0.17032,-0.602358 0.24095,-0.909768 l 0.33649,-1.291955 c 0.0415,-0.211864 0.19109,-0.764372 0.22848,-0.98039 0.14955,-0.552508 0.14955,-0.56497 0.43619,-1.03024 0.46527,-0.710368 1.1881,-1.532898 2.31804,-1.532898 0.81007,0 0.85576,0.668825 0.85576,1.013623 0,0.868227 -0.62313,2.471746 -0.85576,3.082413 -0.15786,0.407111 -0.21602,0.540045 -0.21602,0.785143 0,0.764372 0.63559,1.23795 1.37088,1.23795 1.45397,0 2.08956,-1.989859 2.08956,-2.210032 0,-0.191093 -0.18694,-0.191093 -0.23263,-0.191093 -0.20356,0 -0.22017,0.08724 -0.27418,0.249252 -0.33649,1.154866 -0.97208,1.75307 -1.53705,1.75307 -0.30741,0 -0.36557,-0.203555 -0.36557,-0.506812 0,-0.332335 0.0789,-0.519274 0.33649,-1.175637 0.17032,-0.448653 0.76853,-1.981551 0.76853,-2.791619 0,-0.232635 0,-0.843301 -0.54005,-1.258721 -0.2451,-0.186938 -0.66052,-0.390494 -1.3418,-0.390494 -1.05932,0 -1.8237,0.577433 -2.35959,1.200562 z"
+           inkscape:connector-curvature="0"
+           style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;letter-spacing:normal;word-spacing:normal;text-anchor:start;fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10.43299961;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+      </g>
+      <path
+         d="m -221.9664,134.53787 c 0,47.377 -38.40586,85.78279 -85.78286,85.78279 -47.37558,0 -85.78286,-38.40579 -85.78286,-85.78279 0,-47.376899 38.40728,-85.782758 85.78286,-85.782758 47.377,0 85.78286,38.405859 85.78286,85.782758 z"
+         style="fill:none;stroke:#0000ed;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path74-6-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.94732,92.899355 h 23.38378 m 0,0 v 23.383875 m 0,0 h -23.38378 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path14-5"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.94732,116.84661 h 23.38378 m 0,0 v 23.38388 m 0,0 h -23.38378 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path16-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.94732,140.79564 h 23.38378 m 0,0 v 23.38387 m 0,0 h -23.38378 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path18-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.94732,164.74466 h 23.38378 m 0,0 v 23.38388 m 0,0 h -23.38378 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path20-9"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.99839,68.950327 h 23.38391 m 0,0 v 23.383878 m 0,0 h -23.38391 m 0,0 V 68.950327"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path24-1"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.99839,92.899355 h 23.38391 m 0,0 v 23.383875 m 0,0 h -23.38391 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path26-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.99839,116.84661 h 23.38391 m 0,0 v 23.38388 m 0,0 h -23.38391 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path28-70"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.99839,140.79564 h 23.38391 m 0,0 v 23.38387 m 0,0 h -23.38391 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path30-9"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.99839,164.74466 h 23.38391 m 0,0 v 23.38388 m 0,0 h -23.38391 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path32-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.99839,188.69334 h 23.38391 m 0,0 v 23.38423 m 0,0 h -23.38391 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path34-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.87289,68.950254 h 23.38391 m 0,0 v 23.383951 m 0,0 h -23.38391 m 0,0 V 68.950254"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path36-0"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,45.001301 h 23.38388 m 0,0 v 23.383876 m 0,0 h -23.38388 m 0,0 V 45.001301"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path38-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,68.950327 h 23.38388 m 0,0 v 23.383878 m 0,0 h -23.38388 m 0,0 V 68.950327"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path40-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,92.899355 h 23.38388 m 0,0 v 23.383875 m 0,0 h -23.38388 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path42-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,116.84661 h 23.38388 m 0,0 v 23.38388 m 0,0 h -23.38388 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path44-1"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,140.79564 h 23.38388 m 0,0 v 23.38387 m 0,0 h -23.38388 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path46-8"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,164.74466 h 23.38388 m 0,0 v 23.38388 m 0,0 h -23.38388 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path48-7"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,188.69334 h 23.38388 m 0,0 v 23.38423 m 0,0 h -23.38388 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path50-9"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -356.0511,212.64088 h 23.38388 m 0,0 v 23.38395 m 0,0 h -23.38388 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path52-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,45.001301 h 23.38387 m 0,0 v 23.383876 m 0,0 h -23.38387 m 0,0 V 45.001301"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path54-0"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,68.950327 h 23.38387 m 0,0 v 23.383878 m 0,0 h -23.38387 m 0,0 V 68.950327"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path56-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,92.899355 h 23.38387 m 0,0 v 23.383875 m 0,0 h -23.38387 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path58-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,116.84661 h 23.38387 m 0,0 v 23.38388 m 0,0 h -23.38387 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path60-7"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,140.79564 h 23.38387 m 0,0 v 23.38387 m 0,0 h -23.38387 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path62-5"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,164.74466 h 23.38387 m 0,0 v 23.38388 m 0,0 h -23.38387 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path64-9"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,188.69334 h 23.38387 m 0,0 v 23.38423 m 0,0 h -23.38387 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path66-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -332.10207,212.64088 h 23.38387 m 0,0 v 23.38395 m 0,0 h -23.38387 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path68-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,45.001301 h 23.38387 m 0,0 v 23.383876 m 0,0 h -23.38387 m 0,0 V 45.001301"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path70-8"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,68.950327 h 23.38387 m 0,0 v 23.383878 m 0,0 h -23.38387 m 0,0 V 68.950327"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path72-9"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,92.899355 h 23.38387 m 0,0 v 23.383875 m 0,0 h -23.38387 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path74-7"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,116.84661 h 23.38387 m 0,0 v 23.38388 m 0,0 h -23.38387 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path76-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,140.79564 h 23.38387 m 0,0 v 23.38387 m 0,0 h -23.38387 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path78-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,164.74466 h 23.38387 m 0,0 v 23.38388 m 0,0 h -23.38387 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path80-1"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,188.69334 h 23.38387 m 0,0 v 23.38423 m 0,0 h -23.38387 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path82-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -308.15305,212.64088 h 23.38387 m 0,0 v 23.38395 m 0,0 h -23.38387 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path84-9"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,45.001301 h 23.38423 m 0,0 v 23.383876 m 0,0 h -23.38423 m 0,0 V 45.001301"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path86-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,68.950327 h 23.38423 m 0,0 v 23.383878 m 0,0 h -23.38423 m 0,0 V 68.950327"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path88-1"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,92.899355 h 23.38423 m 0,0 v 23.383875 m 0,0 h -23.38423 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path90-9"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,116.84661 h 23.38423 m 0,0 v 23.38388 m 0,0 h -23.38423 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path92-4"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,140.79564 h 23.38423 m 0,0 v 23.38387 m 0,0 h -23.38423 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path94-7"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,164.74466 h 23.38423 m 0,0 v 23.38388 m 0,0 h -23.38423 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path96-8"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,188.69334 h 23.38423 m 0,0 v 23.38423 m 0,0 h -23.38423 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path98-4"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -284.20438,212.64088 h 23.38423 m 0,0 v 23.38395 m 0,0 h -23.38423 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path100-5"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.25676,68.950327 h 23.38387 m 0,0 v 23.383878 m 0,0 h -23.38387 m 0,0 V 68.950327"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path102-0"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.25676,92.899355 h 23.38387 m 0,0 v 23.383875 m 0,0 h -23.38387 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path104-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.25676,116.84661 h 23.38387 m 0,0 v 23.38388 m 0,0 h -23.38387 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path106-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.25676,140.79564 h 23.38387 m 0,0 v 23.38387 m 0,0 h -23.38387 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path108-1"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.25676,164.74466 h 23.38387 m 0,0 v 23.38388 m 0,0 h -23.38387 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path110-0"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.25676,188.69334 h 23.38387 m 0,0 v 23.38423 m 0,0 h -23.38387 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path112-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.30774,92.899355 h 23.38388 m 0,0 v 23.383875 m 0,0 h -23.38388 m 0,0 V 92.899355"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path116-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.30774,116.84661 h 23.38388 m 0,0 v 23.38388 m 0,0 h -23.38388 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path118-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.30774,140.79564 h 23.38388 m 0,0 v 23.38387 m 0,0 h -23.38388 m 0,0 v -23.38387"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path120-0"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.30774,164.74466 h 23.38388 m 0,0 v 23.38388 m 0,0 h -23.38388 m 0,0 v -23.38388"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path122-6"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.29098,45.530468 h 23.38423 m 0,0 v 23.383876 m 0,0 h -23.38423 m 0,0 V 45.530468"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path86-3-2"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.90675,45.001301 h 23.38391 m 0,0 v 23.383951 m 0,0 h -23.38391 m 0,0 V 45.001301"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path36-0-5"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.34372,188.69334 h 23.38387 m 0,0 v 23.38423 m 0,0 h -23.38387 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path112-6-4"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -260.25712,212.60674 h 23.38423 m 0,0 v 23.38395 m 0,0 h -23.38423 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path100-5-7"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -236.34372,212.60674 h 23.38387 m 0,0 v 23.38423 m 0,0 h -23.38387 m 0,0 v -23.38423"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path112-6-4-4"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.99839,212.60674 h 23.38388 m 0,0 v 23.38395 m 0,0 h -23.38388 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path52-2-4"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.91144,189.22279 h 23.38388 m 0,0 v 23.38395 m 0,0 h -23.38388 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path52-2-3"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.91144,212.60674 h 23.38388 m 0,0 v 23.38395 m 0,0 h -23.38388 m 0,0 v -23.38395"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path52-2-4-0"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.91134,68.950327 h 23.38378 m 0,0 v 23.383877 m 0,0 h -23.38378 m 0,0 V 68.950327"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path14-5-7"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -379.96418,45.001301 h 23.38391 m 0,0 v 23.383878 m 0,0 h -23.38391 m 0,0 V 45.001301"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path24-1-8"
+         inkscape:connector-curvature="0" />
+      <path
+         d="m -403.34809,45.001301 h 23.38391 m 0,0 v 23.383878 m 0,0 h -23.38391 m 0,0 V 45.001301"
+         style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path24-1-8-6"
+         inkscape:connector-curvature="0" />
+    </g>
+  </g>
+  <g
+     inkscape:groupmode="layer"
+     id="layer2"
+     inkscape:label="Layer 2"
+     transform="translate(404.44732,-203.25495)" />
+</svg>
diff --git a/doc/doxygen/images/step-85-ghost-faces.svg b/doc/doxygen/images/step-85-ghost-faces.svg
new file mode 100644 (file)
index 0000000..396725c
--- /dev/null
@@ -0,0 +1,239 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+   xmlns:dc="http://purl.org/dc/elements/1.1/"
+   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+   xmlns:svg="http://www.w3.org/2000/svg"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   width="57.677299mm"
+   height="67.801544mm"
+   viewBox="0 0 57.677299 67.801544"
+   version="1.1"
+   id="svg8"
+   inkscape:version="0.92.5 (2060ec1f9f, 2020-04-08)"
+   sodipodi:docname="step-85-ghost-faces.svg">
+  <defs
+     id="defs2" />
+  <sodipodi:namedview
+     id="base"
+     pagecolor="#ffffff"
+     bordercolor="#666666"
+     borderopacity="1.0"
+     inkscape:pageopacity="0.0"
+     inkscape:pageshadow="2"
+     inkscape:zoom="1"
+     inkscape:cx="-169.03363"
+     inkscape:cy="95.021623"
+     inkscape:document-units="mm"
+     inkscape:current-layer="layer1"
+     showgrid="false"
+     inkscape:window-width="1366"
+     inkscape:window-height="713"
+     inkscape:window-x="0"
+     inkscape:window-y="0"
+     inkscape:window-maximized="1"
+     inkscape:snap-global="false"
+     fit-margin-top="0"
+     fit-margin-left="0"
+     fit-margin-right="0"
+     fit-margin-bottom="0"
+     showguides="false" />
+  <metadata
+     id="metadata5">
+    <rdf:RDF>
+      <cc:Work
+         rdf:about="">
+        <dc:format>image/svg+xml</dc:format>
+        <dc:type
+           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
+        <dc:title></dc:title>
+      </cc:Work>
+    </rdf:RDF>
+  </metadata>
+  <g
+     inkscape:label="Layer 1"
+     inkscape:groupmode="layer"
+     id="layer1"
+     transform="translate(263.73917,-203.17274)">
+    <g
+       id="g2730"
+       transform="matrix(0.30193828,0,0,0.30193828,-282.32919,189.8183)">
+      <g
+         transform="matrix(1.4829248,0,0,1.4829248,20.021617,115.70991)"
+         id="g993">
+        <g
+           style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;letter-spacing:normal;word-spacing:normal;text-anchor:start;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10.43299961;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+           id="g991"
+           transform="matrix(0.3,0,0,-0.3,-52.5,222.75)"
+           xml:space="preserve"
+           stroke-miterlimit="10.433"
+           font-style="normal"
+           font-variant="normal"
+           font-weight="normal"
+           font-stretch="normal"
+           font-size-adjust="none"
+           letter-spacing="normal"
+           word-spacing="normal"><path
+             style="fill:#000000;stroke-width:0"
+             inkscape:connector-curvature="0"
+             d="m 496.88,438.34 v 0.11 l -0.01,0.1 -0.02,0.09 -0.01,0.1 -0.03,0.09 -0.03,0.08 -0.03,0.09 -0.04,0.08 -0.04,0.07 -0.05,0.08 -0.05,0.07 -0.06,0.07 -0.06,0.06 -0.06,0.06 -0.07,0.06 -0.07,0.05 -0.07,0.06 -0.08,0.05 -0.16,0.09 -0.17,0.08 -0.19,0.08 -0.19,0.06 -0.2,0.06 -0.21,0.04 -0.22,0.04 -0.21,0.04 -0.23,0.03 -0.22,0.02 -0.23,0.02 -0.22,0.01 -0.23,0.01 h -0.22 l -0.22,0.01 h -0.21 l -0.21,-0.01 h -0.2 -0.2 l -0.18,-0.01 h -0.18 -0.16 -0.15 -18.78 c -2.49,0 -6.92,-1.44 -7.72,-4.14 0.09,-0.3 0.2,-0.35 0.5,-0.35 1.09,0 2.59,0.89 3.28,1.64 1.55,0 3.05,0.11 4.59,0.11 h 3.08 c 0,0 -0.2,-0.9 -0.2,-0.95 -1.94,-9.11 -4.88,-16.84 -7.06,-21.72 -0.3,-0.66 -3,-6.87 -3.69,-7.47 -2.2,0 -3.85,1.19 -4.74,3.13 -0.15,0.31 -0.15,0.5 -0.5,0.5 -1.14,0 -3.89,-1.6 -3.89,-2.44 0,-0.09 0.05,-0.2 0.11,-0.3 0.89,-2.34 3.13,-3.68 5.58,-3.68 2.53,0 5.47,1.84 7.06,3.68 1.66,1.89 4.94,10.07 6.24,13.05 h 9.61 c -0.05,-0.05 -0.1,-0.09 -0.1,-0.14 0,-0.2 0.3,-0.36 0.5,-0.36 1.1,0 3.94,1.41 3.94,2.95 0,0.35 -0.34,0.3 -0.95,0.3 h -11.96 c 1.3,4.48 2.99,8.86 3.9,13.45 h 6.57 c 1.24,0 4.74,0.25 5.63,-0.86 0.09,-0.29 -0.05,-0.79 0.09,-1.04 0.11,-0.1 0.25,-0.14 0.41,-0.14 1.19,0 3.89,1.5 3.89,2.84 z"
+             id="path987" /><path
+             style="fill:#000000;stroke-width:0"
+             inkscape:connector-curvature="0"
+             d="m 502.19,421.9 v 0.01 0 0 0.01 h 0.01 v 0.01 0.01 0 l 0.01,0.01 v 0.01 0.01 l 0.01,0.02 v 0.02 l 0.01,0.02 v 0.02 l 0.01,0.02 0.01,0.02 v 0.03 l 0.01,0.02 0.01,0.05 0.01,0.04 0.01,0.03 v 0.02 l 0.01,0.02 v 0.02 l 0.01,0.02 v 0.02 0.02 l 0.01,0.01 v 0.01 0.01 0 0.01 0 l 0.01,0.01 v 0 0.01 0 0 0.01 0 c 0,0.19 -0.14,0.5 -0.56,0.5 -0.69,0 -3.58,-0.28 -4.46,-0.36 -0.28,-0.03 -0.76,-0.06 -0.76,-0.8 0,-0.48 0.48,-0.48 0.9,-0.48 1.67,0 1.67,-0.25 1.67,-0.53 0,-0.24 -0.07,-0.46 -0.14,-0.77 l -4.92,-19.72 c -0.17,-0.62 -0.17,-0.68 -0.17,-0.76 0,-0.52 0.42,-1.11 1.19,-1.11 0.39,0 1.04,0.17 1.44,0.9 0.09,0.21 0.4,1.46 0.57,2.19 l 0.82,3.11 c 0.09,0.52 0.45,1.84 0.54,2.36 0.36,1.33 0.36,1.36 1.05,2.48 1.12,1.71 2.86,3.69 5.58,3.69 1.95,0 2.06,-1.61 2.06,-2.44 0,-2.09 -1.5,-5.95 -2.06,-7.42 -0.38,-0.98 -0.52,-1.29 -0.52,-1.89 0,-1.84 1.53,-2.98 3.3,-2.98 3.5,0 5.03,4.79 5.03,5.33 0,0.45 -0.45,0.45 -0.56,0.45 -0.49,0 -0.53,-0.2 -0.66,-0.59 -0.81,-2.79 -2.34,-4.22 -3.7,-4.22 -0.74,0 -0.88,0.48 -0.88,1.22 0,0.79 0.19,1.25 0.82,2.82 0.4,1.08 1.84,4.77 1.84,6.72 0,0.56 0,2.03 -1.3,3.03 -0.59,0.46 -1.59,0.94 -3.23,0.94 -2.55,0 -4.39,-1.39 -5.67,-2.89 z"
+             id="path989" /></g>      </g>
+      <g
+         transform="translate(-11.833664,1.0690805)"
+         id="g1466">
+        <path
+           inkscape:connector-curvature="0"
+           id="path14-8"
+           style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 96.786388,67.108982 v 23.383876 m 0,0 H 73.402605"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path16-7"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="M 96.786388,91.056597 V 114.44047"
+           sodipodi:nodetypes="cc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path18-9"
+           style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 73.402605,115.00562 h 23.383783 m 0,0 v 23.38388 m 0,0 H 73.402605"
+           sodipodi:nodetypes="cccccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path20-2"
+           style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 96.786388,138.95429 v 23.38423 m 0,0 H 73.402605"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path24-2"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 120.73545,43.159955 v 23.38423 m 0,0 H 97.350158"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path26-3"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 120.73545,67.108982 v 23.383876 m 0,0 H 97.350158"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path28-7"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 97.350158,162.90332 h 23.385292 m 0,0 v 23.38388 m 0,0 H 97.350158 m 0,0 v -23.38388" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path32-9"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 144.68271,43.159955 v 23.38423 m 0,0 h -23.38387"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path34-2"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 121.29884,186.85094 h 23.38387 m 0,0 v 23.38391 m 0,0 h -23.38387 m 0,0 v -23.38391" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path38-8"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 168.63174,43.159955 v 23.38423 m 0,0 h -23.38388"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path40-9"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 145.24786,210.79996 h 23.38388 m 0,0 v 23.38381 m -23.38388,0 v -23.38381"
+           sodipodi:nodetypes="cccccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path42-7"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 192.58077,43.159955 v 23.38423 m 0,0 h -23.38388"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path44-3"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 169.19689,210.79996 h 23.38388"
+           sodipodi:nodetypes="cc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path46-6"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 216.52979,43.159955 v 23.38423 m 0,0 H 193.1445"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path48-1"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 193.1445,186.85094 h 23.38529 m 0,0 v 23.38391 m 0,0 H 193.1445 m 0,0 v -23.38391" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path50-2"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="M 193.1445,234.18377 V 210.79996"
+           sodipodi:nodetypes="cc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path52-9"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="M 240.47705,66.544185 H 217.09318"
+           sodipodi:nodetypes="cc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path54-3"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 240.47705,67.108982 v 23.383876 m 0,0 h -23.38387 m 0,0 V 67.108982"
+           sodipodi:nodetypes="cccccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path56-1"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="m 217.09318,162.90332 h 23.38387 m 0,0 v 23.38388 m 0,0 h -23.38387 m 0,0 v -23.38388" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path60-4"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="M 264.42608,90.492858 H 241.0422"
+           sodipodi:nodetypes="cc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path62-7"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="M 264.42608,114.44047 H 241.0422 m 0,0 V 91.056597"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path64-8"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="M 264.42608,138.3895 H 241.0422 m 0,0 v -23.38388"
+           sodipodi:nodetypes="cccc" />
+        <path
+           inkscape:connector-curvature="0"
+           id="path66-4"
+           style="fill:none;stroke:#000000;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+           d="M 264.42608,162.33852 H 241.0422 m 0,0 v -23.38423"
+           sodipodi:nodetypes="cccc" />
+      </g>
+      <path
+         d="m 243.5496,133.76575 c 0,47.377 -38.40621,85.78286 -85.78286,85.78286 -47.37701,0 -85.782751,-38.40586 -85.782751,-85.78286 0,-47.376622 38.405741,-85.782832 85.782751,-85.782832 47.37665,0 85.78286,38.40621 85.78286,85.782832 z"
+         style="fill:none;stroke:#0000ff;stroke-width:1.25;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1"
+         id="path124"
+         inkscape:connector-curvature="0" />
+    </g>
+  </g>
+  <g
+     inkscape:groupmode="layer"
+     id="layer2"
+     inkscape:label="Layer 2"
+     transform="translate(263.73917,-203.17274)" />
+</svg>
diff --git a/doc/doxygen/images/step-85-level-set.png b/doc/doxygen/images/step-85-level-set.png
new file mode 100644 (file)
index 0000000..68cb21e
Binary files /dev/null and b/doc/doxygen/images/step-85-level-set.png differ
diff --git a/doc/doxygen/images/step-85-solution.png b/doc/doxygen/images/step-85-solution.png
new file mode 100644 (file)
index 0000000..3b8932a
Binary files /dev/null and b/doc/doxygen/images/step-85-solution.png differ
index ceb8ad7ec628a6ccdd3558aa609791a8d8f32db4..ca27997f697af4e031ccfada65bb71a485b9b12f 100644 (file)
   url = {https://arxiv.org/abs/2106.13877}
 }
 
+% ------------------------------------
+% Step 85
+% ------------------------------------
+
+@article{saye_2015,
+       title = {High-{Order} {Quadrature} {Methods} for {Implicitly} {Defined} {Surfaces} and {Volumes} in {Hyperrectangles}},
+       volume = {37},
+       issn = {1064-8275, 1095-7197},
+       url = {http://epubs.siam.org/doi/10.1137/140966290},
+       doi = {10.1137/140966290},
+       language = {en},
+       number = {2},
+       urldate = {2016-01-27},
+       journal = {SIAM Journal on Scientific Computing},
+       author = {Saye, R. I.},
+       month = jan,
+       year = {2015},
+       pages = {A993--A1019}
+}
+
+@article{burman_hansbo_2012,
+       title = {Fictitious domain finite element methods using cut elements: {II}. {A} stabilized {Nitsche} method},
+       volume = {62},
+       issn = {01689274},
+       url = {http://linkinghub.elsevier.com/retrieve/pii/S0168927411000298},
+       doi = {10.1016/j.apnum.2011.01.008},
+       language = {en},
+       number = {4},
+       urldate = {2015-12-17},
+       journal = {Applied Numerical Mathematics},
+       author = {Burman, Erik and Hansbo, Peter},
+       month = apr,
+       year = {2012},
+       pages = {328--341}
+}
+
+@article{cutfem_2015,
+       title = {{CutFEM}: {Discretizing} geometry and partial differential equations},
+       volume = {104},
+       issn = {00295981},
+       url = {http://doi.wiley.com/10.1002/nme.4823},
+       doi = {10.1002/nme.4823},
+       language = {en},
+       number = {7},
+       urldate = {2016-01-27},
+       journal = {International Journal for Numerical Methods in Engineering},
+       author = {Burman, Erik and Claus, Susanne and Hansbo, Peter and Larson, Mats G. and Massing, André},
+       month = nov,
+       year = {2015},
+       pages = {472--501}
+}
+
+
 % ------------------------------------
 % References used elsewhere
 % ------------------------------------
index 851bda147696c01a2e49dac9765356f6c17e6eb6..fed553831567a8e9a3788289d5cef29801fd1943 100644 (file)
  *        <td> Solving the fourth-order biharmonic equation using a lifting operator approach.
  *       </td></tr>
  *
+ *   <tr valign="top">
+ *       <td>step-85</td>
+ *       <td> Solving the Poisson equation using the cut finite element method.
+ *       </td></tr>
+ *
  * </table>
  *
  *
  *     </td>
  *     <td>
  *       step-60,
- *       step-70
+ *       step-70,
+ *       step-85
  *     </td>
  *   </tr>
  *
diff --git a/examples/step-85/CMakeLists.txt b/examples/step-85/CMakeLists.txt
new file mode 100644 (file)
index 0000000..dbc0f7d
--- /dev/null
@@ -0,0 +1,39 @@
+##
+#  CMake script for the step-85 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-85")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#    FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#    FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#    SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0)
+
+FIND_PACKAGE(deal.II 10.0.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-85/doc/builds-on b/examples/step-85/doc/builds-on
new file mode 100644 (file)
index 0000000..a226759
--- /dev/null
@@ -0,0 +1 @@
+step-12 step-46
diff --git a/examples/step-85/doc/intro.dox b/examples/step-85/doc/intro.dox
new file mode 100644 (file)
index 0000000..dbd5fec
--- /dev/null
@@ -0,0 +1,251 @@
+<i>
+This program was contributed by Simon Sticko.
+
+The material is based upon work partially supported by
+eSSENCE of e-Science and the Swedish Research Council
+under grants 2014-6088 (Kreiss) and 2017-05038 (Massing).
+</i>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+<h3>The Cut Finite Element Method</h3>
+
+In this example, we show how to use the cut finite element method (CutFEM) in deal.II.
+For illustration, we want to solve the simplest possible problem,
+so we again consider Poisson's equation:
+@f{align*}
+  -\Delta u &= f \qquad && \text{in }\, \Omega,
+  \\
+  u &= u_D       \qquad && \text{on }\, \Gamma = \partial \Omega,
+@f}
+where we choose $f(x) = 4$ and $u_D(x) = 1$.
+CutFEM is an immersed method.
+In this context,
+"immersed" means that the mesh is unfitted to the geometry of the domain, $\Omega$.
+Instead, $\Omega$ floats freely on top of a uniform background mesh, $\mathcal{T}^h$.
+@image html step-85-background-mesh.svg
+Since we no longer use the mesh to describe the geometry of the domain,
+we need some other way to represent it.
+This can be done in several ways but here we assume that $\Omega$ is described by a level set function,
+$\psi : \mathbb{R}^{\text{dim}} \to \mathbb{R}$ such that
+@f{align*}
+  \Omega &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) < 0 \}, \\
+  \Gamma &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) = 0 \}.
+@f}
+For simplicity, we choose $\Omega$ to be a unit disk, so that
+@f{equation*}
+  \psi(x) = \| x \| - 1.
+@f}
+As can be seen from the figure below,
+the level set function is negative for points in $\Omega$,
+zero on the boundary, and positive everywhere else.
+@image html step-85-level-set.png
+To solve this problem,
+we want to distribute degrees of freedom over the smallest submesh, $\mathcal{T}_\Omega^h$,
+that completely covers the domain:
+@f{equation*}
+  \mathcal{T}_\Omega^h = \{ T \in \mathcal{T}^{h} : T \cap \Omega \neq \emptyset \}.
+@f}
+This is usually referred to as the "active mesh".
+@image html step-85-active-mesh.svg
+The finite element space where we want to find our numerical solution, $u_h$, is now
+@f{equation*}
+  V_\Omega^h = \{ v \in C(\mathcal{N}_\Omega^h) : v \in Q_p(T), \, T \in \mathcal{T}_\Omega^h \},
+@f}
+where
+@f{equation*}
+  \mathcal{N}_\Omega^h = \bigcup_{T \in \mathcal{T}_\Omega^h} \overline{T},
+@f}
+and $\overline{T}$ denotes the closure of $T$.
+The set $\mathcal{N}_\Omega^h$ is sometimes referred to as the "fictitious domain".
+Since $\Omega \subset \mathcal{N}_\Omega^h$,
+we see that the numerical solution is defined over a slightly larger region than the analytical solution.
+
+In this type of immersed finite element method,
+the standard way to apply boundary conditions is using Nitsche's method.
+Multiplying the PDE with a test function, $v_h \in V_\Omega^h$,
+and integrating by parts over $\Omega$, as usual, gives us
+@f{equation*}
+  (\nabla u_h, \nabla v_h)_\Omega - (\partial_n u_h, v_h)_\Gamma = (f,v)_\Omega.
+@f}
+Let $\gamma_D > 0$ be a scalar penalty parameter and let $h$ be some measure of the local cell size.
+We now note that the following terms are consistent with the Dirichlet boundary condition:
+@f{align*}
+  -(u_h, \partial_n v_h)_\Gamma &= -(u_D, \partial_n v_h)_\Gamma, \\
+  \left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma &= \left (\frac{\gamma_D}{h}u_D, v_h \right )_\Gamma.
+@f}
+Thus, we can add these to the weak formulation to enforce the boundary condition.
+This leads to the following weak formulation:
+Find $u_h \in V_\Omega^h$ such that
+@f{equation*}
+  a_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h,
+@f}
+where
+@f{align*}
+  a_h(u_h, v_h) &=  (\nabla u_h, \nabla v_h)_\Omega
+                  - (\partial_n u_h, v_h)_\Gamma
+                  - (u_h, \partial_n v_h)_\Gamma
+                  + \left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma,
+  \\
+  L_h(v_h)      &=  (f,v)_\Omega
+                  + \left (u_D, \frac{\gamma_D}{h} v_h -\partial_n v_h \right )_\Gamma.
+@f}
+In this formulation, there is one big difference,
+compared to a standard boundary-fitted finite element method.
+On each cell,
+we need to integrate over the part of the domain and the part of the boundary that falls within the cell.
+Thus, on each cell intersected by $\Gamma$,
+we need special quadrature rules that only integrate over these parts of the cell,
+that is, over $T \cap \Omega$ and $T \cap \Gamma$.
+@image html immersed_quadratures.svg
+Since $\Omega \cap T$ is the part of the cell that lies inside the domain,
+we shall refer to the following regions
+@f{align*}
+  \{x \in T : \psi(x) < 0 \}, \\
+  \{x \in T : \psi(x) > 0 \}, \\
+  \{x \in T : \psi(x) = 0 \},
+@f}
+as the "inside", "outside" and the "surface region" of the cell $T$.
+
+The above finite element method that uses the bilinear form $a_h(\cdot, \cdot)$
+is sometimes referred to as the "naive weak formulation"
+because it suffers from the so-called "small cut problem".
+Depending on how $\Omega$ is located relative to $\mathcal{T}_h$,
+a cut between a cell, $T \in \mathcal{T}_h$, and $\Omega$ can become arbitrarily small:
+$|\Omega \cap T | \rightarrow 0$.
+For Neumann boundary conditions,
+the consequence is that the stiffness matrix can become arbitrarily ill-conditioned
+as the cut-size approaches zero.
+For a Dirichlet condition, the situation is even worse.
+For any finite choice of Nitsche constant, $\gamma_D$,
+the bilinear form $a_h(\cdot,\cdot)$ loses coercivity as the size of a cell cut approaches zero.
+This makes the above weak formulation essentially useless
+because as we refine we typically can not control how the cells intersect $\Gamma$.
+One way to avoid this problem is to add a so-called ghost penalty term, $g_h$,
+to the weak formulation (see e.g. @cite burman_hansbo_2012 and @cite cutfem_2015).
+This leads to the stabilized cut finite element method,
+which reads: Find $u_h \in V_\Omega^h$ such that
+@f{equation*}
+  A_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h,
+@f}
+where
+@f{equation*}
+  A_h(u_h,v_h) = a_h(u_h,v_h) + g_h(u_h, v_h).
+@f}
+The point of this ghost penalty is that it makes the numerical method essentially independent
+of how $\Omega$ relates to the background mesh.
+In particular, $A_h$ can be shown to be continuous and coercive,
+with constants that do not depend on how $\Omega$ intersects $\mathcal{T}^h$.
+To define the ghost penalty, let $\mathcal{T}_\Gamma^h$ be the set of intersected cells:
+@f{equation*}
+  \mathcal{T}_{\Gamma}^h = \{ T \in \mathcal{T}_{\Omega}^{h} : T \cap \Gamma \neq \emptyset \},
+@f}
+and let $\mathcal{F}_h$ denote the interior faces of the intersected cells in the active mesh:
+@f{equation*}
+  \mathcal{F}_h = \{ F = \overline{T}_+ \cap \overline{T}_- : \,
+                     T_+ \in \mathcal{T}_{\Gamma}^h, \,
+                     T_- \in \mathcal{T}_{\Omega}^h
+                  \}.
+@f}
+@image html step-85-ghost-faces.svg
+The ghost penalty acts on these faces and reads
+@f{equation*}
+  g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} g_F(u_h, v_h),
+@f}
+where $g_F$ is the face-wise ghost penalty:
+@f{equation*}
+  g_F(u_h, v_h) = \gamma_A \sum_{k=0}^p \left(\frac{h_F^{2k-1}}{k!^2}[\partial_n^k u_h], [\partial_n^k v_h] \right)_F.
+@f}
+Here, $\gamma_A$ is a penalty parameter and $h_F$ is some measure of the face size.
+We see that $g_F$ penalizes the jumps in the face-normal derivatives, $\partial_n^k$,
+over $F = \overline{T}_+ \cap \overline{T}_-$.
+Since we include all normal derivatives up to the polynomial degree,
+we weakly force the piecewise polynomial to behave as a single polynomial over $\overline{T}_+ \cup \overline{T}_-$.
+Hand-wavingly speaking,
+this is the reason why we obtain a cut-independent method when we enforce $g_F(u_h, v_h) = 0$ over the faces in $\mathcal{F}_h$.
+Here, we shall use a continuous space of $Q_1$-elements,
+so the ghost penalty is reduced to
+@f{equation*}
+  g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} (h_F [\partial_n u_h], [\partial_n v_h])_F.
+@f}
+
+<h3>Discrete Level Set Function</h3>
+A typical use case of a level set method is a problem where the domain is advected in a velocity field,
+such that the domain deforms with time.
+For such a problem,
+one would typically solve for an approximation of the level set function, $\psi_h \in V^h$,
+in a separate finite element space over the whole background mesh:
+@f{equation*}
+  V^h = \{ v \in C(\mathcal{N}^h) : v \in Q_p(T), \, T \in \mathcal{T}^h \},
+@f}
+where $\mathcal{N}^h = \bigcup_{T \in \mathcal{T}^h} \overline{T}$.
+Even if we solve a much simpler problem with a stationary domain in this tutorial,
+we shall, just to illustrate, still use a discrete level set function for the Poisson problem.
+Technically,
+this is a so-called "variational crime" because we are actually not using the bilinear form $a_h$ but instead
+@f{equation*}
+  a_h^\star(u_h, v_h) = (\nabla u_h, \nabla v_h)_{\Omega_h}
+                      - (\partial_n u_h, v_h)_{\Gamma_h} + \ldots
+@f}
+This is an approximation of $a_h$ since we integrate over the approximations of the geometry that we get via the discrete level set function:
+@f{align*}
+  \Omega_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) < 0 \}, \\
+  \Gamma_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) = 0 \}.
+@f}
+Using $\Omega_h$ instead of $\Omega$ in the method will give rise to a larger error in the numerical solution.
+This is often referred to as the "geometrical error".
+However, when the same element order, $p$, is used in $V^h$ and $V_\Omega^h$,
+one can often show that the method gives the same order of convergence
+as if the exact domain would have been used.
+Furthermore, deal.II allows us to independently choose a more accurate geometry representation
+with a higher-order level set function, compared to the function space for solving the Poisson equation.
+
+<h3>The MeshClassifier Class</h3>
+Even if we have used $\mathcal{T}_\Omega^h$ to define the finite element space,
+we will not create this submesh in practice. As in step-46, we shall
+instead use the hp-framework. To create $V_\Omega^h$, we first add an FE_Q and an
+FE_Nothing element to an hp::FECollection. We then iterate over each cell,
+$T$, and depending on whether $T$ belongs to $\mathcal{T}_\Omega^h$ or not,
+we set the active_fe_index to either 0 or 1. To do so, we need to
+determine if a given cell is in $\mathcal{T}_\Omega^h$ or not.
+For this purpose, we will use the class NonMatching::MeshClassifier.
+The MeshClassifier takes the discrete level set function,
+described as a (DoFHandler, Vector)-pair, as arguments to its constructor:
+@code
+  MeshClassifier(const DoFHandler<dim> &level_set_dof_handler,
+                 const VectorType &     level_set);
+@endcode
+When we call the reclassify() function on an object of this class,
+each active cell and face is associated with one of the values
+{inside, outside, intersected} of the enum NonMatching::LocationToLevelSet.
+Here, "inside" means that the level set function is negative over the whole cell
+so that it lies completely inside the domain.
+Analogously, "outside" means that $\psi$ is positive over the whole cell,
+and "intersected" means that $\psi$ varies in sign over $T$
+so that the zero-contour of $\psi$ goes through $T$.
+
+| LocationToLevelSet | $\psi(x)$ for $x \in T$  | Relation to $\Omega$           |
+|:------------------:|:------------------------:|:------------------------------:|
+| inside             | $\psi(x) < 0$            | $T \cap \Omega = T$            |
+| outside            | $0 < \psi(x)$            | $T \cap \Omega = \emptyset$    |
+| intersected        | $\psi(x)$ varies in sign | $T \cap \Gamma \neq \emptyset$ |
+
+Each active face is classified in the same way, according to how the sign of $\psi$ varies over the face.
+MeshClassifier lets you query this information for a given cell/face via its MeshClassifier::location_to_level_set() methods:
+@code
+  NonMatching::MeshClassifier<dim> mesh_classifier(dof_handler, level_set);
+  mesh_classifier.reclassify();
+
+  for (const auto &cell : triangulation.active_cell_iterators())
+    {
+      NonMatching::LocationToLevelSet cell_location =
+        mesh_classifier.location_to_level_set(cell);
+
+      for (const unsigned int f : cell->face_indices())
+        {
+          NonMatching::LocationToLevelSet face_location =
+            mesh_classifier.location_to_level_set(cell, f);
+        }
+    }
+@endcode
diff --git a/examples/step-85/doc/kind b/examples/step-85/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-85/doc/results.dox b/examples/step-85/doc/results.dox
new file mode 100644 (file)
index 0000000..b8cfe75
--- /dev/null
@@ -0,0 +1,20 @@
+<h1>Results</h1>
+
+The numerical solution for one of the refinements is shown in the below figure.
+The zero-contour of the level set function is shown as a white line.
+On the intersected cells,
+we see that the numerical solution has a value also outside $\overline{\Omega}$.
+As mentioned earlier, this extension of the solution is artificial.
+
+The results of the convergence study is shown in the table below.
+We see that the $L^2$ error decreases as we refine and that the estimated
+order of convergence, EOC, is close to 2.
+
+@image html step-85-solution.png
+
+| Cycle | Mesh size | $L^2$-Error | EOC  |
+|:-----:|:---------:|:-----------:|:----:|
+|   0   |  0.3025   |  8.0657e-02 |  -   |
+|   1   |  0.1513   |  1.8711e-02 | 2.11 |
+|   2   |  0.0756   |  4.1624e-03 | 2.17 |
+|   3   |  0.0378   |  9.3979e-04 | 2.15 |
diff --git a/examples/step-85/doc/tooltip b/examples/step-85/doc/tooltip
new file mode 100644 (file)
index 0000000..927b4dc
--- /dev/null
@@ -0,0 +1 @@
+Solving Poisson's equation using the cut finite element method.
diff --git a/examples/step-85/step-85.cc b/examples/step-85/step-85.cc
new file mode 100644 (file)
index 0000000..a3a8e8b
--- /dev/null
@@ -0,0 +1,714 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2021 - 2021 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// @sect3{Include files}
+
+// The first include files have all been treated in previous examples.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <vector>
+
+// The first new header contains some common level set functions.
+// For example, the spherical geometry that we use here.
+#include <deal.II/base/function_signed_distance.h>
+
+// We also need 3 new headers from the NonMatching namespace.
+#include <deal.II/non_matching/fe_immersed_values.h>
+#include <deal.II/non_matching/fe_values.h>
+#include <deal.II/non_matching/mesh_classifier.h>
+
+// @sect3{The LaplaceSolver class Template}
+// We then define the main class that solves the Laplace problem.
+
+namespace Step85
+{
+  using namespace dealii;
+
+  template <int dim>
+  class LaplaceSolver
+  {
+  public:
+    LaplaceSolver();
+
+    void run();
+
+  private:
+    void make_grid();
+
+    void setup_discrete_level_set();
+
+    void distribute_dofs();
+
+    void initialize_matrices();
+
+    void assemble_system();
+
+    void solve();
+
+    void output_results() const;
+
+    double compute_L2_error() const;
+
+    bool face_has_ghost_penalty(
+      const typename Triangulation<dim>::active_cell_iterator &cell,
+      const unsigned int face_index) const;
+
+    const unsigned int fe_degree;
+
+    const Functions::ConstantFunction<dim> rhs_function;
+    const Functions::ConstantFunction<dim> boundary_condition;
+
+    Triangulation<dim> triangulation;
+
+    // We need two separate DoFHandlers. The first manages the DoFs for the
+    // discrete level set function that describes the geometry of the domain.
+    const FE_Q<dim> fe_level_set;
+    DoFHandler<dim> level_set_dof_handler;
+    Vector<double>  level_set;
+
+    // The second DoFHandler manages the DoFs for the solution of the Poisson
+    // equation.
+    hp::FECollection<dim> fe_collection;
+    DoFHandler<dim>       dof_handler;
+    Vector<double>        solution;
+
+    NonMatching::MeshClassifier<dim> mesh_classifier;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> stiffness_matrix;
+    Vector<double>       rhs;
+  };
+
+
+
+  template <int dim>
+  LaplaceSolver<dim>::LaplaceSolver()
+    : fe_degree(1)
+    , rhs_function(4.0)
+    , boundary_condition(1.0)
+    , fe_level_set(fe_degree)
+    , level_set_dof_handler(triangulation)
+    , dof_handler(triangulation)
+    , mesh_classifier(level_set_dof_handler, level_set)
+  {}
+
+
+
+  // @sect3{Setting up the Background Mesh}
+  // We generate a background mesh with perfectly Cartesian cells. Our domain is
+  // a unit disc centered at the origin, so we need to make the background mesh
+  // a bit larger than $[-1, 1]^{\text{dim}}$ to completely cover $\Omega$.
+  template <int dim>
+  void LaplaceSolver<dim>::make_grid()
+  {
+    std::cout << "Creating background mesh" << std::endl;
+
+    GridGenerator::hyper_cube(triangulation, -1.21, 1.21);
+    triangulation.refine_global(2);
+  }
+
+
+
+  // @sect3{Setting up the Discrete Level Set Function}
+  // The discrete level set function is defined on the whole background mesh.
+  // Thus, to set up the DoFHandler for the level set function, we distribute
+  // DoFs over all elements in $\mathcal{T}_h$. We then set up the discrete
+  // level set function by interpolating onto this finite element space.
+  template <int dim>
+  void LaplaceSolver<dim>::setup_discrete_level_set()
+  {
+    std::cout << "Setting up discrete level set function" << std::endl;
+
+    level_set_dof_handler.distribute_dofs(fe_level_set);
+    level_set.reinit(level_set_dof_handler.n_dofs());
+
+    const Functions::SignedDistance::Sphere<dim> signed_distance_sphere;
+    VectorTools::interpolate(level_set_dof_handler,
+                             signed_distance_sphere,
+                             level_set);
+  }
+
+
+
+  // @sect3{Setting up the Finite Element Space}
+  // To set up the finite element space $V_\Omega^h$, we will use 2 different
+  // elements: FE_Q and FE_Nothing. For better readability we define an enum for
+  // the indices in the order we store them in the hp::FECollection.
+  enum ActiveFEIndex
+  {
+    lagrange = 0,
+    nothing  = 1
+  };
+
+  // We then use the MeshClassifier to check LocationToLevelSet for each cell in
+  // the mesh and tell the DoFHandler to use FE_Q on elements that are inside or
+  // intersected, and FE_Nothing on the elements that are outside.
+  template <int dim>
+  void LaplaceSolver<dim>::distribute_dofs()
+  {
+    std::cout << "Distributing degrees of freedom" << std::endl;
+
+    fe_collection.push_back(FE_Q<dim>(fe_degree));
+    fe_collection.push_back(FE_Nothing<dim>());
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        const NonMatching::LocationToLevelSet cell_location =
+          mesh_classifier.location_to_level_set(cell);
+
+        if (cell_location == NonMatching::LocationToLevelSet::outside)
+          cell->set_active_fe_index(ActiveFEIndex::nothing);
+        else
+          cell->set_active_fe_index(ActiveFEIndex::lagrange);
+      }
+
+    dof_handler.distribute_dofs(fe_collection);
+  }
+
+
+
+  // @sect3{Sparsity Pattern}
+  // The added ghost penalty results in a sparsity pattern similar to a DG
+  // method with a symmetric-interior-penalty term. Thus, we can use the
+  // make_flux_sparsity_pattern() function to create it. However, since the
+  // ghost-penalty terms only act on the faces in $\mathcal{F}_h$, we can pass
+  // in a lambda function that tells make_flux_sparsity_pattern() over which
+  // faces the flux-terms appear. This gives us a sparsity pattern with minimal
+  // number of entries. When passing a lambda function,
+  // make_flux_sparsity_pattern requires us to also pass cell and face coupling
+  // tables to it. If the problem was vector-valued, these tables would allow us
+  // to couple only some of the vector components. This is discussed in step-46.
+  template <int dim>
+  void LaplaceSolver<dim>::initialize_matrices()
+  {
+    std::cout << "Initializing matrices" << std::endl;
+
+    const auto face_has_flux_coupling = [&](const auto &       cell,
+                                            const unsigned int face_index) {
+      return this->face_has_ghost_penalty(cell, face_index);
+    };
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+
+    const unsigned int           n_components = fe_collection.n_components();
+    Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components);
+    Table<2, DoFTools::Coupling> face_coupling(n_components, n_components);
+    cell_coupling[0][0] = DoFTools::always;
+    face_coupling[0][0] = DoFTools::always;
+
+    const AffineConstraints<double> constraints;
+    const bool                      keep_constrained_dofs = true;
+
+    DoFTools::make_flux_sparsity_pattern(dof_handler,
+                                         dsp,
+                                         constraints,
+                                         keep_constrained_dofs,
+                                         cell_coupling,
+                                         face_coupling,
+                                         numbers::invalid_subdomain_id,
+                                         face_has_flux_coupling);
+    sparsity_pattern.copy_from(dsp);
+
+    stiffness_matrix.reinit(sparsity_pattern);
+    solution.reinit(dof_handler.n_dofs());
+    rhs.reinit(dof_handler.n_dofs());
+  }
+
+
+
+  // The following function describes which faces are part of the set
+  // $\mathcal{F}_h$. That is, it returns true if the face of the incoming cell
+  // belongs to the set $\mathcal{F}_h$.
+  template <int dim>
+  bool LaplaceSolver<dim>::face_has_ghost_penalty(
+    const typename Triangulation<dim>::active_cell_iterator &cell,
+    const unsigned int                                       face_index) const
+  {
+    if (cell->at_boundary(face_index))
+      return false;
+
+    const NonMatching::LocationToLevelSet cell_location =
+      mesh_classifier.location_to_level_set(cell);
+
+    const NonMatching::LocationToLevelSet neighbor_location =
+      mesh_classifier.location_to_level_set(cell->neighbor(face_index));
+
+    if (cell_location == NonMatching::LocationToLevelSet::intersected &&
+        neighbor_location != NonMatching::LocationToLevelSet::outside)
+      return true;
+
+    if (neighbor_location == NonMatching::LocationToLevelSet::intersected &&
+        cell_location != NonMatching::LocationToLevelSet::outside)
+      return true;
+
+    return false;
+  }
+
+
+
+  // @sect3{Assembling the System}
+  template <int dim>
+  void LaplaceSolver<dim>::assemble_system()
+  {
+    std::cout << "Assembling" << std::endl;
+
+    const unsigned int n_dofs_per_cell = fe_collection[0].dofs_per_cell;
+    FullMatrix<double> local_stiffness(n_dofs_per_cell, n_dofs_per_cell);
+    Vector<double>     local_rhs(n_dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(n_dofs_per_cell);
+
+    const double ghost_parameter   = 0.5;
+    const double nitsche_parameter = 5 * (fe_degree + 1) * fe_degree;
+
+    // Since the ghost penalty is similar to a DG flux term, the simplest way to
+    // assemble it is to use an FEInterfaceValues object.
+    const QGauss<dim - 1>  face_quadrature(fe_degree + 1);
+    FEInterfaceValues<dim> fe_interface_values(fe_collection[0],
+                                               face_quadrature,
+                                               update_gradients |
+                                                 update_JxW_values |
+                                                 update_normal_vectors);
+
+
+    // As we iterate over the cells in the mesh, we would in principle have to
+    // do the following on each cell, $T$,
+    //
+    // 1. Construct one quadrature rule to integrate over the intersection with
+    // the domain, $T \cap \Omega$, and one quadrature rule to integrate over
+    // the intersection with the boundary, $T \cap \Gamma$.
+    // 2. Create FEValues-like objects with the new quadratures.
+    // 3. Assemble the local matrix using the created FEValues-objects.
+    //
+    // To make the assembly easier, we use the class NonMatching::FEValues,
+    // which does the above steps 1 and 2 for us. The algorithm @cite saye_2015
+    // that is used to generate the quadrature rules on the intersected cells
+    // uses a 1-dimensional quadrature rule as base. Thus, we pass a 1D
+    // Gauss--Legendre quadrature to the constructor of NonMatching::FEValues.
+    // On the non-intersected cells, a tensor product of this 1D-quadrature will
+    // be used.
+    //
+    // As stated in the introduction, each cell has 3 different regions: inside,
+    // surface, and outside, where the level set function in each region is
+    // negative, zero, and positive. We need an UpdateFlags variable for each
+    // such region. These are stored on an object of type
+    // NonMatching::RegionUpdateFlags, which we pass to NonMatching::FEValues.
+    const QGauss<1> quadrature_1D(fe_degree + 1);
+
+    NonMatching::RegionUpdateFlags region_update_flags;
+    region_update_flags.inside = update_values | update_gradients |
+                                 update_JxW_values | update_quadrature_points;
+    region_update_flags.surface = update_values | update_gradients |
+                                  update_JxW_values | update_quadrature_points |
+                                  update_normal_vectors;
+
+    NonMatching::FEValues<dim> non_matching_fe_values(fe_collection,
+                                                      quadrature_1D,
+                                                      region_update_flags,
+                                                      mesh_classifier,
+                                                      level_set_dof_handler,
+                                                      level_set);
+
+    // As we iterate over the cells, we don't need to do anything on the cells
+    // that have FENothing elements. To disregard them we use an iterator
+    // filter.
+    for (const auto &cell :
+         dof_handler.active_cell_iterators() |
+           IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange))
+      {
+        local_stiffness = 0;
+        local_rhs       = 0;
+
+        const double cell_side_length = cell->minimum_vertex_distance();
+
+        // First, we call the reinit function of our NonMatching::FEValues
+        // object. In the background, NonMatching::FEValues uses the
+        // MeshClassifier passed to its constructor to check if the incoming
+        // cell is intersected. If that is the case, NonMatching::FEValues calls
+        // the NonMatching::QuadratureGenerator in the background to create the
+        // immersed quadrature rules.
+        non_matching_fe_values.reinit(cell);
+
+        // After calling reinit, we can retrieve a dealii::FEValues object with
+        // quadrature points that corresponds to integrating over the inside
+        // region of the cell. This is the object we use to do the local
+        // assembly. This is similar to how hp::FEValues builds dealii::FEValues
+        // objects. However, one difference here is that the dealii::FEValues
+        // object is returned as an optional. This is a type that wraps an
+        // object that may or may not be present. This requires us to add an
+        // if-statement to check if the returned optional contains a value,
+        // before we use it. This might seem odd at first. Why does the function
+        // not just return a reference to a const FEValues<dim>? The reason is
+        // that in an immersed method, we have essentially no control of how the
+        // cuts occur. Even if the cell is formally intersected: $T \cap \Omega
+        // \neq \emptyset$, it might be that the cut is only of floating point
+        // size $|T \cap \Omega| \sim \epsilon$. When this is the case, we can
+        // not expect that the algorithm that generates the quadrature rule
+        // produces anything useful. It can happen that the algorithm produces 0
+        // quadrature points. When this happens, the returned optional will not
+        // contain a value, even if the cell is formally intersected.
+        const std_cxx17::optional<FEValues<dim>> &inside_fe_values =
+          non_matching_fe_values.get_inside_fe_values();
+
+        if (inside_fe_values)
+          for (const unsigned int q :
+               inside_fe_values->quadrature_point_indices())
+            {
+              const Point<dim> &point = inside_fe_values->quadrature_point(q);
+              for (const unsigned int i : inside_fe_values->dof_indices())
+                {
+                  for (const unsigned int j : inside_fe_values->dof_indices())
+                    {
+                      local_stiffness(i, j) +=
+                        inside_fe_values->shape_grad(i, q) *
+                        inside_fe_values->shape_grad(j, q) *
+                        inside_fe_values->JxW(q);
+                    }
+                  local_rhs(i) += rhs_function.value(point) *
+                                  inside_fe_values->shape_value(i, q) *
+                                  inside_fe_values->JxW(q);
+                }
+            }
+
+        // In the same way, we can use NonMatching::FEValues to retrieve an
+        // FEFaceValues-like object to integrate over $T \cap \Gamma$. The only
+        // thing that is new here is the type of the object. The transformation
+        // from quadrature weights to JxW-values is different for surfaces, so
+        // we need a new class: NonMatching::FEImmersedSurfaceValues. In
+        // addition to the ordinary functions shape_value(..), shape_grad(..),
+        // etc., one can use its normal_vector(..)-function to get an outward
+        // normal to the immersed surface, $\Gamma$. In terms of the level set
+        // function, this normal reads
+        // @f{equation*}
+        //   n = \frac{\nabla \psi}{\| \nabla \psi \|}.
+        // @f}
+        // An additional benefit of std::optional is that we do not need any
+        // other check for whether we are on intersected cells: In case we are
+        // on an inside cell, we get an empty object here.
+        const std_cxx17::optional<NonMatching::FEImmersedSurfaceValues<dim>>
+          &surface_fe_values = non_matching_fe_values.get_surface_fe_values();
+
+        if (surface_fe_values)
+          {
+            for (const unsigned int q :
+                 surface_fe_values->quadrature_point_indices())
+              {
+                const Point<dim> &point =
+                  surface_fe_values->quadrature_point(q);
+                const Tensor<1, dim> &normal =
+                  surface_fe_values->normal_vector(q);
+                for (const unsigned int i : surface_fe_values->dof_indices())
+                  {
+                    for (const unsigned int j :
+                         surface_fe_values->dof_indices())
+                      {
+                        local_stiffness(i, j) +=
+                          (-normal * surface_fe_values->shape_grad(i, q) *
+                             surface_fe_values->shape_value(j, q) +
+                           -normal * surface_fe_values->shape_grad(j, q) *
+                             surface_fe_values->shape_value(i, q) +
+                           nitsche_parameter / cell_side_length *
+                             surface_fe_values->shape_value(i, q) *
+                             surface_fe_values->shape_value(j, q)) *
+                          surface_fe_values->JxW(q);
+                      }
+                    local_rhs(i) +=
+                      boundary_condition.value(point) *
+                      (nitsche_parameter / cell_side_length *
+                         surface_fe_values->shape_value(i, q) -
+                       normal * surface_fe_values->shape_grad(i, q)) *
+                      surface_fe_values->JxW(q);
+                  }
+              }
+          }
+
+        cell->get_dof_indices(local_dof_indices);
+
+        stiffness_matrix.add(local_dof_indices, local_stiffness);
+        rhs.add(local_dof_indices, local_rhs);
+
+        // The assembly of the ghost penalty term is straight forward. As we
+        // iterate over the local faces, we first check if the current face
+        // belongs to the set $\mathcal{F}_h$. The actual assembly is simple
+        // using FEInterfaceValues. Assembling in this we will traverse each
+        // internal face in the mesh twice, so in order to get the penalty
+        // constant we expect, we multiply the penalty term with a factor 1/2.
+        for (unsigned int f : cell->face_indices())
+          if (face_has_ghost_penalty(cell, f))
+            {
+              const unsigned int invalid_subface =
+                numbers::invalid_unsigned_int;
+
+              fe_interface_values.reinit(cell,
+                                         f,
+                                         invalid_subface,
+                                         cell->neighbor(f),
+                                         cell->neighbor_of_neighbor(f),
+                                         invalid_subface);
+
+              const unsigned int n_interface_dofs =
+                fe_interface_values.n_current_interface_dofs();
+              FullMatrix<double> local_stabilization(n_interface_dofs,
+                                                     n_interface_dofs);
+              for (unsigned int q = 0;
+                   q < fe_interface_values.n_quadrature_points;
+                   ++q)
+                {
+                  const Tensor<1, dim> normal = fe_interface_values.normal(q);
+                  for (unsigned int i = 0; i < n_interface_dofs; ++i)
+                    for (unsigned int j = 0; j < n_interface_dofs; ++j)
+                      {
+                        local_stabilization(i, j) +=
+                          .5 * ghost_parameter * cell_side_length * normal *
+                          fe_interface_values.jump_in_shape_gradients(i, q) *
+                          normal *
+                          fe_interface_values.jump_in_shape_gradients(j, q) *
+                          fe_interface_values.JxW(q);
+                      }
+                }
+
+              const std::vector<types::global_dof_index>
+                local_interface_dof_indices =
+                  fe_interface_values.get_interface_dof_indices();
+
+              stiffness_matrix.add(local_interface_dof_indices,
+                                   local_stabilization);
+            }
+      }
+  }
+
+
+  // @sect3{Solving the System}
+  template <int dim>
+  void LaplaceSolver<dim>::solve()
+  {
+    std::cout << "Solving system" << std::endl;
+
+    const unsigned int max_iterations = solution.size();
+    SolverControl      solver_control(max_iterations);
+    SolverCG<>         solver(solver_control);
+    solver.solve(stiffness_matrix, solution, rhs, PreconditionIdentity());
+  }
+
+
+
+  // @sect3{Data Output}
+  // Since both DoFHandler instances use the same triangulation, we can add both
+  // the level set function and the solution to the same vtu-file. Further, we
+  // do not want to output the cells that have LocationToLevelSet value outside.
+  // To disregard them, we write a small lambda function and use the
+  // set_cell_selection function of the DataOut class.
+  template <int dim>
+  void LaplaceSolver<dim>::output_results() const
+  {
+    std::cout << "Writing vtu file" << std::endl;
+
+    DataOut<dim> data_out;
+    data_out.add_data_vector(dof_handler, solution, "solution");
+    data_out.add_data_vector(level_set_dof_handler, level_set, "level_set");
+
+    data_out.set_cell_selection(
+      [this](const typename Triangulation<dim>::cell_iterator &cell) {
+        return cell->is_active() &&
+               mesh_classifier.location_to_level_set(cell) !=
+                 NonMatching::LocationToLevelSet::outside;
+      });
+
+    data_out.build_patches();
+    std::ofstream output("step-85.vtu");
+    data_out.write_vtu(output);
+  }
+
+
+
+  // @sect3{$L^2$-Error}
+  // To test that the implementation works as expected, we want to compute the
+  // error in the solution in the $L^2$-norm. The analytical solution to the
+  // Poisson problem stated in the introduction reads
+  // @f{align*}
+  //  u(x) = 1 - \frac{2}{\text{dim}}(\| x \|^2 - 1) , \qquad x \in
+  //  \overline{\Omega}.
+  // @f}
+  // We first create a function corresponding to the analytical solution:
+  template <int dim>
+  class AnalyticalSolution : public Function<dim>
+  {
+  public:
+    double value(const Point<dim> & point,
+                 const unsigned int component = 0) const override;
+  };
+
+
+
+  template <int dim>
+  double AnalyticalSolution<dim>::value(const Point<dim> & point,
+                                        const unsigned int component) const
+  {
+    AssertIndexRange(component, this->n_components);
+    (void)component;
+
+    return 1. - 2. / dim * (point.norm_square() - 1.);
+  }
+
+
+
+  // Of course, the analytical solution, and thus also the error, is only
+  // defined in $\overline{\Omega}$. Thus, to compute the $L^2$-error we must
+  // proceed in the same way as when we assembled the linear system. We first
+  // create an NonMatching::FEValues object.
+  template <int dim>
+  double LaplaceSolver<dim>::compute_L2_error() const
+  {
+    std::cout << "Computing L2 error" << std::endl;
+
+    const QGauss<1> quadrature_1D(fe_degree + 1);
+
+    NonMatching::RegionUpdateFlags region_update_flags;
+    region_update_flags.inside =
+      update_values | update_JxW_values | update_quadrature_points;
+
+    NonMatching::FEValues<dim> non_matching_fe_values(fe_collection,
+                                                      quadrature_1D,
+                                                      region_update_flags,
+                                                      mesh_classifier,
+                                                      level_set_dof_handler,
+                                                      level_set);
+
+    // We then iterate iterate over the cells that have LocationToLevelSetValue
+    // value inside or intersected again. For each quadrature point, we compute
+    // the pointwise error and use this to compute the integral.
+    const AnalyticalSolution<dim> analytical_solution;
+    double                        error_L2_squared = 0;
+
+    for (const auto &cell :
+         dof_handler.active_cell_iterators() |
+           IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange))
+      {
+        non_matching_fe_values.reinit(cell);
+
+        const std_cxx17::optional<FEValues<dim>> &fe_values =
+          non_matching_fe_values.get_inside_fe_values();
+
+        if (fe_values)
+          {
+            std::vector<double> solution_values(fe_values->n_quadrature_points);
+            fe_values->get_function_values(solution, solution_values);
+
+            for (const unsigned int q : fe_values->quadrature_point_indices())
+              {
+                const Point<dim> &point = fe_values->quadrature_point(q);
+                const double      error_at_point =
+                  solution_values.at(q) - analytical_solution.value(point);
+                error_L2_squared +=
+                  std::pow(error_at_point, 2) * fe_values->JxW(q);
+              }
+          }
+      }
+
+    return std::sqrt(error_L2_squared);
+  }
+
+
+
+  // @sect3{A Convergence Study}
+  // Finally, we do a convergence study to check that the $L^2$-error decreases
+  // with the expected rate. We refine the background mesh a few times. In each
+  // refinement cycle, we solve the problem, compute the error, and add the
+  // $L^2$-error and the mesh size to a ConvergenceTable.
+  template <int dim>
+  void LaplaceSolver<dim>::run()
+  {
+    ConvergenceTable   convergence_table;
+    const unsigned int n_refinements = 3;
+
+    make_grid();
+    for (unsigned int cycle = 0; cycle <= n_refinements; cycle++)
+      {
+        std::cout << "Refinement cycle " << cycle << std::endl;
+        triangulation.refine_global(1);
+        setup_discrete_level_set();
+        std::cout << "Classifying cells" << std::endl;
+        mesh_classifier.reclassify();
+        distribute_dofs();
+        initialize_matrices();
+        assemble_system();
+        solve();
+        if (cycle == 1)
+          output_results();
+        const double error_L2 = compute_L2_error();
+        const double cell_side_length =
+          triangulation.begin_active()->minimum_vertex_distance();
+
+        convergence_table.add_value("Cycle", cycle);
+        convergence_table.add_value("Mesh size", cell_side_length);
+        convergence_table.add_value("L2-Error", error_L2);
+
+        convergence_table.evaluate_convergence_rates(
+          "L2-Error", ConvergenceTable::reduction_rate_log2);
+        convergence_table.set_scientific("L2-Error", true);
+
+        std::cout << std::endl;
+        convergence_table.write_text(std::cout);
+        std::cout << std::endl;
+      }
+  }
+
+} // namespace Step85
+
+
+
+// @sect3{The main() function}
+int main()
+{
+  const int dim = 2;
+
+  Step85::LaplaceSolver<dim> laplace_solver;
+  laplace_solver.run();
+}

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.